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

Для алгоритма поиска нулевого пространства см. блок-алгоритм Ланцоша .

Информацию о методе интерполяции см. в разделе Передискретизация Ланцоша .

Информацию об аппроксимации гамма-функции см. в разделе «Аппроксимация Ланцоша» .

Алгоритм Ланцоша — это итерационный метод, разработанный Корнелиусом Ланцошем , который представляет собой адаптацию степенных методов для поиска $m$ «наиболее полезные» (стремящиеся к крайне высокому/самому низкому) собственным значениям и собственным векторам $n\times n$ Эрмитова матрица , где $m$ часто, но не обязательно, намного меньше, чем $n$. Несмотря на свою вычислительную эффективность в принципе, метод в его первоначальной формулировке оказался бесполезным из-за его численной нестабильности .

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

В своей оригинальной работе эти авторы также предложили, как выбрать начальный вектор (т.е. использовать генератор случайных чисел для выбора каждого элемента начального вектора) и предложили эмпирически определяемый метод определения $m$, уменьшенное количество векторов (т.е. оно должно быть выбрано примерно в 1,5 раза больше желаемого количества точных собственных значений). Вскоре после этого за их работой последовала Пейдж, которая также провела анализ ошибок. В 1988 году Охалво представил более подробную историю этого алгоритма и эффективный тест на ошибку собственных значений.

Алгоритм

Введите матрицу эрмитову $A$ размера $n\times n$и, возможно, несколько итераций $m$ (по умолчанию пусть $m=n$).

• Строго говоря, алгоритму не нужен доступ к явной матрице, а только функция $v\mapsto Av$ который вычисляет произведение матрицы на произвольный вектор. Эта функция вызывается не более чем $m$ раз.

Выведите ​ $n\times m$ матрица $V$ с ортонормированными столбцами и трехдиагональной вещественной симметричной матрицей $T=V^{*}AV$ размера $m\times m$. Если $m=n$, затем $V$ является унитарным , и $A=VTV^{*}$.

Предупреждение. Итерация Ланцоша склонна к числовой нестабильности. При выполнении неточной арифметики следует принять дополнительные меры (как описано в последующих разделах) для обеспечения достоверности результатов.

1. Позволять $v_{1}\in \mathbb {C} ^{n}$ — произвольный вектор с евклидовой нормой $1$.
2. Сокращенный начальный шаг итерации:
   1. Позволять $w_{1}'=Av_{1}$.
   2. Позволять $\alpha _{1}=w_{1}'^{*}v_{1}$.
   3. Позволять $w_{1}=w_{1}'-\alpha _{1}v_{1}$.
3. Для $j=2,\dots ,m$ делать:
   1. Позволять $\beta _{j}=\|w_{j-1}\|$ (также евклидова норма ).
   2. Если $\beta _{j}\neq 0$, тогда пусть $v_{j}=w_{j-1}/\beta _{j}$,

      еще выбери как $v_{j}$ произвольный вектор с евклидовой нормой $1$ что ортогонально всем $v_{1},\dots ,v_{j-1}$.

   3. Позволять $w_{j}'=Av_{j}$.
   4. Позволять $\alpha _{j}=w_{j}'^{*}v_{j}$.
   5. Позволять $w_{j}=w_{j}'-\alpha _{j}v_{j}-\beta _{j}v_{j-1}$.
4. Позволять $V$ быть матрицей со столбцами $v_{1},\dots ,v_{m}$. Позволять $T={\begin{pmatrix}\alpha _{1}&\beta _{2}&&&&0\\\beta _{2}&\alpha _{2}&\beta _{3}&&&\\&\beta _{3}&\alpha _{3}&\ddots &&\\&&\ddots &\ddots &\beta _{m-1}&\\&&&\beta _{m-1}&\alpha _{m-1}&\beta _{m}\\0&&&&\beta _{m}&\alpha _{m}\\\end{pmatrix}}$.

Примечание $Av_{j}=w_{j}'=\beta _{j+1}v_{j+1}+\alpha _{j}v_{j}+\beta _{j}v_{j-1}$ для $2<j<m$.

В принципе существует четыре способа написания итерационной процедуры. Пейдж и другие работы показывают, что указанный выше порядок операций является наиболее численно стабильным. На практике исходный вектор $v_{1}$ можно рассматривать как еще один аргумент процедуры, при этом $\beta _{j}=0$ и индикаторы числовой неточности, включенные в качестве дополнительных условий завершения цикла.

Не считая умножения матрицы на вектор, каждая итерация $O(n)$ арифметические операции. Умножение матрицы на вектор можно выполнить в $O(dn)$ арифметические операции, где $d$ — среднее количество ненулевых элементов в строке. Таким образом, общая сложность $O(dmn)$, или $O(dn^{2})$ если $m=n$; алгоритм Ланцоша может быть очень быстрым для разреженных матриц. Схемы улучшения численной стабильности обычно оцениваются по этой высокой производительности.

