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

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

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

Если алгоритм применяется на эрмитовых матрицах, то он сводится к алгоритму Ланцоша. Итерация Арнольди была придумана Уолтером Эдвином Арнольди в 1951 г.

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

Эта последовательность сходится к собственному вектору соответствующего собственного значения ${\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 q_{1},q_{2},q_{3},\dots }$, называемых векторами Арнольди, таких, что для каждого ${\displaystyle n}$ векторы ${\displaystyle q_{1},\dots ,q_{n}}$ являются базисом подпространства Крылова ${\displaystyle {\mathcal {K}}_{n}}$. Алгоритм выглядит следующим образом:

• Начать с произвольного вектора q₁ с нормой 1.
• Повторить для k = 2, 3, …
  • ${\displaystyle q_{k}\leftarrow Aq_{k-1}\,}$
  • for j = 1 ... k − 1
    • ${\displaystyle h_{j,k-1}\leftarrow q_{j}^{*}q_{k}\,}$
    • ${\displaystyle q_{k}\leftarrow q_{k}-h_{j,k-1}q_{j}\,}$
  • ${\displaystyle h_{k,k-1}\leftarrow \|q_{k}\|\,}$
  • ${\displaystyle q_{k}\leftarrow {\frac {q_{k}}{h_{k,k-1}}}\,}$

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

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

Каждый шаг цикла по ${\displaystyle k}$ производит одно умножение матрицы на вектор и около ${\displaystyle 4mk}$ операций с дробными числами.

На языке программирования 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