Метод ланцоша собственные значения

Алгоритм Ланцоша

В линейной алгебре , то Ланцош алгоритм (или Ланцош метод ) представляет собой итерационный алгоритм для определения значений и собственных векторов матрицы А квадратной матрицы , или сингулярного разложения прямоугольной матрицы.

Этот алгоритм не имеет никакой связи с Lanczos оконного (используемым, например , для изменения размера изображения), за исключение того, что оба получают свои имена из того же самого изобретателя, венгерского физика и математика Ланцош .

Резюме

• 1 Описание алгоритма
  • 1.1 Метод повторной мощности
  • 1.2 Алгоритм Ланцоша
  • 4.1 Связанные статьи
  • 4.2 Ссылки
  • 4.3 Внешние ссылки

  Описание алгоритма

  Итерированный метод мощности

  Метод повторной мощности для нахождения наибольшего собственного значения симметричной квадратной матрицы порядка p выглядит следующим образом. В

  Если — вектор и , то . Икс 0 > Икс нет + 1 знак равно В Икс нет = Ax_ > Икс нет знак равно В нет Икс 0 <\ displaystyle x_ = A ^ x_ >

  Мы знаем, что симметричная матрица диагонализуема, а члены диагональной матрицы являются собственными значениями матрицы. Так что либо σ я , я знак равно 1 . п , я = 1 . p>

  В знак равно U ⋅ ( σ 1 0 ⋱ σ я ⋱ ⋮ ⋱ 0 σ п ) ⋅ U ′ знак равно U ⋅ диагональ ⁡ ( σ я ) ⋅ U ′ \ sigma _ &&&&& 0 \\ & \ ddots &&&& \\ && \ sigma _ &&& \\ &&&& \ ddots && \\ && \ vdots && \ ddots & \\ 0 &&&&& \ sigma _

  \\\ end > \ cdot U ‘= U \ cdot \ operatorname (\ sigma _ ) \ cdot U’>

  собственное разложение из , где U представляет собой ортогональную матрицу ; так В

  В нет знак равно U диагональ ⁡ ( σ я нет ) U ′ = U \ operatorname (\ sigma _ ^ ) U ‘>

  Для очень больших значений в диагональной матрице собственных значений будет преобладать наибольшее из них — за исключением случая, когда есть два больших равных значения: действительно, пусть будет собственным значением большего модуля. нет σ м в Икс >

  Lim нет → ∞ В нет знак равно Lim нет → ∞ | σ м в Икс | нет ⋅ U диагональ ⁡ ( σ я нет | σ м в Икс | нет ) U ′ знак равно U диагональ ⁡ ( 0 . 0 , ( σ м в Икс | σ м в Икс | ) нет , 0 . 0 ) U ′ A ^ = \ lim _ | \ sigma _ | ^ \ cdot U \ operatorname ( ^ > <| \ sigma _ | ^ >>) U ‘= U \ operatorname (0 . 0, \ left ( <\ frac <\ sigma _ > <| \ sigma _ |>> \ right) ^ , 0 . 0) U ‘>

  Затем сходится к наибольшему собственному значению и соответствующему собственному вектору. | Икс нет + 1 | / | Икс нет | | / | x_ |> Икс нет / | Икс нет | <\ displaystyle x_ / | x_ |>

  Если есть несколько максимальных собственных значений, сходится к вектору в подпространстве, созданном собственными векторами, связанными с этими значениями. После определения первого можно последовательно ограничить алгоритм ядром известных собственных векторов, чтобы определить остальные. Икс нет >

  На практике у этого алгоритма есть два недостатка: уязвимость к ошибкам округления и иногда слишком низкая скорость сходимости.

  Алгоритм Ланцоша

  Алгоритм Ланцоша улучшает предыдущий метод, в котором каждое из них ограничено ортогональностью всех предыдущих значений. Во время построения этих векторов нормировочные константы группируются в трехдиагональную матрицу , наиболее значимые собственные значения которой быстро сходятся к собственным значениям исходной матрицы. ты

  Тогда умножение на — единственная крупномасштабная операция, делающая этот метод интересным. В

  Приложения, варианты и обобщения

  Алгоритм Ланцоша особенно полезен для определения разложения очень больших матриц в метеорологии или при обработке естественного языка, где эти матрицы могут содержать сотни тысяч терминов.

  Это обобщает алгоритм Арнольди (en) для несимметричных матриц.

  Вариант этого алгоритма используется для разрешения линейных систем (в частности, для полых линейных систем) и поиска элементов ядра матрицы. В этой форме он часто используется в рамках таких алгоритмов, как квадратное решето или алгебраическое решето для факторизации целых чисел или индексный метод для задачи дискретного логарифмирования . Это связано с методом сопряженных градиентов .

  Примечания и ссылки

  1. ↑ Для демонстрации ср. П.Г. Чиарле, Введение в численный матричный анализ и оптимизацию , под ред. Массон, сб. «Математика. Прил. на степень магистра », 1985 г. ( перепечатка 2001 г.) ( ISBN2-225-68893-1 ) , гл. 1, теорема 1.2.1
  2. ↑ (in) Р. Беллман, Введение в матричный анализ , SIAM al. «Классика в приложении. Математика. «, 1960 г. ( перепечатка 1970, 1997) ( ISBN0-89871-399-4 ) , «3» , стр. 35–36

  Смотрите также

  Статьи по Теме

  • Алгоритм поиска собственных значений
  • Метод сопряженных градиентов
  • Метод двусопряженного градиента
  • Исключение Гаусса-Джордана

  Рекомендации

  • (fr) Эта статья частично или полностью взята из статьи в Википедии на английском языке под названием « Алгоритм Ланцоша » ( см. список авторов ) .
  • (en) Джин Х. Голуб и Чарльз Ф. Ван Лоан(en) , Matrix Computations , JHU Press, 1996 г. , 694 с. ( ISBN978-0-8018-5414-9 , онлайн-презентация )

  Внешние ссылки

  • (ru) Эндрю Й. Нг , Алиса Х. Чжэн и Майкл И. Джордан , «Анализ связей, собственные векторы и стабильность» , в IJCAI-01 , 2001 г. ( читать онлайн ) , стр. 903-910 : сравнение методов ранжирования HITS и PageRank (алгоритм Google) и их сходимость и устойчивость собственных векторов перед лицом изменений в наборах ссылок, организованных поисковыми системами.
  • (ru) Брайан А. Ламаккиа и Эндрю М. Одлызко , «Решение больших разреженных линейных систем над конечными полями» , Proceeding CRYPTO ’90 , p. 109-133 , § 3: Ланцоша и методы сопряженных градиентов.

  1 Свойства и структура алгоритмов

  Алгоритм Ланцоша — итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный Корнелием Ланцошом, адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности n с ограниченным числом итераций, m, где m много меньше n. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman [1] показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций m. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности — Алгоритм Ланцоша с выборочной ортогонализацией. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

  1.2 Математическое описание алгоритма

  Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица A=A T и вектор начального приближения b. Мы пытаемся найти трехдиагональную симметричную матрицу T_k=Q_k T AQ_k, собственные значения которой приближают собственные значения матрицы A. Иными словами на k-м шаге из ортонормированных векторов Ланцоша строится матрица Q_k = [q₁,q₂. q_k,] и в качестве приближенных собственных значений матрицы A принимаются числа Ритца. Пусть T_k=V [math]\Lambda[/math] V T есть спектральное разложение матрицы T_k, столбцы матрицы Q_kV рассматриваются как приближения к соответсвующим собственным векторам матрицы A [2] : [math] \begin q_1 = & \frac,\;\beta_0=0,\;q_0=0\\ for \; & j = 1 \; to \; k \\ & z = Aq_j\\ & \alpha_j = q^T_j z\\ & z=z-\alpha_jq_j-\beta_q_\\ & z = z — \sum^_\left(z^Tq_i\right)q_i,\;z = z — \sum^_\left(z^Tq_i\right)q_i \\ & \beta_j = \|z\|_2\\ & if \; \beta_j =0\; break\\ & q_ = \frac\\ end \; & for \end [/math]

  Диагональные элементы обозначены как [math]\alpha_j=t_[/math] , а элементы побочной диагонали [math]\beta_j=t_=t_[/math] . После каждой итерации мы вычисляем [math]\alpha_j, \beta_j[/math] , из которых строится матрица [math]T_j = \begin \alpha_1 & \beta_2 & & & & 0 \\ \beta_2 & \alpha_2 & \beta_3 & & & \\ & \beta_3 & \alpha_3 & \ddots & & \\ & & \ddots & \ddots & \beta_ & \\ & & & \beta_ & \alpha_ & \beta_j \\ 0 & & & & \beta_j & \alpha_j \\ \end[/math]

  Полная переортогонализация соответствует повторному проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что [math]z[/math] будет ортогонален векторам q₁,q₂. q_j-1. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем повторные копии сошедшихся чисел Ритца. То есть для больших k матрица T_k может иметь не одно собственное значение, очень близкое к [math]\lambda_i(A)[/math] , но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи MRRR алгоритма они могут быть вычислены за O(k 2 ) флопов.

  1.3 Вычислительное ядро алгоритма

  За один шаг алгоритма можно выделить следующие вычислительные ядра:

  • умножение исходной матрицы на вектор [math]z=Aq_j[/math]
  • процесс ортогонализации Грама-Шмидта [math]z = z — \sum^_\left(z^Tq_i\right)q_i[/math]

  1.4 Макроструктура алгоритма

  Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

  1.5 Схема реализации последовательного алгоритма

  В разделе #Математическое описание алгоритма уже представлено описание алгоритма.

  Пусть задана некая симметрическая матрица [math]A[/math]

  Шаг 1: Вычислить [math]q_1, \beta_0, q_0[/math]

  Шаг 2:В цикле последовательно вычислять матрицу [math]T_j,\alpha_j, \beta_j[/math]

  • Найти собственные значения матрицы [math]T_j[/math]
  • Если [math]\beta_j=0[/math] Выход

  

  Рисунок 1. Блок-схема алгоритма

  1.6 Последовательная сложность алгоритма

  Для вычисления k-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:

  • n+2k умножений вектора на вектор
  • 2k умножение вектора на константу
  • 2k сложения векторов
  • O(k 2 ) операций на разложение матрицы T_k

  Легко проверить, что m шагов этого алгоритма потребуют O(m 2 n) флопов и O(mn) слов памяти.

  1.7 Информационный граф

  

  Рисунок 2. Вычислительная схема j-й итерации. Синим цветом помечены промежуточные данные

  1.8 Ресурс параллелизма алгоритма

  Алгоритм Ланцоша с полной переортогонализацией является итерационным. На j-м шаге вычисляются числа [math]\alpha_j, \beta_j[/math] и ортонормированный вектор [math]q_[/math] , которые служат входными данными для j+1-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации j-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции. Умножение вещественной симметричной матрицы A размера [math]n \times n[/math] на вектор b длины n требует последовательного выполнения n ярусов умножений и сложений.

  1.9 Входные и выходные данные алгоритма

  Входные данные: матрица [math]A\in\mathbb^[/math] . Дополнительные ограничения:

  • [math]A[/math] – симметричная матрица, т. е. [math]a_= a_, i, j = 1, \ldots, n[/math] .

  Объём входных данных: [math]\frac[/math] .

  Выходные данные:

  • вектор приближенных собственных значений матрицы [math]\Lambda[/math] (элементы [math]\lambda_[/math] ).
  • приближения к собственным векторам [math]v_i,i=1,\ldots, m [/math] .

  Объём выходных данных: [math] m(n+1) [/math] .

  1.10 Свойства алгоритма

  • Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом k итераций оценивается в
  • Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в [math]\frac\approx 2\,k[/math] так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо [math]k\lt \lt n[/math]
  • Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
  • Часто уже при небольшом [math]( j \sim 2 \sqrt)[/math] числе шагов метода крайние пары Ритца [math](\lambda_i,v_i|[v_1,v_2. v_j]=Q_jV)[/math] оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
  • Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из l, l > 1, векторов. Аналогом векторов [math]z,q_i[/math] простого алгоритма в блочном являются прямоугольные матрицы размера n * l, аналогом скаляров [math]\beta_i, \alpha_i[/math] — квадратные матрицы порядка l, аналогом трехдиагональной матрицы T_i — ленточная матрица с полушириной ленты l. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до l. Поэтому блочный алгоритм с l начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность l. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и Алгоритм Ланцоша для точной арифметики может быть дополнен выборочной ортогонализацией.

  2 Программная реализация алгоритма

  2.1 Особенности реализации последовательного алгоритма

  2.2 Локальность данных и вычислений

  2.3 Возможные способы и особенности параллельной реализации алгоритма

  2.4 Масштабируемость алгоритма и его реализации

  Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки TRLan на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[1] Набор и границы значений изменяемых параметров реализации алгоритма:

  • размер матрицы от 5000 до 25000 с шагом 500.
  • число процессоров 1, 2, 4, [8:144] с шагом 8.

  Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:

  Алгоритм Ланцоша для точной арифметики (без переортогонализации)

  Алгоритм Ланцоша был опубликован венгерским математиком и физиком Корнелием Ланцошем [1] в 1950 году [2] . Этот метод является частным случаем алгоритма Арнольда в случае, если исходная матрица [math]A[/math] — симметрична, и был представлен как итерационный метод вычисления собственных значений симметричной матрицы. Этот метод позволяет за [math]k[/math] итераций вычислять [math]k[/math] приближений собственных значений и собственных векторов исходной матрицы. Хотя алгоритм и был эффективным в вычислительном смысле, но он на некоторое время был предан забвению из-за численной неустойчивости. Только в 1970 Ojalvo и Newman модифицировали алгоритм для использования в арифметике с плавающей точкой [3] . Новый метод получил название алгоритма Ланцоша с полной переортогонализацией.

  В данной статье рассматривается исходная версия алгоритма [4] .

  1.2 Математическое описание алгоритма

  Математическое описание алгоритма Ланцоша для вычисления собственных значений симметричной матрицы А требует 1) математического описания метода Ланцоща для построения крыловского подпространства и 2) математического описания т.н. процедуры Рэлея-Ритца.

  1.2.1 Метод Ланцоша для построения крыловского подпространства

  Пусть дано уравнение [math]Ax = b[/math] , причём вектор [math]b[/math] известен и имеется способ вычисления [math]Ax[/math] .

  Вычислим [math]y_ = b[/math] , [math]y_ = Ab[/math] , [math]y_ = A^b[/math] . [math]y_ = A^<(n - 1)>b[/math] , где [math]n[/math] — порядок матрицы [math]A[/math] . Определим матрицу [math]K = [y_. y_][/math] .

  Тогда [math]AK = [Ay_. Ay_, Ay_] = [y_. y_,A^n][/math] . В предположении, что матрица [math]K[/math] — невырождена, можно вычислить вектор [math]c = -K^A^[/math] . Поскольку первые [math]n — 1[/math] столбцов в [math]AK[/math] совпадают с последними [math]n — 1[/math] столбцами в [math]K[/math] , то справедливо: [math]AK = K[e_,e_. e_,-c] = KC[/math] , то есть [math]K^AK = C[/math] , где [math]C[/math] — верхняя хессенбергова матрица.

  Заменим [math]K[/math] ортогональной матрицей [math]Q[/math] так, что при любом [math]k[/math] линейные оболочки первых [math]k[/math] столбцов в [math]K[/math] и [math]Q[/math] — одно и то же подпространство. Это подпространство называется крыловским и точно определяется как [math]\kappa_(A, b)[/math] — линейная оболочка векторов [math]b[/math] , [math]Ab[/math] . [math]A^<(k - 1)>b[/math] . Суть алгоритма Ланцоша заключается в вычислении стольких первых столбцов в [math]Q[/math] , сколько необходимо для получения требуемого приближения к решению [math]Ax = b [/math] ( [math]Ax = <\lambda>x[/math] ).

  Используем QR-разложение матрицы [math]K = QQ[/math]
  Тогда [math]K^AK = (R^Q^)A(QR) = C[/math] , откуда [math]Q^AQ = RCR^ = H[/math]

  Матрица [math]H[/math] является верхней хессенберговой в силу того, что верхней хессенберговой является матрица [math]C[/math] , а матрицы [math]R[/math] и [math]R^[/math] — верхние треугольные.

  Однако алгоритм Ланцоша подразумевает симметричность матрицы [math]A[/math] . Положим [math]Q = [Q_, Q_][/math] , где [math]Q_ = [q_. q_][/math] и [math]Q_ = [q_. q_][/math] . В итоге получается, что [math]H = Q^AQ = [Q_,Q_]^A[Q_,Q_][/math] Из симметричности [math]A[/math] , таким образом, следует симметричность и трехдиагональность матрицы [math]H = T[/math] . Теперь можно вывести две основные формулы, используемые в методе Ланцоша:

  • [math]AQ_ = \beta_q_ + \alpha_q_ + \beta_q_[/math] — получается приравниванием столбцов [math]j[/math] в обеих частях равенства [math]AQ = QT[/math] .
  • [math]q_^Aq_ = \alpha_[/math] — получается домножением обеих частей предыдущего соотношения на [math]q_^[/math] и учетом ортогональности [math]Q[/math] .

  1.2.2 Процедура Рэлея-Ритца

  Пусть [math]Q = [Q_, Q_][/math] — произвольная ортогональная матрица порядка [math]n[/math] , причём [math]Q_[/math] и [math]Q_[/math] имеют размеры [math]nk[/math] и [math]n(n — k)[/math] .

  Пусть [math]T = Q^AQ = [Q_,Q_]^A[Q_,Q_][/math] . Обратим внимание, что в случае [math]k = 1[/math] выражение [math]T_ = Q_^AQ_[/math] — отношение Рэлея [math]\rho(Q_,A)[/math] , то есть число [math](Q_^AQ_) / (Q_^Q_)[/math] . Определим теперь суть процедуры Рэлея-Ритца:

  • Процедура Рэлея-Ритца — интерпретация собственных значений матрицы [math]T = Q^AQ[/math] как приближений к собственным значениям матрицы [math]A[/math] . Эти приближения называются числами Ритца. Пусть [math]T_ = V<\lambda>V^[/math] — спектральное разложение матрицы [math]T_[/math] . Столбцы матрицы [math]Q_V[/math] рассматриваются как приближения к соответствующим собственным векторам и называются векторами Ритца.

  Существует несколько важных фактов о процедуре Рэлея-Ритца, характеризующих получаемые приближения собственных значений матрицы [math]A[/math] .

  • Минимум величины [math]||AQ_ — Q_R||_[/math] по всем симметричным [math]k^2[/math] матрицам [math]R[/math] достигается при [math]R = T_[/math] . При этом [math]||AQ_ — Q_R||_ = ||T_||_[/math] . Пусть [math]T_ = V<\lambda>V^[/math] — спектральное разложение матрицы [math]T_[/math] . Минимум величины [math]||AP_ — P_D||_[/math] , когда [math]P_[/math] пробегает множество [math]nk [/math] матрицы с ортонормированными столбцами, таких, что [math]span(P_) = span(Q_)[/math] , а [math]D[/math] пробегает множество диагональных [math]k^2[/math] матриц также равен [math]||T_||_[/math] и достигается для [math]P_ = Q_V[/math] и [math]D = \lambda[/math] .

  Пусть [math]T_[/math] , [math]T_[/math] , [math]Q_[/math] — те же самые матрицы. Пусть [math]T_ = V<\lambda>V^[/math] — всё то же спектральное разложение, причём [math]V = [v_. v_][/math] — ортогональная матрица, а [math]\lambda = diag(\theta_. \theta_)[/math] . Тогда:

  • Найдутся [math]k[/math] необязательно наибольших собственных значений [math]\alpha_. \alpha_[/math] матрицы [math]A[/math] , такие, что [math]|\theta_ — \alpha_| \lt = ||T_||_[/math] для [math]i = 1. k[/math] . Если [math]Q_[/math] вычислена методом Ланцоша,то правая часть этого выражения равняется [math]\beta_[/math] , то есть единственному элементу блока [math]T_[/math] , который может быть отличен от нуля. Указанный эдемент находится в правом верхнем углу блока.
  • [math]||A(Q_v_) — (Q_v_)\theta_||_ = ||T_v_||_[/math] . Таким образом, разность между числом Ритца [math]\theta_[/math] и некоторым собственным значением [math]\alpha[/math] матрицы [math]A[/math] не превышает величины [math]||T_v_||_[/math] , которая может быть много меньше, чем [math]||T_||_[/math] . Если матрица [math]Q_[/math] вычислена алгоритмом Ланцоша, то эта величина равна [math]\beta_|v_(k)|[/math] , где [math]v_(k)[/math] — k-ая компонента вектора [math]v_[/math] . Таким образом, можно дещево вычислить норму невязки — левой части выражения, не выполняя ни одного умножения вектора на матрицы [math]Q_[/math] или [math]A[/math] .
  • Не имея информации о спектре матрицы [math]T_[/math] , нельзя дать какую-либо содержательную оценку для погрешности в векторе Ритца [math]Q_v_[/math] . Если известно, что [math]\theta_[/math] отделено расстоянием, не меньшим [math]g[/math] , от прочих собственных значений матриц [math]T_[/math] , [math]T_[/math] и матрица [math]Q_[/math] вычислена алгоритмом Ланцоша, то угол [math]\theta[/math] между [math]Q_v_[/math] и точным собственным вектором [math]A[/math] можно оценить формулой [math]\frac\sin \lt = \frac[/math]

  В алгоритме Ланцоша из ортонормированных векторов строится матрица [math]Q_ = [q_. q_][/math] и в качестве приближенных собственных значений [math]A[/math] принимаются числа Ритца — собственные значения симметричной трёхдиагональной матрицы [math]T_ = Q_^AQ_[/math] .

  1.3 Вычислительное ядро алгоритма

  Вычислительное ядро алгоритма Ланцоща состоит в получении на каждой итерации очередного промежуточного вектора [math]z[/math] , получаемого путём умножения исходной матрицы [math]A[/math] на вектор Ланцоша [math]q_[/math] , полученный на предыдущей итерации.

  Также при значениях [math]k[/math] , сопоставимых с [math]n[/math] , процедура вычисления собственных значений и собственных векторов симметричной трёхдиагональной матрицы по некоторому выбранному алгоритму подразумевает значительный объём расчётов и может считаться вычислительным ядром. Однако на практике метод используется прежде всего при малых [math]k[/math] .

  1.4 Макроструктура алгоритма

  Метод Ланцоша соединяет метод Ланцоша для построения крыловского подпространства в процедурой Рэлея-Ритца, подразумевающей вычисление собственных значений симметричной трёхдиагональной матрицы. Первая часть алгоритма представляет собой строго последовательные итерации, на каждой из которых вначале строится очередной столбец матрицы [math]Q_[/math] из первых [math]j[/math] ортонормированных векторов Ланцоша и матрица [math]T_ = Q^T_AQ_[/math] порядка [math]j[/math] . Итерационный процесс завершается, когда [math]j = k[/math] . Вторая часть алгоритма представляет собой поиск собственных значений и собственных векторов симметричной трёхдиагональной матрицы [math]T_[/math] , что реализуется отдельным алгоритмом, на практике используется QR-алгоритм.

  Особенность методов крыловского подпространства состоит в предположении, что матрица [math]A[/math] доступна в виде «черного ящика», а именно подпрограммы, вычисляющей по заданному вектору [math]z[/math] произведение [math]y = Az[/math] (а также, возможно, произведение [math]y = A^z[/math] , если матрица [math]A[/math] несимметрична). Другими словами, прямой доступ к элементам матрицы и их изменение не используются, т.к.

  • Умножение матрицы на вектор — наиболее дешевая нетривиальная операция, которую можно проделать с разреженной матрицей: в случае разреженной [math]A[/math] , содержащей [math]m[/math] ненулевых элементов, для матрично-векторного умножения нужны [math]m[/math] скалярных умножений и не более [math]m[/math] сложений.
  • [math]A[/math] может быть не представлена в виде матрицы явно, а доступна именно через подпрограмму для вычисления произведений [math]Ax[/math] .

  Таким образом, как QR-алгоритм поиска собственных значений, собственных векторов матрицы [math]T_[/math] и оценки погрешности в них результирующей трехдиагональной матрицы, так и умножение матрицы на вектор внутри каждой итерации могут быть рассмотрены в качестве отдельных макроопераций. Проблемой является то, что именно эти манипуляции порождают основную вычислительную сложность алгоритма. Поскольку матрица, к которой применяется QR-алгоритм, имеет порядок [math]k[/math] и на практике бывает мала, сделаем допущение о плотности матрицы [math]A[/math] . Тогда умножение этой матрицы на вектор внутри каждой итерации нужно будет реализовывать в рамках самого алгоритма. Еще одним важным замечанием является момент вычисления собственных значений матрицы [math]T_[/math] . Её собственные значения, полученные на итерации [math]p[/math] , никак не используются на итерации [math]p + 1[/math] , поэтому целесообразно вынести эту процедуру за основной цикл. Однако на практике в силу малых размеров матрицы [math]T_[/math] вычисление соответствующих величин производится непосредственно в теле основного цикла, что позволяет сразу оценить достигнутую точность вычислений.

  1.5 Схема последовательной реализации алгоритма

  Алгоритм итерационный, на каждой итерации выполняется вычисление очередного столбца матрицы [math]Q[/math] , а также вычисление элементов очередной пары (диагонального и околодиагонального) элементов матрицы [math]T[/math] с использованием значений столбца матрицы [math]Q[/math] , полученного на предыдущей итерации и значения околодиагонального элемента матрицы [math]T[/math] с предыдущей итерации. Вектор [math]z[/math] носит вспомогательный характер. Коэффициенты [math]\alpha_[/math] соответствуют диагональным элементам матрицы [math]T_[/math] , коэффициенты [math]\beta_[/math] – наддиагональным и поддиагональным элементам той же матрицы. Далее приводится пример псевдокода. По окончании итераций вычислеяются собственные значения и собственные векторы матрицы [math]T_[/math] , что в псевдокоде не отражается (т.к. это отдельный алгоритм).

  Вычислить столбец q1 матрицы Qk q1 = b / ||b|| q0 - нулевой вектор beta_0 = 0 foreach ( j = 1 .. k ) < z = A * qj alpha_j = qj * z z = z – alpha_j * qj – beta_(j-1) * q(j – 1) beta_j = ||z|| if ( beta_j == 0 ) < break // Все ненулевые собственные значения найдены. Выход может быть досрочным. >else < q(j + 1) = z/beta_j >>

  Вычислить собственные значения и собственные векторы симметричной трехдиагональной матрицы Tj наиболее подходящим методом. Полученный вектор Lj есть искомые собственные значения.

  1.6 Последовательная сложность алгоритма

  Все дальнейшие выкладки верны для наиболее быстрого последовательного варианта выполнения указываемых операций.

  • Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] требует [math]n^2[/math] умножений и сложений.
  • Перемножение векторов длины [math]n[/math] требует по [math]n[/math] умножений и сложений.
  • Поэлементное сложение векторов длины [math]n[/math] требует [math]n[/math] сложений.
  • Умножение вектора длины [math]n[/math] на число требует [math]n[/math] умножений.
  • Нахождение квадратичной нормы вектора длины [math]n[/math] требует по [math]n[/math] умножений и сложений, а также одну операцию извлечения квадратного корня.
  • Вычисление собственных значений матрицы порядка [math]k[/math] QR-алгоритмом требует [math]O(k^2)[/math] операций, вычисление также и собственных векторов требует примерно [math]6k^3[/math] , то есть [math]O(k^3)[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]O(k^)[/math] операций.

  Таким образом, для выполнения одной итерации метода Ланцоша требуется 1 операция вычисления квадратного корня, [math]n[/math] умножений, а также по [math]n ^ 2 + n + 2n + n[/math] сложений и умножений, а также суммарно [math]O(k ^ 3)[/math] операций, требуемых для поиска собственных значений и собственных векторов матрицы [math]T_[/math] . Таким образом последовательная сложность алгоритма Ланцоша составляет [math]O(n ^ 2) + O(k ^ 3)[/math] . Это, очевидно, меньше чем [math]O(n^3)[/math] операций, требуемых в алгоритмах вычисления всех собственных значений произвольных симметричных матриц, в чём и заключается эффективность применения метода Ланцоша.

  1.7 Информационный граф

  Информационный граф алгоритма, определяемый в [5] требует следующего замечания: подразумевается, что, например, умножение квадратной матрицы на вектор распараллеливается на [math]n[/math] потоков, где [math]n[/math] — порядок матрицы, а не на [math]n\log[/math] потоков, чего можно было бы достичь, используя суммирование сдваиванием. Это допущение аналогично допущению из статьи Умножение матрицы на вектор и Это же допущение будет использовано и при оценке ресурса параллелизма алгоритма. Приводится граф без отображения входных и выходных данных.

  

  Рис. 1: Граф алгоритма для матрицы А порядка 5 и k = 3. Op — одна итерация алгоритма, Eigen — процедура вычисления собственных значений

  

  Рис. 2: Внутренний граф итерации алгоритма для матрицы A порядка 5. Vector sum — нахождение суммы элементов вектора-поэлементного произведения векторов z и qj. Op — вычитания из z векторов qj и q(j — 1), домноженных на коэффициенты. Norm — вычисление нормы вектора z.

  1.8 Ресурс параллелизма алгоритма

  Каждая итерация алгоритма Ланцоша выполняется строго последовательно. Однако внутри итерации возможно распараллеливание алгоритма.

  • Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] требует последовательного выполнения [math]n[/math] ярусов умножений и сложений.
  • Все остальные операции в рамках одной итерации, ведущие к вычислению матрицы [math]T_[/math] выполняются строго последовательно. При этом вычисление значений векторов может быть выполнено за 1 ярус умножений / сложений, а вычисление чисел — за [math]\log[/math] таких операций.
  • Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма и рассматривается в соответствующей статье.

  Таким образом, при классификации по ширине ЯПФ алгоритм обладает линейной сложностью. В случае классификации по высоте ЯПФ алгоритм также имеет линейную сложность.

  1.9 Входные и выходные данные алгоритма

  Входные данные: квадратная матрица [math]A[/math] (элементы [math]a_[/math] ), вектор [math]b[/math] , число [math]k[/math]
  Объём входных данных: [math]n^2 + n + 1[/math]
  Выходные данные: вектор [math]L[/math] , матрица [math]E[/math] (элементы [math]e_[/math] )
  Объём выходных данных: [math]k + kn[/math]

  1.10 Свойства алгоритма

  В случае неограниченности ресурсов соотношение последовательной и параллельной сложности алгоритма Ланцоща по сумме всех итераций представляет собой [math]k ^ 3 + kn ^ 2[/math] к [math]X + kn[/math] , где [math]X[/math] — параллельная сложность алгоритма, используемого для вычисления собственных значений. В случае использования метода вычисления собственных значений с линейной сложностью при классисифкации по высоте ЯПФ это соотношение будет равно [math](k ^ 2 + n ^ 2) / (k + n)[/math] .

  Алгоритм Ланцоша не является полностью детерминированным в том смысле, что возможно выполнение числа итераций алгоритма, меньшего, чем заданное (из-за того, что вычислены все ненулевые собственные значения матрицы [math]A[/math] ).

  Важное свойство метода Ланцоша состоит в том, что первыми в матрице [math]T_[/math] появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы [math]A[/math] , находящихся на краях её спектра.

  2 Программная реализация алгоритма

  2.1 Особенности реализации последователього алгоритма

  2.2 Возможные способы и особенности параллельной реализации алгоритма

  Подробный анализ алгоритмов нахождения собственных значений, включая алгоритм Ланцоша, был проведен в [6] . Существует несколько как последовательных, так и параллельных реализаций алгоритма Ланцоша, включенных в программные библиотеки для поиска собственных значений матриц. Свойства этих реализаций можно увидеть в таблицах ниже. Первая таблица — относительно недавние реализации, вторая — «музейные экспонаты». Следует уделить особенное внимание столбцам «тип метода» и «параллелизация». В первом из них значение только N соответствует описанному варианту без реортогонализации, F – полной реортогонализации, P – частичной реортогонализации, S – выборочной реортогонализации. Во втором значние «none» соответствует отсутствию параллельной реализации, M – параллельной реализации посредством MPI, O – параллельной реализации посредством OpenMP. Их приведение целесообразно, так как на практике алгоритм Ланцоша без переортогонализации неустойчив. Также указываются библиотеки, в которых реализованы более глубокие модификации метода Ланцоша, с указанием изменений в графе «тип метода».

  Таблица 1 — «современные» реализации.

  Название библиотеки	Язык программирования	Дата появления, версия библиотки	Тип метода	Параллелизация
  BLKLAN	C/Matlab	2003	P	none
  BLZPACK	F77	2000, 04/00	P + S	M
  IETL	C++	2006, 2.2	N	none
  SLEPc	C/F77	2009, 3.0.0	All	M
  TRLAN	F90	2006	Dynamic thick-restart	M
  PROPACK	F77/Matlab	2005, 2.1 / 1.1	P, finds SVD	O
  IRBLEIGS	Matlab	2000, 1.0	Indefinitie symmetric	none

  Таблица 2 — «исторические» реализации.

  Название библиотеки	Язык программирования	Дата появления, версия библиотки	Тип метода	Параллелизация
  ARPACK	F77	1995, 2.1	Arnoldi iterations, impliicit restart	M
  ARPACK++	C++	1998, 1.1	Arnoldi iterations, impliicit restart	none
  LANCZOS	F77	1992	N	none
  LANZ	F77	1991, 1.0	P	none
  LASO	F77	1983, 2	S	none
  NAPACK	F77	1987	N	none
  QMRPACK	F77	1996	Designed for nonsymmetrical matrices (lookahead)	none
  SVDPACK	C/F77	1992	P, finds SVD	none
  Underwood	F77	1975	F, block version	none

  В настоящее время алгоритм Ланцоша поиска собственных значений квадратной симметричной матрицы включён в несколько библиотек и реализован на различных языках, среди которых такие наиболее распространенные как C, C++, FORTRAN77/90, MathLab. Ссылки на наиболее проверенные и часто используемые реализации алгоритма:

  • Так, например, в языке MathLab во встроенном пакете ARPACK найти собственные значения методом Ланцоша можно при помощи вызова функции eigs().
  • На языке С существует библиотека PRIMME, название которой расшифровывается как PReconditioned Iterative MultiMethod Eigensolver (итерационные методы поиска собственных значений с предусловием).
  • Библиотека NAG Numerical Library также содержит обширный набор функций, позволяющих проводить математический анализ, среди которых есть и реализация алгоритма Ланцоша. Библиотека доступна в трех пакетах: NAG C Library, NAG Fortran Library и NAG Library for .NET, совместна со многими языками и с основными операционными системами (Windows, Linux и OS X, а также Solaris, AIX и HP-UX).

  Заслуживают внимания также реализации в проектах IETL [7] , ARPACK [8] и SLEPc [9] .

  2.3 Результаты прогонов

  2.4 Выводы для классов архитектур

  3 Литература

  1. ↑https://ru.wikipedia.org/wiki/Ланцош,_Корнелий
  2. ↑ Lanczos, C. «An iteration method for the solution of the eigenvalue problem of linear differential and integral operators», J. Res. Nat’l Bur. Std. 45, 255-282 (1950).
  3. ↑ Ojalvo, I.U. and Newman, M., «Vibration modes of large structures by an automatic matrix-reduction method», AIAA J., 8 (7), 1234–1239 (1970).
  4. ↑ Деммель Дж. Вычислительная линейная алгебра. Теория и приложения — М.: Мир, 2001. (см.)
  5. ↑ Параллельные вычисления (Воеводин В.В., Воеводин Вл.В.) — Спб, изд-во «БХВ-Петербург», 2002 (см.)
  6. ↑ V. Hernandez, J. E. Roman, A. Tomas, V. Vidal. A Survey of Software for Sparse Eigenvalue Problems (см.)
  7. ↑http://www.comp-phys.org/software/ietl/lanczos.html
  8. ↑http://www.caam.rice.edu/software/ARPACK/
  9. ↑ Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.
  • Уровень алгоритма
  • Статьи в работе

  Методы собственных чисел в RF-STABILITY

  Das Zusatzmodul RF-STABIL ermittelt die Verzweigungslastfaktoren, Knicklängen und Eigenformen von RFEM-Modellen. Die Stabilitätsuntersuchungen können dabei nach verschiedenen Eigenwertmethoden erfolgen, die je nach System und Rechnerkonfiguration ihre Vorteile haben.

  Метод Ланцоша

  Собственные числа определяются напрямую. С помощью этого алгоритма обычно можно достичь быстрой сходимости. Данный метод подходит для стандартных моделей и поэтому установлен по умолчанию.

  Корни нормативного многочлена

  Данный метод также основан на методе прямого расчета. Для более крупных конструктивных систем этот метод может быть быстрее, чем метод Ланцоша. Основное преимущество — точность вычисления высших собственных чисел.

  Метод итерации субпространства

  Все собственные значения определяются за один шаг. Спектр матрицы жесткости оказывает сильное влияние на длительность расчета. Поскольку матрица жесткости хранится в оперативной памяти, данный метод не подходит для сложных систем. Кроме того, могут отображаться отрицательные коэффициенты критической нагрузки.

  Метод итерации ICG

  Метод неполного сопряженного градиента требует небольшого объема оперативной памяти. Поскольку собственные значения определяются одно за другим, для расчета малых и средних конструктивных систем требуется больше времени по сравнению с прямым методом. Однако спектр не влияет на длительность расчета. Метод ICG подходит для анализа очень больших систем с небольшим количеством собственных значений. Данный метод не дает никаких отрицательных коэффициентов критической нагрузки.

  Похожие публикации:

  1. К чему приводит токовая перегрузка проводки
  2. Курсы 3д моделирования в цивил геология москва
  3. Магистральная и радиальная схема электроснабжения
  4. Создание лестницы в ревит