Векторы $v_{j}$ называются векторами Ланцоша . Вектор $w_{j}'$ не используется после $w_{j}$ вычисляется, и вектор $w_{j}$ не используется после $v_{j+1}$ вычисляется. Следовательно, можно использовать одно и то же хранилище для всех трех. Аналогично, если только трехдиагональная матрица $T$ ищется, то необработанная итерация не требуется $v_{j-1}$ после того, как вычислил $w_{j}$, хотя некоторым схемам улучшения численной стабильности это понадобится позже. Иногда последующие векторы Ланцоша пересчитываются из $v_{1}$ когда это необходимо.

Применение к собственной проблеме

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

Если $\lambda$ является собственным значением $T$, и $x$ его собственный вектор ( $Tx=\lambda x$), затем $y=Vx$ является соответствующим собственным вектором $A$ с тем же собственным значением:

$${\begin{aligned}Ay&=AVx\\&=VTV^{*}Vx\\&=VTIx\\&=VTx\\&=V(\lambda x)\\&=\lambda Vx\\&=\lambda y.\end{aligned}}$$

Таким образом, алгоритм Ланцоша преобразует задачу собственного разложения для $A$ в задачу собственного разложения для $T$.

1. Для трехдиагональных матриц существует ряд специализированных алгоритмов, часто с большей вычислительной сложностью, чем алгоритмы общего назначения. Например, если $T$ это $m\times m$ трехдиагональная симметричная матрица тогда:
   • Непрерывная рекурсия позволяет вычислить характеристический полином в $O(m^{2})$ операций и оценивая их в определенный момент $O(m)$ операции.
   • Алгоритм собственных значений «разделяй и властвуй» можно использовать для вычисления всего собственного разложения $T$ в $O(m^{2})$ операции.
   • Метод быстрых мультиполей может вычислить все собственные значения всего за $O(m\log m)$ операции.
2. некоторые общие алгоритмы собственного разложения, в частности алгоритм QR Известно, что , сходятся быстрее для трехдиагональных матриц, чем для общих матриц. Асимптотическая сложность трехдиагонального QR равна $O(m^{2})$ так же, как и для алгоритма «разделяй и властвуй» (хотя постоянный коэффициент может быть и другим); поскольку собственные векторы вместе имеют $m^{2}$ элементов, это асимптотически оптимально.
3. Даже алгоритмы, на скорость сходимости которых не влияют унитарные преобразования, такие как степенной метод и обратная итерация , могут получить небольшой выигрыш в производительности от применения к трехдиагональной матрице. $T$ вместо исходной матрицы $A$. С $T$ очень разреженный, все ненулевые элементы находятся в очень предсказуемых позициях, он обеспечивает компактное хранение с превосходной производительностью по сравнению с кэшированием . Так же, $T$ является вещественной матрицей со всеми собственными векторами и собственными значениями, тогда как $A$ вообще может иметь комплексные элементы и собственные векторы, поэтому реальной арифметики достаточно для нахождения собственных векторов и собственных значений $T$.
4. Если $n$ очень велика, то уменьшается $m$ так что $T$ имеет управляемый размер, все равно позволит найти более крайние собственные значения и собственные векторы $A$; в $m\ll n$ Алгоритм Ланцоша можно рассматривать как схему сжатия с потерями для эрмитовых матриц, в которой особое внимание уделяется сохранению крайних собственных значений.

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

Приложение к трехдиагонализации

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

Аспекты, в которых различаются два алгоритма, включают:

• Ланцош воспользовался $A$ является разреженной матрицей, тогда как Householder этого не делает и будет генерировать fill-in .
• Ланчош полностью работает с оригинальной матрицей. $A$ (и у него нет проблем с тем, что она известна только неявно), тогда как необработанный Хаусхолдер хочет изменить матрицу во время вычисления (хотя этого можно избежать).
• Каждая итерация алгоритма Ланцоша создает еще один столбец окончательной матрицы преобразования. $V$, тогда как итерация Householder создает еще один фактор в унитарной факторизации $Q_{1}Q_{2}\dots Q_{n}$ из $V$. Однако каждый фактор определяется одним вектором, поэтому требования к хранению одинаковы для обоих алгоритмов. $V=Q_{1}Q_{2}\dots Q_{n}$ можно вычислить в $O(n^{3})$ время.
• Домохозяин численно стабилен, а необработанный Ланцош — нет.
• Ланцош во многом аналогичен, только $O(n)$ точки синхронизации (вычисления $\alpha _{j}$ и $\beta _{j}$). Домохозяин менее параллелен, имеет последовательность $O(n^{2})$ вычисляются скалярные величины, каждая из которых зависит от предыдущей величины в последовательности.

Вывод алгоритма

