Русская Википедия:Итерация Арнольди
Шаблон:Нет источников В численной линейной алгебре итерация Арнольди является алгоритмом вычисления собственных значений. Арнольди находит приближение собственных значений и собственных векторов матриц общего вида(возможно не эрмитовой) с помощью построения ортонормированного базиса подпространства Крылова.
Метод Арнольди относится к алгоритмам линейной алгебры, которые позволяют получить частичное решение после небольшого количества итераций, в отличие от так называемых прямых методов, которые должны полностью завершиться для получения каких-либо удовлетворительных результатов(например отражения Хаусхолдера).
Если алгоритм применяется на эрмитовых матрицах, то он сводится к алгоритму Ланцоша. Итерация Арнольди была придумана Уолтером Эдвином Арнольди в 1951 г.
Подпространство Крылова и степенной метод
Итуитивным методом нахождения наибольшего (по модулю) собственного значения данной матрицы <math>A</math> размерами <math>m \times m</math> является степенной метод: начать с произвольного начального вектора <math>b</math>, вычислить <math>Ab, A^2b, A^3b, \dots</math>, нормирую результат после каждого вычисления.
Эта последовательность сходится к собственному вектору соответствующего собственного значения <math>\lambda_1</math> с максимальным модулем. Это наводит на мысль, что много вычислений тратится впустую, т.к. в итоге используется лишь конечный результат <math>A^{n - 1}b</math>. Тогда давайте вместо этого составим так называемую матрицу Крылова:
<math>K_{n} = \begin{bmatrix}b & Ab & A^{2}b & \cdots & A^{n-1}b \end{bmatrix}.</math>
Столбцы этой матрицы в общем случае не являются ортогональными, но мы можем получить из них ортогональный базис с помощью ортогонализации Грама-Шмидта. Полученное множество векторов будет являться ортогональным базисом подпространства Крылова <math>\mathcal{K}_{n}</math>. Можно ожидать, что вектора этого базиса будут хорошим приближением к векторам, соответствующим <math>n</math> наибольшим по модулю собственным значениям.
Итерация Арнольди
Итерация Арнольди использует стабилизированный процесс Грама-Шмидта для получения последовательности ортонормированных векторов <math>q_1, q_2, q_3, \dots</math>, называемых векторами Арнольди, таких, что для каждого <math>n</math> векторы <math>q_1, \dots, q_n</math> являются базисом подпространства Крылова <math>\mathcal{K}_{n}</math>. Алгоритм выглядит следующим образом:
- Начать с произвольного вектора q1 с нормой 1.
- Повторить для k = 2, 3, …
- <math> q_k \leftarrow Aq_{k-1} \,</math>
- for j = 1 ... k − 1
- <math> h_{j,k-1} \leftarrow q_j^* q_k \, </math>
- <math> q_k \leftarrow q_k - h_{j,k-1} q_j \, </math>
- <math> h_{k,k-1} \leftarrow \|q_k\| \, </math>
- <math> q_k \leftarrow \frac{q_k}{h_{k,k-1}} \, </math>
Цикл по <math>j</math> проецирует компоненту <math>q_k</math> на <math>q_1, \dots, q_{k - 1}.</math> Это обеспечивает ортогональность всех построенных векторов.
Алгоритм останавливается, когда <math>q_k</math> является нулевым вектором. Это происходит, когда минимальный многочлен матрицы <math>A</math> будет степени <math>k.</math>
Каждый шаг цикла по <math>k</math> производит одно умножение матрицы на вектор и около <math>4mk</math> операций с дробными числами.
На языке программирования python:
import numpy as np
def arnoldi_iteration(A, b, n: int):
"""Computes a basis of the (n + 1)-Krylov subspace of A: the space
spanned by {b, Ab, ..., A^n b}.
Arguments
A: m × m array
b: initial vector (length m)
n: dimension of Krylov subspace, must be >= 1
Returns
Q: m x (n + 1) array, the columns are an orthonormal basis of the
Krylov subspace.
h: (n + 1) x n array, A on basis Q. It is upper Hessenberg.
"""
m = A.shape[0]
h = np.zeros((n + 1, n))
Q = np.zeros((m, n + 1))
q = b / np.linalg.norm(b) # Normalize the input vector
Q[:, 0] = q # Use it as the first Krylov vector
for k in range(n):
v = A.dot(q) # Generate a new candidate vector
for j in range(k + 1): # Subtract the projections on previous vectors
h[j, k] = np.dot(Q[:, j].conj(), v)
v = v - h[j, k] * Q[:, j]
h[k + 1, k] = np.linalg.norm(v)
eps = 1e-12 # If v is shorter than this threshold it is the zero vector
if h[k + 1, k] > eps: # Add the produced vector to the list, unless
q = v / h[k + 1, k] # the zero vector is produced.
Q[:, k + 1] = q
else: # If that happens, stop iterating.
return Q, h
return Q, h
