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

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

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

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

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

[math]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].[/math]

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

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

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

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

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

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

Выходные данные: трёхдиагональная симметричная вещественная матрица [math]T_k = \left[ \begin{array}{cccc} \alpha_1 & \beta_1 & & \\ \beta_1 & \ddots & \ddots & \\ & \ddots & \ddots & \beta_{k - 1}\\ & & \beta_{k - 1} & \alpha_k \end{array} \right][/math]

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

[math] \begin{array}{l} q_1 = b / \Vert b \Vert_2\\ 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} [/math]

Собственные значения и собственные векторы матрицы [math]T_k[/math] ищутся любым прямым алгоритмом.

Вычисление [math]z_j', z_j''[/math] — это полная переортогонализация [math]z_j[/math] методом Грама — Шмидта. Двойной запуск практически гарантирует, что [math]z_j''[/math] будет ортогонален [math]q_1, \dots, q_j[/math]. Заметим, что [math]\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z[/math], где [math]Q_j = [q_1, \dots, q_j][/math], поэтому переортогонализацию можно переписать в виде [math]z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'[/math].

Если [math]\beta_j = 0[/math] для какого-либо [math]j \lt k[/math], то [math]\Vert T_{ju}\Vert_2 = 0[/math], а значит, собственные значения [math]T_j[/math] в точности совпадают с каким-то собственными значениями [math]A[/math]. В этом случае дальнейшие вычисления прекращаются и либо используется полученная [math]T_j[/math] размерности меньшей, чем [math]k[/math], либо процедура запускается заново с другим начальным вектором [math]b[/math].

1.3 Вычислительное ядро алгоритма

У алгоритма построения матрицы [math]T_k[/math] в том виде, как он описан выше, можно выделить два ядра:

1. умножение матрицы [math]A[/math] на вектор [math]q_j[/math]: [math]z_j = A q_j[/math];
2. переортогонализация: [math]z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'[/math]

1.4 Макроструктура алгоритма

Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:

• арифметические операции над матрицами и векторами (сложение и умножение на число);
• умножение матрицы на вектор и скалярное произведение векторов;
• вычисление нормы вектора;
• добавление к матрице столбца.

1.5 Схема реализации последовательного алгоритма

Если переортогонализацию производить на каждом шаге строго по приведённым выше формулам, то получается следующий алгоритм:

[math] \begin{array}{l} \mathtt{input:}\; A, b, k\\ q_1 = b / \Vert b \Vert_2\\ Q = q_1\\ \mathtt{for}\; j = \overline{1, k}:\\ \quad z = A q_j \\ \quad \alpha_j = q_j^T z \\ \quad z = z - Q \left( Q^T z \right) \\ \quad z = z - Q \left( Q^T z \right)\\ \quad \beta_j = \Vert z \Vert_2\\ \quad \mathtt{if }\; \beta_j = 0\\ \quad \quad \mathtt{exit}\\ \quad \mathtt{end \; if}\\ \quad q_{j+1} = z/\beta_j\\ \quad Q = [Q, q_{j + 1}]\\ \mathtt{end \; for}\\ \mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k]. \end{array} [/math]