Существует несколько направлений рассуждений, которые приводят к алгоритму Ланцоша.

Более экономный метод питания

Основная статья: Итерация мощности

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

1. Выберите случайный вектор $u_{1}\neq 0$.
2. Для $j\geqslant 1$ (до направления $u_{j}$ сошлось) делаем:
   1. Позволять $u_{j+1}'=Au_{j}.$
   2. Позволять $u_{j+1}=u_{j+1}'/\|u_{j+1}'\|.$

• В большом $j$ предел, $u_{j}$ приближается к нормированному собственному вектору, соответствующему собственному значению наибольшей величины.

Критика, которую можно выдвинуть против этого метода, заключается в том, что он расточителен: он тратит много работы (матрично-векторные продукты на этапе 2.1) на извлечение информации из матрицы. $A$, но обращает внимание только на самый последний результат; реализации обычно используют одну и ту же переменную для всех векторов $u_{j}$, при этом каждая новая итерация перезаписывает результаты предыдущей. Вместо этого может быть желательно сохранить все промежуточные результаты и систематизировать данные.

Одна часть информации, которую тривиально можно получить из векторов $u_{j}$ есть цепочка подпространств Крылова . Один из способов заявить, что без введения множеств в алгоритм — это заявить, что он вычисляет

подмножество $\{v_{j}\}_{j=1}^{m}$ основе $\mathbb {C} ^{n}$ такой, что $Ax\in \operatorname {span} (v_{1},\dotsc ,v_{j+1})$ для каждого $x\in \operatorname {span} (v_{1},\dotsc ,v_{j})$ и все $1\leqslant j<m;$

это тривиально удовлетворяется $v_{j}=u_{j}$ пока $u_{j}$ линейно не зависит от $u_{1},\dotsc ,u_{j-1}$ (и в случае наличия такой зависимости можно продолжить последовательность, выбрав в качестве $v_{j}$ произвольный вектор, линейно независимый от $u_{1},\dotsc ,u_{j-1}$). Основа, содержащая $u_{j}$ векторов, однако, вероятно, будет численно плохо обусловленным , поскольку эта последовательность векторов по замыслу предназначена для сходимости к собственному вектору $A$. Чтобы избежать этого, можно объединить степенную итерацию с процессом Грама – Шмидта , чтобы вместо этого создать ортонормированный базис этих подпространств Крылова.

1. Выберите случайный вектор $u_{1}$ евклидовой нормы $1$. Позволять $v_{1}=u_{1}$.
2. Для $j=1,\dotsc ,m-1$ делать:
   1. Позволять $u_{j+1}'=Au_{j}$.
   2. Для всех $k=1,\dotsc ,j$ позволять $g_{k,j}=v_{k}^{*}u_{j+1}'$. (это координаты $Au_{j}=u_{j+1}'$ относительно базисных векторов $v_{1},\dotsc ,v_{j}$.)
   3. Позволять $w_{j+1}=u_{j+1}'-\sum _{k=1}^{j}g_{k,j}v_{k}$. (Отменить компонент $u_{j+1}'$ это в $\operatorname {span} (v_{1},\dotsc ,v_{j})$.)
   4. Если $w_{j+1}\neq 0$ тогда пусть $u_{j+1}=u_{j+1}'/\|u_{j+1}'\|$ и $v_{j+1}=w_{j+1}/\|w_{j+1}\|$,

      в противном случае выберите как $u_{j+1}=v_{j+1}$ произвольный вектор евклидовой нормы $1$ что ортогонально всем $v_{1},\dotsc ,v_{j}$.

Связь между векторами итераций мощности $u_{j}$ и ортогональные векторы $v_{j}$ это что

$Au_{j}=\|u_{j+1}'\|u_{j+1}=u_{j+1}'=w_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}=\|w_{j+1}\|v_{j+1}+\sum _{k=1}^{j}g_{k,j}v_{k}$.

Здесь можно заметить, что на самом деле нам не нужна $u_{j}$ векторы для вычисления этих $v_{j}$, потому что $u_{j}-v_{j}\in \operatorname {span} (v_{1},\dotsc ,v_{j-1})$ и, следовательно, разница между $u_{j+1}'=Au_{j}$ и $w_{j+1}'=Av_{j}$ находится в $\operatorname {span} (v_{1},\dotsc ,v_{j})$, что компенсируется процессом ортогонализации. Таким образом, тот же базис для цепочки подпространств Крылова вычисляется формулой

