Главная
Новости
Строительство
Ремонт
Дизайн и интерьер
Полезные советы



















Яндекс.Метрика





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

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

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

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

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

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

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

K n = [ b A b A 2 b ⋯ A n − 1 b ] . {displaystyle K_{n}={egin{bmatrix}b&Ab&A^{2}b&cdots &A^{n-1}bend{bmatrix}}.}

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

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

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

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

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

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

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

Имя:*
E-Mail:
Комментарий: