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

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

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

Применительно к эрмитовым матрицам он сводится к алгоритму Ланцоша . Итерация Арнольди была изобретена WE Arnoldi в 1951 году.

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

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

Эта последовательность сходится к собственному вектору, соответствующему собственному значению с наибольшим абсолютным значением: $\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$ наибольшие собственные значения по той же причине, что $A^{n-1}b$ аппроксимирует доминирующий собственный вектор.

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

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

Start with an arbitrary vector q1 with norm 1.
Repeat for k = 2, 3, ...
  qk := A qk−1
  for j from 1 to k − 1
    hj,k−1 :=  qj* qk
    qk := qk − hj,k−1 qj
  hk,k−1 := ‖qk‖
  qk := qk / hk,k−1

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

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

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

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

import numpy as np

def arnoldi_iteration(A, b, n: int):
    """Compute a basis of the (n + 1)-Krylov subspace of the matrix A.

    This is the space spanned by the vectors {b, Ab, ..., A^n b}.

    Parameters
    ----------
    A : array_like
        An m × m array.
    b : array_like
        Initial vector (length m).
    n : int
        One less than the dimension of the Krylov subspace, or equivalently the *degree* of the Krylov space. Must be >= 1.
    
    Returns
    -------
    Q : numpy.array
        An m x (n + 1) array, where the columns are an orthonormal basis of the Krylov subspace.
    h : numpy.array
        An (n + 1) x n array. A on basis Q. It is upper Hessenberg.
    """
    eps = 1e-12
    h = np.zeros((n + 1, n))
    Q = np.zeros((A.shape[0], n + 1))
    # Normalize the input vector
    Q[:, 0] = b / np.linalg.norm(b, 2)  # Use it as the first Krylov vector
    for k in range(1, n + 1):
        v = np.dot(A, Q[:, k - 1])  # Generate a new candidate vector
        for j in range(k):  # Subtract the projections on previous vectors
            h[j, k - 1] = np.dot(Q[:, j].conj(), v)
            v = v - h[j, k - 1] * Q[:, j]
        h[k, k - 1] = np.linalg.norm(v, 2)
        if h[k, k - 1] > eps:  # Add the produced vector to the list, unless
            Q[:, k] = v / h[k, k - 1]
        else:  # If that happens, stop iterating.
            return Q, h
    return Q, h

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

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

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

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

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

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

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

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

где

${\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 можно применять за один этап.

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

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

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

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

Перезапущена итерация Арнольди.

Из практических соображений хранения общие реализации методов Арнольди обычно перезапускаются после фиксированного количества итераций. Одним из подходов является метод неявного перезапуска Арнольди (IRAM). от Lehoucq и Sorensen, который был популяризирован в бесплатном пакете программного обеспечения с открытым исходным кодом ARPACK . Другой подход — алгоритм Крылова-Шура Г. В. Стюарта, который более стабилен и прост в реализации, чем IRAM.

См. также

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

Ссылки

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