Участник: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] в точности совпадают с каким-то собственными значениями <mathA</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 Макроструктура алгоритма

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

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

Можно предложить также и другой набор макрооперация, соответствующий использованию дополнительной памяти для выполнения переортогонализации:

• арифметические операции над матрицами и векторами;
• умножение матрицы на вектор и скалярное произведение векторов;
• вычисление нормы вектора.
• вычисление произведений вида [math]QQ^T[/math] — поскольку итоговая матрица является симметричной, можно такие произведения можно считать вдвое быстрее, чем обычные.

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]