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

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

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

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

Подпространство Крылова и степенной метод

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

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

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

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

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

Итерация Арнольди использует стабилизированный процесс Грама-Шмидта для получения последовательности ортонормированных векторов $q_{1},q_{2},q_{3},\dots$, называемых векторами Арнольди, таких, что для каждого $n$ векторы $q_{1},\dots ,q_{n}$ являются базисом подпространства Крылова ${\mathcal {K}}_{n}$. Алгоритм выглядит следующим образом:

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

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

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

Каждый шаг цикла по $k$ производит одно умножение матрицы на вектор и около $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