1. Выберите случайный вектор $v_{1}$ евклидовой нормы $1$.
2. Для $j=1,\dotsc ,m-1$ делать:
   1. Позволять $w_{j+1}'=Av_{j}$.
   2. Для всех $k=1,\dotsc ,j$ позволять $h_{k,j}=v_{k}^{*}w_{j+1}'$.
   3. Позволять $w_{j+1}=w_{j+1}'-\sum _{k=1}^{j}h_{k,j}v_{k}$.
   4. Позволять $h_{j+1,j}=\|w_{j+1}\|$.
   5. Если $h_{j+1,j}\neq 0$ тогда пусть $v_{j+1}=w_{j+1}/h_{j+1,j}$,

      в противном случае выберите как $v_{j+1}$ произвольный вектор евклидовой нормы $1$ что ортогонально всем $v_{1},\dotsc ,v_{j}$.

Априори коэффициенты $h_{k,j}$ удовлетворить

$Av_{j}=\sum _{k=1}^{j+1}h_{k,j}v_{k}$ для всех $j<m$;

определение $h_{j+1,j}=\|w_{j+1}\|$ может показаться немного странным, но соответствует общей схеме $h_{k,j}=v_{k}^{*}w_{j+1}'$ с

$v_{j+1}^{*}w_{j+1}'=v_{j+1}^{*}w_{j+1}=\|w_{j+1}\|v_{j+1}^{*}v_{j+1}=\|w_{j+1}\|.$

Поскольку векторы итераций мощности $u_{j}$ которые были исключены из этой рекурсии, удовлетворяют $u_{j}\in \operatorname {span} (v_{1},\ldots ,v_{j}),$ векторы $\{v_{j}\}_{j=1}^{m}$ и коэффициенты $h_{k,j}$ содержать достаточно информации из $A$ что все $u_{1},\ldots ,u_{m}$ можно вычислить, поэтому при переключении векторов ничего не потеряется. (Действительно, оказывается, что собранные здесь данные дают значительно лучшее приближение наибольшего собственного значения, чем те, которые можно получить при равном количестве итераций степенного метода, хотя на данный момент это не обязательно очевидно.)

Эта последняя процедура является итерацией Арнольди . Алгоритм Ланцоша возникает как упрощение, получаемое за счет исключения шагов вычислений, которые оказываются тривиальными, когда $A$ является эрмитовым, в частности, большая часть $h_{k,j}$ коэффициенты оказываются равными нулю.

Элементарно, если $A$ тогда является эрмитовым

$h_{k,j}=v_{k}^{*}w_{j+1}'=v_{k}^{*}Av_{j}=v_{k}^{*}A^{*}v_{j}=(Av_{k})^{*}v_{j}.$

Для $k<j-1$ мы знаем это $Av_{k}\in \operatorname {span} (v_{1},\ldots ,v_{j-1})$, и поскольку $v_{j}$ по построению ортогонально этому подпространству, это скалярное произведение должно быть равно нулю. (По сути, это также причина, по которой последовательностям ортогональных полиномов всегда можно задать трехчленное рекуррентное соотношение .) Для $k=j-1$ каждый получает

$h_{j-1,j}=(Av_{j-1})^{*}v_{j}={\overline {v_{j}^{*}Av_{j-1}}}={\overline {h_{j,j-1}}}=h_{j,j-1}$

поскольку последнее действительно, поскольку является нормой вектора. Для $k=j$ каждый получает

$h_{j,j}=(Av_{j})^{*}v_{j}={\overline {v_{j}^{*}Av_{j}}}={\overline {h_{j,j}}},$

значит, это тоже реально.

Более абстрактно, если $V$ это матрица со столбцами $v_{1},\ldots ,v_{m}$ тогда цифры $h_{k,j}$ могут быть идентифицированы как элементы матрицы $H=V^{*}AV$, и $h_{k,j}=0$ для $k>j+1;$ матрица $H$ находится верхний Хессенберг . С

$H^{*}=\left(V^{*}AV\right)^{*}=V^{*}A^{*}V=V^{*}AV=H$

матрица $H$ является эрмитовым. Это подразумевает, что $H$ также является нижним Хессенбергом, поэтому на самом деле он должен быть трехдиагональным. Будучи эрмитовой, ее главная диагональ действительна, а поскольку ее первая поддиагональ действительна по построению, то же самое верно и для ее первой наддиагонали. Поэтому, $H$ является реальной симметричной матрицей — матрицей $T$ спецификации алгоритма Ланцоша.

Одновременная аппроксимация крайних собственных значений

Один из способов характеристики собственных векторов эрмитовой матрицы. $A$ является стационарными точками фактора Рэлея

$r(x)={\frac {x^{*}Ax}{x^{*}x}},\qquad x\in \mathbb {C} ^{n}.$

В частности, наибольшее собственное значение $\lambda _{\max }$ это глобальный максимум $r$ и наименьшее собственное значение $\lambda _{\min }$ это глобальный минимум $r$.

