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

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

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

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

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

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

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

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

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

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

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

Выходные данные: трёхдиагональная симметричная вещественная матрица [math]\displaystyle{ 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]\displaystyle{ q_1, \dots, q_k }[/math]: [math]\displaystyle{ T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k] }[/math].

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

[math]\displaystyle{ \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]\displaystyle{ T_k }[/math] ищутся любым прямым алгоритмом.

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

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

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

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

1. умножение матрицы [math]\displaystyle{ A }[/math] на вектор [math]\displaystyle{ q_j }[/math]: [math]\displaystyle{ z_j = A q_j }[/math] ([math]\displaystyle{ O(n^2) }[/math] умножений вещественных чисел);
2. переортогонализация: [math]\displaystyle{ z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j' }[/math] ([math]\displaystyle{ O(nk^2) }[/math] умножений вещественных чисел).

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

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

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

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

[math]\displaystyle{ \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]

Последовательная сложность алгоритма

На [math]\displaystyle{ j }[/math]-ой итерации алгоритм выполняет

• [math]\displaystyle{ n^2 + n + 4nj + n + n = n^2 + 3n + 4nj }[/math] умножений и делений;
• [math]\displaystyle{ n(n - 1) + (n - 1) + (4n - 2)j + (n - 1) = n^2 + (n - 2) + (4n - 2)j }[/math] сложений и вычитаний;
• одно извлечение квадратного корня.

Всего за [math]\displaystyle{ k }[/math] итераций алгоритм выполняет

• [math]\displaystyle{ (n + 3)nk + 2nk(k + 1) }[/math] умножений и делений;
• [math]\displaystyle{ (n^2 + n - 2)k + (2n - 1)k(k + 1) }[/math] сложений и вычитаний.

Таким образом, алгоритм имеет последовательную сложность [math]\displaystyle{ O(n^2k + nk^2) }[/math].

Информационный граф

Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.