Итерация Арнольди

Итерационный метод аппроксимации собственных векторов

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

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

Применительно к эрмитовым матрицам он сводится к алгоритму Ланцоша . Итерация Арнольди была изобретена У. Э. Арнольди в 1951 году. ^[1]

Подпространства Крылова и степенная итерация

Интуитивным методом нахождения наибольшего (по абсолютной величине) собственного значения заданной матрицы размера m × m является степенная итерация : начиная с произвольного начального вектора b , вычисляем Ab , A ² b , A ³ b , ..., нормализуя результат после каждого применения матрицы A. ${\displaystyle A}$

Эта последовательность сходится к собственному вектору, соответствующему собственному значению с наибольшим абсолютным значением, . Однако много потенциально полезных вычислений тратится впустую, если использовать только конечный результат, . Это говорит о том, что вместо этого мы формируем так называемую матрицу Крылова : ${\displaystyle \lambda _{1}}$ ${\displaystyle A^{n-1}b}$

${\displaystyle K_{n}={\begin{bmatrix}b&Ab&A^{2}b&\cdots &A^{n-1}b\end{bmatrix}}.}$

Столбцы этой матрицы в общем случае не ортогональны , но мы можем извлечь ортогональный базис , с помощью такого метода, как ортогонализация Грама–Шмидта . Таким образом, полученный набор векторов является ортогональным базисом подпространства Крылова , . Мы можем ожидать, что векторы этого базиса будут охватывать хорошие приближения собственных векторов, соответствующих наибольшим собственным значениям, по той же причине, по которой приближается доминирующий собственный вектор. ${\displaystyle {\mathcal {K}}_{n}}$ ${\displaystyle n}$ ${\displaystyle A^{n-1}b}$

Итерация Арнольди

Итерация Арнольди использует модифицированный процесс Грама–Шмидта для создания последовательности ортонормированных векторов q ₁ , q ₂ , q ₃ , ..., называемых векторами Арнольди , таких, что для каждого n векторы q ₁ , ..., q _n охватывают подпространство Крылова . Явно алгоритм выглядит следующим образом: ${\displaystyle {\mathcal {K}}_{n}}$

Начнем с произвольного вектора q 1 с нормой 1. Повторим для  k = 2, 3, ... q k  := A  q k −1  для  j  от 1 до  k − 1 h j , k −1  := q j *  q k  q k  := q k − h j , k −1  q j  h k , k −1  := ‖ q k ‖ q k  := q k / h k , k −1

j - контур проецирует компоненту в направлениях . Это обеспечивает ортогональность всех сгенерированных векторов. ${\displaystyle q_{k}}$ ${\displaystyle q_{1},\dots ,q_{k-1}}$

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

Каждый шаг k -цикла требует одного произведения матрицы на вектор и приблизительно 4 mk операций с плавающей точкой.

На языке программирования Python с поддержкой библиотеки NumPy :

импортировать  numpy  как  npdef  arnoldi_iteration ( A ,  b ,  n :  int ): """Вычислить базис (n + 1)-подпространства Крылова матрицы A.  Это пространство, охватываемое векторами {b, Ab, ..., A^nb}. Параметры  ----------  A : array_like  Массив размером m × m.  b : array_like  Начальный вектор (длина m).  n : int  На единицу меньше размерности подпространства Крылова или, что эквивалентно, *степени* пространства Крылова. Должно быть >= 1.  Возвращает  -------  Q : numpy.array  Массив размером mx (n + 1), где столбцы являются ортонормированным базисом подпространства Крылова.  h : numpy.array  Массив размером (n + 1) xn. A на базисе Q. Это верхний Хессенберг.  """ eps = 1e-12 h = np . zeros (( n + 1 , n )) Q = np . zeros (( A . shape [ 0 ], n + 1 )) # Нормализуем входной вектор Q [:, 0 ] = b / np . linalg . norm ( b , 2 ) # Используем его как первый вектор Крылова для k в диапазоне ( 1 , n + 1 ): v = np . dot ( A , Q [:, k - 1 ]) # Генерируем новый вектор-кандидат для j в диапазоне ( k ): # Вычитаем проекции на предыдущие векторы h [ j , k - 1 ] = np . dot ( Q [:, j ] . conj (), v ) v = v - h [ j , k - 1 ] * Q [:, j ] h [ k , k - 1 ] = np . линалг.норма ( v , 2 ) если h [ k , k - 1​​                                                                            ]  >  eps :  # Добавить полученный вектор в список, если только  Q [:,  k ]  =  v  /  h [ k ,  k  -  1 ]  иначе :  # Если это произошло, прекратить итерацию.  return  Q ,  h  return  Q ,  h