В низкомерном подпространстве ${\mathcal {L}}$ из $\mathbb {C} ^{n}$ возможно найти максимум $x$ и минимум $y$ из $r$. Повторяя это для возрастающей цепочки ${\mathcal {L}}_{1}\subset {\mathcal {L}}_{2}\subset \cdots$ создает две последовательности векторов: $x_{1},x_{2},\ldots$ и $y_{1},y_{2},\dotsc$ такой, что $x_{j},y_{j}\in {\mathcal {L}}_{j}$ и

${\begin{aligned}r(x_{1})&\leqslant r(x_{2})\leqslant \cdots \leqslant \lambda _{\max }\\r(y_{1})&\geqslant r(y_{2})\geqslant \cdots \geqslant \lambda _{\min }\end{aligned}}$

Тогда возникает вопрос, как выбрать подпространства так, чтобы эти последовательности сходились с оптимальной скоростью.

От $x_{j}$, оптимальное направление поиска больших значений $r$ это градиент $\nabla r(x_{j})$, и аналогично из $y_{j}$ оптимальное направление, в котором следует искать меньшие значения $r$ это отрицательный градиент $-\nabla r(y_{j})$. В общем

$\nabla r(x)={\frac {2}{x^{*}x}}(Ax-r(x)x),$

поэтому интересующие направления достаточно легко вычислить с помощью матричной арифметики, но если кто-то хочет улучшить оба $x_{j}$ и $y_{j}$ то есть два новых направления, которые следует принять во внимание: $Ax_{j}$ и $Ay_{j};$ с $x_{j}$ и $y_{j}$ могут быть линейно независимыми векторами (действительно, близкими к ортогональным), вообще нельзя ожидать $Ax_{j}$ и $Ay_{j}$ быть параллельным. Нет необходимости увеличивать размерность ${\mathcal {L}}_{j}$ к $2$ на каждом шагу, если $\{{\mathcal {L}}_{j}\}_{j=1}^{m}$ считаются подпространствами Крылова, поскольку тогда $Az\in {\mathcal {L}}_{j+1}$ для всех $z\in {\mathcal {L}}_{j},$ таким образом, в частности для обоих $z=x_{j}$ и $z=y_{j}$.

Другими словами, мы можем начать с произвольного начального вектора $x_{1}=y_{1},$ построить векторные пространства

${\mathcal {L}}_{j}=\operatorname {span} (x_{1},Ax_{1},\ldots ,A^{j-1}x_{1})$

а затем искать $x_{j},y_{j}\in {\mathcal {L}}_{j}$ такой, что

$r(x_{j})=\max _{z\in {\mathcal {L}}_{j}}r(z)\qquad {\text{and}}\qquad r(y_{j})=\min _{z\in {\mathcal {L}}_{j}}r(z).$

Поскольку $j$итерация метода мощности $u_{j}$ принадлежит ${\mathcal {L}}_{j},$ из этого следует, что итерация для создания $x_{j}$ и $y_{j}$ не может сходиться медленнее, чем у степенного метода, и достигнет большего за счет аппроксимации обоих крайних значений собственных значений. Для подзадачи оптимизации $r$ на некоторых ${\mathcal {L}}_{j}$, удобно иметь ортонормированный базис $\{v_{1},\ldots ,v_{j}\}$ для этого векторного пространства. Таким образом, мы снова приходим к задаче итеративного вычисления такого базиса последовательности подпространств Крылова.

Конвергенция и другая динамика

При анализе динамики алгоритма удобно брать собственные значения и собственные векторы $A$ как задано, даже если они явно не известны пользователю. Чтобы исправить обозначения, пусть $\lambda _{1}\geqslant \lambda _{2}\geqslant \dotsb \geqslant \lambda _{n}$ — собственные значения (всем известно, что они вещественны и, следовательно, их можно упорядочить) и пусть $z_{1},\dotsc ,z_{n}$ — ортонормированный набор собственных векторов такой, что $Az_{k}=\lambda _{k}z_{k}$ для всех $k=1,\dotsc ,n$.

Также удобно зафиксировать обозначения для коэффициентов исходного вектора Ланцоша $v_{1}$ относительно этого собственного базиса; позволять $d_{k}=z_{k}^{*}v_{1}$ для всех $k=1,\dotsc ,n$, так что $\textstyle v_{1}=\sum _{k=1}^{n}d_{k}z_{k}$. Стартовый вектор $v_{1}$ обеднение некоторого собственного компонента приведет к задержке сходимости к соответствующему собственному значению, и хотя это просто проявляется как постоянный фактор в границах ошибки, истощение остается нежелательным. Одним из распространенных способов избежать постоянного поражения им является выбор $v_{1}$ сначала рисуя элементы случайным образом в соответствии с тем же нормальным распределением со средним значением $0$ а затем измените масштаб вектора до нормы $1$. До масштабирования это приводит к тому, что коэффициенты $d_{k}$ также быть независимыми нормально распределенными стохастическим переменными из одного и того же нормального распределения (поскольку замена координат унитарна), а после масштабирования вектора $(d_{1},\dotsc ,d_{n})$ будет иметь равномерное распределение на единичной сфере в $\mathbb {C} ^{n}$. Это позволяет ограничить вероятность того, что, например, $|d_{1}|<\varepsilon$.

Тот факт, что алгоритм Ланцоша не зависит от координат (операции рассматривают только внутренние произведения векторов, а не отдельные элементы векторов), позволяет легко создавать примеры с известной собственной структурой для запуска алгоритма: make $A$ диагональная матрица с искомыми собственными значениями на диагонали; пока стартовый вектор $v_{1}$ имеет достаточно ненулевых элементов, алгоритм выведет общую трехдиагональную симметричную матрицу как $T$.

Теория конвергенции Каниэля – Пейджа

После $m$ итерационные шаги алгоритма Ланцоша, $T$ это $m\times m$ реальная симметричная матрица, которая, как и выше, имеет $m$ собственные значения $\theta _{1}\geqslant \theta _{2}\geqslant \dots \geqslant \theta _{m}.$ Под конвергенцией понимают прежде всего конвергенцию $\theta _{1}$ к $\lambda _{1}$ (и симметричная сходимость $\theta _{m}$ к $\lambda _{n}$) как $m$ растет, и во вторую очередь сближение некоторого диапазона $\theta _{1},\ldots ,\theta _{k}$ собственных значений $T$ своим коллегам $\lambda _{1},\ldots ,\lambda _{k}$ из $A$. Сходимость алгоритма Ланцоша часто на несколько порядков быстрее, чем сходимость алгоритма степенной итерации.

Границы для $\theta _{1}$ исходят из приведенной выше интерпретации собственных значений как крайних значений коэффициента Рэлея. $r(x)$. С $\lambda _{1}$ априори является максимумом $r$ в целом $\mathbb {C} ^{n},$ тогда как $\theta _{1}$ это просто максимум на $m$-мерное подпространство Крылова, тривиально получаем $\lambda _{1}\geqslant \theta _{1}$. И наоборот, любая точка $x$ в этом подпространстве Крылова имеется нижняя оценка $r(x)$ для $\theta _{1}$, поэтому, если можно указать точку, для которой $\lambda _{1}-r(x)$ мал, то это обеспечивает жесткую границу $\theta _{1}$.

Размер $m$ Krylov subspace is

$\operatorname {span} \left\{v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}\right\},$

поэтому любой его элемент может быть выражен как $p(A)v_{1}$ для некоторого полинома $p$ степени максимум $m-1$; коэффициенты этого полинома - это просто коэффициенты линейной комбинации векторов $v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}$. У полинома, который мы хотим, будут действительные коэффициенты, но на данный момент нам следует учесть и комплексные коэффициенты, и мы напишем $p^{*}$ для многочлена, полученного комплексным сопряжением всех коэффициентов $p$. В этой параметризации подпространства Крылова имеем

$r(p(A)v_{1})={\frac {(p(A)v_{1})^{*}Ap(A)v_{1}}{(p(A)v_{1})^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p(A)^{*}Ap(A)v_{1}}{v_{1}^{*}p(A)^{*}p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A^{*})Ap(A)v_{1}}{v_{1}^{*}p^{*}(A^{*})p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A)Ap(A)v_{1}}{v_{1}^{*}p^{*}(A)p(A)v_{1}}}$

Используя теперь выражение для $v_{1}$ как линейную комбинацию собственных векторов, мы получаем

$Av_{1}=A\sum _{k=1}^{n}d_{k}z_{k}=\sum _{k=1}^{n}d_{k}\lambda _{k}z_{k}$

и вообще

$q(A)v_{1}=\sum _{k=1}^{n}d_{k}q(\lambda _{k})z_{k}$

для любого многочлена $q$.

Таким образом

$\lambda _{1}-r(p(A)v_{1})=\lambda _{1}-{\frac {v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})\lambda _{k}p(\lambda _{k})z_{k}}{v_{1}^{*}\sum _{k=1}^{n}d_{k}p^{*}(\lambda _{k})p(\lambda _{k})z_{k}}}=\lambda _{1}-{\frac {\sum _{k=1}^{n}|d_{k}|^{2}\lambda _{k}p(\lambda _{k})^{*}p(\lambda _{k})}{\sum _{k=1}^{n}|d_{k}|^{2}p(\lambda _{k})^{*}p(\lambda _{k})}}={\frac {\sum _{k=1}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})\left|p(\lambda _{k})\right|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}\left|p(\lambda _{k})\right|^{2}}}.$

Ключевое различие между числителем и знаменателем здесь заключается в том, что $k=1$ член исчезает в числителе, но не в знаменателе. Таким образом, если можно выбрать $p$ быть большим в $\lambda _{1}$ но мал при всех других собственных значениях, можно получить жесткую границу ошибки $\lambda _{1}-\theta _{1}$.

