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

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

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

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

Алгоритм

Входной эрмитова матрица размера , и , необязательно, число итераций (по умолчанию, пусть ). ${\ displaystyle A}$${\ Displaystyle п \ раз п}$${\ displaystyle m}$${\ displaystyle m = n}$

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

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

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

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

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

   3. Пусть . ${\ displaystyle w_ {j} '= Av_ {j}}$
   4. Пусть . ${\ displaystyle \ alpha _ {j} = w_ {j} '^ {*} v_ {j}}$
   5. Пусть . ${\ displaystyle w_ {j} = w_ {j} '- \ alpha _ {j} v_ {j} - \ beta _ {j} v_ {j-1}}$
4. Позвольте быть матрица со столбцами . Пусть . ${\ displaystyle V}$${\ displaystyle v_ {1}, \ dots, v_ {m}}$${\ displaystyle 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}}}$

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

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

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

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

Приложение к задаче о собственных значениях

Алгоритм Ланцоша чаще всего упоминается в контексте поиска собственных значений и собственных векторов матрицы, но в то время как обычная диагонализация матрицы сделала бы собственные векторы и собственные значения очевидными при проверке, то же самое нельзя сказать о трехдиагонализации, выполненной методом Ланцоша. алгоритм; Для вычисления даже одного собственного значения или собственного вектора необходимы нетривиальные дополнительные шаги. Тем не менее, применение алгоритма Ланцоша часто является значительным шагом вперед в вычислении собственного разложения. Если является собственным значением , и если ( является собственным вектором ), то является соответствующим собственным вектором (поскольку ). Таким образом, алгоритм Ланцоша преобразует задачу разложения на собственные числа для в задачу разложения на собственные числа для . ${\ displaystyle \ lambda}$${\ displaystyle A}$${\ displaystyle Tx = \ lambda x}$${\ displaystyle x}$${\ displaystyle T}$${\ displaystyle y = Vx}$${\ displaystyle A}$${\ Displaystyle Ay = AVx = VTV ^ {*} Vx = VTIx = VTx = V (\ lambda x) = \ lambda Vx = \ lambda y}$${\ displaystyle A}$${\ displaystyle T}$

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

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

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

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

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

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

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

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

Более предусмотрительный метод силы

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

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

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

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

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

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

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

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

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

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

${\ Displaystyle 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}}$ .

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

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

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

Коэффициенты априори удовлетворяют ${\ displaystyle h_ {k, j}}$

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

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

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

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

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

Элементарно, если эрмитово, то ${\ displaystyle A}$

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

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

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

поскольку последний реален в силу того, что является нормой вектора. Для одного получает ${\ displaystyle k = j}$

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

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

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

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

матрица эрмитова. Это означает, что это также нижний уровень Гессенберга, поэтому на самом деле он должен быть трехъядерным. Поскольку его главная диагональ является эрмитовой, она действительна, и поскольку ее первая поддиагональ реальна по построению, то же самое верно и для ее первой наддиагонали. Следовательно, это действительная симметричная матрица - матрица спецификации алгоритма Ланцоша. ${\ displaystyle H}$${\ displaystyle H}$${\ displaystyle H}$${\ displaystyle T}$

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

Один из способов характеризации собственных векторов эрмитовой матрицы как стационарные точки в частном Рэлеи${\ displaystyle A}$

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

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

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

${\ displaystyle {\ begin {align} 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 {выровнено}}}$

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

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

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

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

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

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

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

${\ displaystyle 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).}$

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

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

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

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

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

Теория сходимости Кэниела – Пейджа

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

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

Подпространство Крылова размерности равно ${\ displaystyle m}$

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

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

${\ displaystyle 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}}}}$

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

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

и в более общем плане

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

для любого полинома . ${\ displaystyle q}$

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

${\ displaystyle \ 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} ^ {*} \ сумма _ {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}}}.}$

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

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

${\ displaystyle p (x) = c_ {m-1} \ left ({\ frac {2x- \ lambda _ {2} - \ lambda _ {n}} {\ lambda _ {2} - \ lambda _ {n }}}\верно)}$

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

${\ displaystyle \ 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}}}.}$

более того

${\ displaystyle 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);}$

количество

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

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

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

мы можем сделать вывод, что

${\ displaystyle {\ begin {align} \ 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 {выровнено}}}$

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

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

${\ displaystyle u = (1-t ^ {2}) ^ {1/2} z_ {1} + tz_ {2} \ приблизительно z_ {1} + tz_ {2},}$

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

${\ displaystyle {\ 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}}.}$

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

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

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

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

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

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

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

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

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

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

Вариации

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

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

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

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

Приложения

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

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

Реализации

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

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

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

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

Заметки

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

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

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