Свойства итерации Арнольди

Пусть Q _n обозначает матрицу размером m на n, образованную первыми n векторами Арнольди q ₁ , q ₂ , ..., q _n , а H _n — (верхнюю матрицу Хессенберга ), образованную числами h _{j , k ,} вычисленными с помощью алгоритма:

${\displaystyle H_{n}=Q_{n}^{*}AQ_{n}.}$

Метод ортогонализации должен быть специально выбран таким образом, чтобы нижние компоненты Арнольди/Крылова были удалены из высших векторов Крылова. Как можно выразить в терминах q ₁ , ..., q _{i +1} по построению, они ортогональны q _{i +2} , ..., q _n , ${\displaystyle Aq_{i}}$

Тогда у нас есть

${\displaystyle H_{n}={\begin{bmatrix}h_{1,1}&h_{1,2}&h_{1,3}&\cdots &h_{1,n}\\h_{2,1}&h_{2,2}&h_{2,3}&\cdots &h_{2,n}\\0&h_{3,2}&h_{3,3}&\cdots &h_{3,n}\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\cdots &0&h_{n,n-1}&h_{n,n}\end{bmatrix}}.}$

Матрицу H _n можно рассматривать как A в подпространстве с векторами Арнольди в качестве ортогонального базиса; A ортогонально проецируется на . Матрицу H _n можно охарактеризовать следующим условием оптимальности. Характеристический многочлен H _n минимизирует || p ( A ) q ₁ || ₂ среди всех монических многочленов степени n . Эта задача оптимальности имеет единственное решение тогда и только тогда , когда итерация Арнольди не ломается. ${\displaystyle {\mathcal {K}}_{n}}$ ${\displaystyle {\mathcal {K}}_{n}}$

Соотношение между матрицами Q в последующих итерациях определяется выражением

${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}}$

где

${\displaystyle {\tilde {H}}_{n}={\begin{bmatrix}h_{1,1}&h_{1,2}&h_{1,3}&\cdots &h_{1,n}\\h_{2,1}&h_{2,2}&h_{2,3}&\cdots &h_{2,n}\\0&h_{3,2}&h_{3,3}&\cdots &h_{3,n}\\\vdots &\ddots &\ddots &\ddots &\vdots \\\vdots &&0&h_{n,n-1}&h_{n,n}\\0&\cdots &\cdots &0&h_{n+1,n}\end{bmatrix}}}$

представляет собой матрицу размером ( n +1) на n , образованную путем добавления дополнительной строки к H _n .

Нахождение собственных значений с помощью итерации Арнольди

Идея итерации Арнольди как алгоритма собственных значений заключается в вычислении собственных значений в подпространстве Крылова. Собственные значения H _n называются собственными значениями Ритца . Поскольку H _n является матрицей Хессенберга скромного размера, ее собственные значения можно эффективно вычислить, например, с помощью алгоритма QR или несколько связанного с ним алгоритма Фрэнсиса . Также сам алгоритм Фрэнсиса можно считать связанным с степенными итерациями, работающими на вложенном подпространстве Крылова. Фактически, самая базовая форма алгоритма Фрэнсиса, по-видимому, заключается в выборе b равным Ae ₁ и расширении n до полной размерности A . Улучшенные версии включают один или несколько сдвигов, и более высокие степени A могут быть применены за один шаг. ^[2]