С $A$ имеет гораздо больше собственных значений, чем $p$ имеет коэффициенты, это может показаться сложной задачей, но один из способов выполнить ее — использовать полиномы Чебышева . Письмо $c_{k}$ на степень $k$ Полином Чебышева первого рода (тот, который удовлетворяет условию $c_{k}(\cos x)=\cos(kx)$ для всех $x$), у нас есть полином, который остается в диапазоне $[-1,1]$ на известном интервале $[-1,1]$ но быстро растет за его пределами. При некотором масштабировании аргумента мы можем заставить его отображать все собственные значения, кроме $\lambda _{1}$ в $[-1,1]$. Позволять

$p(x)=c_{m-1}\left({\frac {2x-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)$

(в случае $\lambda _{2}=\lambda _{1}$, используйте вместо этого наибольшее собственное значение, строго меньшее $\lambda _{1}$), то максимальное значение $|p(\lambda _{k})|^{2}$ для $k\geqslant 2$ является $1$ и минимальное значение $0$, так

$\lambda _{1}-\theta _{1}\leqslant \lambda _{1}-r(p(A)v_{1})={\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})|p(\lambda _{k})|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}|p(\lambda _{k})|^{2}}}\leqslant {\frac {\sum _{k=2}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})}{|d_{1}|^{2}|p(\lambda _{1})|^{2}}}\leqslant {\frac {(\lambda _{1}-\lambda _{n})\sum _{k=2}^{n}|d_{k}|^{2}}{|p(\lambda _{1})|^{2}|d_{1}|^{2}}}.$

Более того

$p(\lambda _{1})=c_{m-1}\left({\frac {2\lambda _{1}-\lambda _{2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)=c_{m-1}\left(2{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}+1\right);$

количество

$\rho ={\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n}}}$

(т. е. отношение первой собственной щели к диаметру остальной части спектра ) , таким образом, имеет здесь ключевое значение для скорости сходимости. Также пишу

$R=e^{\operatorname {arcosh} (1+2\rho )}=1+2\rho +2{\sqrt {\rho ^{2}+\rho }},$

мы можем заключить, что

${\begin{aligned}\lambda _{1}-\theta _{1}&\leqslant {\frac {(\lambda _{1}-\lambda _{n})\left(1-|d_{1}|^{2}\right)}{c_{m-1}(2\rho +1)^{2}|d_{1}|^{2}}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {1}{\cosh ^{2}((m-1)\operatorname {arcosh} (1+2\rho ))}}\\[6pt]&={\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {4}{\left(R^{m-1}+R^{-(m-1)}\right)^{2}}}\\[6pt]&\leqslant 4{\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n})R^{-2(m-1)}\end{aligned}}$

Таким образом, скорость сходимости контролируется главным образом $R$, поскольку эта граница уменьшается в раз $R^{-2}$ за каждую дополнительную итерацию.

Для сравнения можно рассмотреть, как зависит скорость сходимости степенного метода от $\rho$, но поскольку степенной метод в первую очередь чувствителен к отношению между абсолютными значениями собственных значений, нам нужно $|\lambda _{n}|\leqslant |\lambda _{2}|$ для собственного зазора между $\lambda _{1}$ и $\lambda _{2}$ быть доминирующим. При этом ограничении случай, наиболее благоприятствующий степенному методу, состоит в том, что $\lambda _{n}=-\lambda _{2}$, так что учтите это. В конце степенного метода вектор итераций:

$u=(1-t^{2})^{1/2}z_{1}+tz_{2}\approx z_{1}+tz_{2},$

где каждая новая итерация эффективно умножает $z_{2}$-амплитуда $t$ к

${\frac {\lambda _{2}}{\lambda _{1}}}={\frac {\lambda _{2}}{\lambda _{2}+(\lambda _{1}-\lambda _{2})}}={\frac {1}{1+{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}}}}}={\frac {1}{1+2\rho }}.$

Тогда оценка наибольшего собственного значения равна

$u^{*}Au=(1-t^{2})\lambda _{1}+t^{2}\lambda _{2},$

поэтому приведенную выше оценку скорости сходимости алгоритма Ланцоша следует сравнить с

$\lambda _{1}-u^{*}Au=(\lambda _{1}-\lambda _{2})t^{2},$

который уменьшается в раз $(1+2\rho )^{-2}$ для каждой итерации. Таким образом, разница сводится к тому, что между $1+2\rho$ и $R=1+2\rho +2{\sqrt {\rho ^{2}+\rho }}$. В $\rho \gg 1$ регион, последнее больше похоже на $1+4\rho$, и работает так же, как степенной метод с собственным зазором в два раза больше; заметное улучшение. Однако более сложным случаем является случай $\rho \ll 1,$ в котором $R\approx 1+2{\sqrt {\rho }}$ это еще большее улучшение собственного зазора; тот $\rho \gg 1$ Область — это область, в которой алгоритм Ланцоша с точки зрения сходимости дает наименьшее улучшение по сравнению со степенным методом.

