Итерация Арнольди (Nmyjgenx Gjukl,;n)

Перейти к навигации Перейти к поиску

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

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

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

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

[править | править код]

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

Эта последовательность сходится к собственному вектору соответствующего собственного значения с максимальным модулем. Это наводит на мысль, что много вычислений тратится впустую, т.к. в итоге используется лишь конечный результат . Тогда давайте вместо этого составим так называемую матрицу Крылова:

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

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

[править | править код]

Итерация Арнольди использует стабилизированный процесс Грама-Шмидта для получения последовательности ортонормированных векторов , называемых векторами Арнольди, таких, что для каждого векторы являются базисом подпространства Крылова . Алгоритм выглядит следующим образом:

  • Начать с произвольного вектора q1 с нормой 1.
  • Повторить для k = 2, 3, …
    • for j = 1 ... k − 1

Цикл по проецирует компоненту на Это обеспечивает ортогональность всех построенных векторов.

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

Каждый шаг цикла по производит одно умножение матрицы на вектор и около операций с дробными числами.

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