Это пример метода Рэлея-Ритца .

На практике часто наблюдается, что некоторые собственные значения Ритца сходятся к собственным значениям A . Поскольку H _n имеет размер n на n , он имеет не более n собственных значений, и не все собственные значения A могут быть приближены. Обычно собственные значения Ритца сходятся к наибольшим собственным значениям A . Чтобы получить наименьшие собственные значения A , вместо этого следует использовать обратную (операцию) A . Это можно связать с характеристикой H _n как матрицы, характеристический полином которой минимизирует || p ( A ) q ₁ || следующим образом. Хороший способ сделать p ( A ) малым — выбрать полином p таким образом, чтобы p ( x ) было малым, когда x является собственным значением A . Следовательно, нули p (и, следовательно, собственные значения Ритца) будут близки к собственным значениям A .

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

Итерация Арнольди, демонстрирующая сходимость значений Ритца (красные) к собственным значениям (черные) матрицы 400x400, составленной из равномерных случайных значений в области [-0,5 +0,5]

Возобновленная итерация Арнольди

Из-за практических соображений хранения общие реализации методов Арнольди обычно перезапускаются после фиксированного числа итераций. Одним из подходов является метод неявного перезапуска Арнольди (Implicitly Restarted Arnoldi Method, IRAM) ^[3] Лехука и Соренсена, который был популяризирован в свободном и открытом программном пакете ARPACK . ^[4] Другим подходом является алгоритм Крылова-Шура от GW Stewart, который более стабилен и прост в реализации, чем IRAM. ^[5]

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

Обобщенный метод минимальных остатков (GMRES) — это метод решения уравнения Ax = b, основанный на итерации Арнольди.

Ссылки

1. Арнольди, У. Э. (1951). «Принцип минимизированных итераций в решении проблемы собственных значений матрицы». Quarterly of Applied Mathematics . 9 (1): 17–29. doi : 10.1090/qam/42792 . ISSN  0033-569X.
2. Дэвид С. Уоткинс. Алгоритм Фрэнсиса. Университет штата Вашингтон. Получено 14 декабря 2022 г.
3. RB Lehoucq & DC Sorensen (1996). «Методы дефляции для неявно перезапущенной итерации Арнольди». Журнал SIAM по матричному анализу и приложениям . 17 (4): 789–821. doi :10.1137/S0895479895281484. hdl : 1911/101832 .
4. RB Lehoucq; DC Sorensen & C. Yang (1998). "ARPACK Users Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods". SIAM. Архивировано из оригинала 2007-06-26 . Получено 2007-06-30 .
5. Стюарт, GW (2002). «Алгоритм Крылова--Шура для больших собственных задач». Журнал SIAM по матричному анализу и приложениям . 23 (3): 601–614. doi :10.1137/S0895479800371529. ISSN  0895-4798.

• У. Э. Арнольди, «Принцип минимизации итераций при решении проблемы собственных значений матрицы», Quarterly of Applied Mathematics , том 9, страницы 17–29, 1951.
• Юсеф Саад , Численные методы решения задач на большие собственные значения , Manchester University Press, 1992. ISBN 0-7190-3386-1 . 
• Ллойд Н. Трефетен и Дэвид Бау, III, Численная линейная алгебра , Общество промышленной и прикладной математики, 1997. ISBN 0-89871-361-7 . 
• Яшке, Леонард: Предварительно обусловленные методы Арнольди для систем нелинейных уравнений . (2004). ISBN 2-84976-001-3 
• Реализация: Matlab поставляется со встроенным ARPACK. Как сохраненные, так и неявные матрицы могут быть проанализированы с помощью функции eigs().