Участник:Kozlov Vladimir/Алгоритм Ланцоша для арифметики с плавающей точкой с полной переортогонализацией

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Алгоритм Ланцоша — это итерационный алгоритм поиска $k$ приближённых собственных значений симметричной вещественной матрицы $A$ размера $n \times n$. Алгоритм применяется, когда матрица $A$ слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм метод Рэлея — Ритца поиска приближённых собственных значений и метод Ланцоша построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска $k$ приближённых собственных значений симметричной вещественной матрицы $A$ размера $n \times n$. Если $Q = [Q_k, Q_u]$ — ортонормированная матрица размера $n \times n$, $Q_k$ имеет размер $n \times k$, $Q_u$ имеет размер $n \times n - k$, то можно записать равенство

$T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] = \left[ \begin{array}{cc} Q_k^T A Q_k & Q_k^T A Q_u\\ Q_u^T A Q_k & Q_u^T A Q_u \end{array} \right] = \left[ \begin{array}{cc} T_{k} & T_{ku}^T\\ T_{ku} & T_{u} \end{array} \right].$

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы $T_k = Q_k^T A Q_k$ объявляются приближёнными собственными значениями матрицы $A$. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если $T_k = V \Lambda V^{-1}$ — спектральное разложение $T_k$, то пара $(Q_k V, \Lambda)$ минимизирует функционал $L(P_k, D) = \Vert A P_k - P_k D \Vert_2$, причём $L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2$, то есть $A \approx (Q_k V) \Lambda (Q_k V)^{-1}$. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы $A$. Более того, можно показать, что собственные значения $T$ отличаются от некоторых собственных значений $A$ не более чем на $\Vert T_{ku} \Vert_2$.

Метод Ланцоша — это метод построения матрицы $Q$, при использовании которого, во-первых, матрица $T$ оказывается симметричной трёхдиагональной, во-вторых, столбцы $Q$ и $T$ вычисляются последовательно. Трёхдиагональность $T$ приводит к следующим явлениям:

1. матрица $T_k$ является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
2. матрица $T_{ku}$ имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца $q_{j + 1}$ матрицы $Q$ достаточно знать только $q_{j - 1}$ и $q_{j}$ в силу трёхдиагональности матрицы $T$. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов $q_{1}, \dots, q_{k}$ перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

1.2 Математическое описание алгоритма

Исходные данные: симметрическая матрица $A$ — матрица, для которой будут вычисляться собственные значения, вектор $b$ — начальное приближение для метода Ланцоша.

Формулы метода:

$\begin{array}{l} q_1 = b / \Vert b \Vert_2, \beta_0 = 0\\ j = \overline{1, k}:\\ \quad z_j = A q_j \\ \quad \alpha_j = q_j^T A q_j = q_j^T z_j \\ \quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\ \quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\ \quad \beta_j = \Vert z_j'' \Vert_2\\ \quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j \end{array}$

Вычисление $z_j', z_j''$ — это полная переортогонализация $z_j$ методом Грама — Шмидта. Двойной запуск практически гарантирует, что $z_j''$ будет ортогонален $q_1, \dots, q_j$.