Численная стабильность

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

Для алгоритма Ланцоша можно доказать, что при точной арифметике набор векторов $v_{1},v_{2},\cdots ,v_{m+1}$ строит ортонормированный базис, а решенные собственные значения/векторы являются хорошими аппроксимациями исходной матрицы. Однако на практике (поскольку вычисления выполняются в арифметике с плавающей запятой, где неточность неизбежна) ортогональность быстро теряется и в некоторых случаях новый вектор может даже линейно зависеть от уже построенного множества. В результате некоторые собственные значения результирующей трехдиагональной матрицы могут не быть аппроксимациями исходной матрицы. Следовательно, алгоритм Ланцоша не очень стабилен.

Пользователи этого алгоритма должны иметь возможность находить и удалять эти «ложные» собственные значения. Практическая реализация алгоритма Ланцоша направлена ​​на решение этой проблемы стабильности в трех направлениях:

1. Предотвратить потерю ортогональности,
2. Восстановите ортогональность после создания базиса.
3. После того, как все хорошие и «ложные» собственные значения идентифицированы, удалите ложные.

Вариации

Существуют варианты алгоритма Ланцоша, в которых задействованные векторы представляют собой высокие узкие матрицы вместо векторов, а нормализующие константы представляют собой небольшие квадратные матрицы. Они называются «блочными» алгоритмами Ланцоша и могут работать намного быстрее на компьютерах с большим количеством регистров и длительным временем выборки из памяти.

Многие реализации алгоритма Ланцоша перезапускаются после определенного количества итераций. Одним из наиболее влиятельных вариантов перезапуска является неявно перезапускаемый метод Ланцоша. который реализован в ARPACK . Это привело к ряду других возобновленных вариантов, таких как возобновление бидиагонализации Ланцо. Еще одним успешным вариантом перезапуска является метод Ланцоша с толстым перезапуском. который был реализован в пакете программного обеспечения под названием TRLan.

Нулевое пространство над конечным полем

Основная статья: Блок-алгоритм Ланцоша

В 1995 году Питер Монтгомери опубликовал алгоритм, основанный на алгоритме Ланцоша, для поиска элементов нулевого пространства большой разреженной матрицы над GF(2) ; поскольку группа людей, интересующихся большими разреженными матрицами над конечными полями, и группа людей, интересующихся большими проблемами собственных значений, почти не перекрываются, это часто также называют блочным алгоритмом Ланцоша, не вызывая необоснованной путаницы.

Приложения

Алгоритмы Ланцоша очень привлекательны, поскольку умножение на $A\,$ является единственной крупномасштабной линейной операцией. Поскольку механизмы взвешенного поиска текста реализуют именно эту операцию, алгоритм Ланцоша можно эффективно применять к текстовым документам (см. скрытое семантическое индексирование ). Собственные векторы также важны для крупномасштабных методов ранжирования, таких как алгоритм HITS, разработанный Джоном Кляйнбергом , или алгоритм PageRank, используемый Google.

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

Реализации

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

MATLAB и GNU Octave поставляются со ARPACK встроенным . Как хранимые, так и неявные матрицы можно анализировать с помощью функции eigs() ( Matlab / Octave ).

Аналогично, в Python пакет SciPy ARPACK имеет scipy.sparse.linalg.eigsh , который также является оболочкой для функций SSEUPD и DSEUPD из , которые используют неявно перезапущенный метод Ланцоша.

Реализация алгоритма Ланцоша в Matlab (обратите внимание на проблемы с точностью) доступна как часть пакета Matlab для распространения гауссовского доверия . ГрафЛаб ​ Библиотека совместной фильтрации включает крупномасштабную параллельную реализацию алгоритма Ланцоша (на C++) для многоядерных вычислений.

Библиотека PRIMME также реализует алгоритм, подобный Ланцошу.

Дальнейшее чтение

• Голуб, Джин Х .; Ван Лоан, Чарльз Ф. (1996). «Методы Ланцоша» . Матричные вычисления . Балтимор: Издательство Университета Джонса Хопкинса. стр. 470–507. ISBN  0-8018-5414-8 .
• Нг, Эндрю Ю .; Чжэн, Алиса X.; Джордан, Майкл И. (2001). «Анализ связей, собственные векторы и стабильность» (PDF) . IJCAI'01 Материалы 17-й Международной совместной конференции по искусственному интеллекту . 2 : 903–910.
• Эрик Кох (2019). «Точная диагонализация и метод Ланцоша» (PDF) . У Э. Паварини; Э. Кох; С. Чжан (ред.). Многочастичные методы для реальных материалов . Юлих. ISBN  978-3-95806-400-3 .