Алгоритм Ланцоша для точной арифметики (без переортогонализации)

Алгоритм Ланцоша без переортогонализации
Последовательный алгоритм
Последовательная сложность	[math]O(k\; n^2)[/math]
Объём входных данных	[math]\frac{n(n + 1)}{2}[/math]
Объём выходных данных	[math]k\; (n + 1)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(n)[/math]
Ширина ярусно-параллельной формы	[math]O(n^2)[/math]

Основные авторы описания: А.Ю.Заспа, (раздел 2.2)

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

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

Алгоритм Ланцоша был опубликовн физиком и математиком Корнелием Ланцощем в 1950 году. Этот метод является частным случаем алгоритма Арнольда в случае, если исходная матрица [math]A[/math] - симметрична, и был представлен как итерационный метод вычисления собственных значений симметричной матрицы. Этот метод позволяет за [math]k[/math] итераций вычислять [math]k[/math] приближений собственных значений и собственных векторов исходной матрицы. Хотя алгоритм и был эффективным в вычислительном смысле, но он на некоторое время был предан забвению из-за численной неустойчивости. Только в 1970 Ojalvo и Newman модифицировали алгоритм для использования в арифметике с плавающей точкой. Новый метод получил название алгоритма Ланцоша с полной переортогонализацией. Но эта статья про его исходную версию. На вход алгоритма подается вещественная симметричная матрица [math]A = A^{T}[/math],

[math] A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} [/math]

Поэтому достаточно хранить только чуть больше половины элементов исходной матрицы.

Сам алгоритм соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. На каждой итерации строится матрица [math]Q_k = [q_1, q_2, \dots, q_k][/math] размерности [math]n \times k[/math], состоящая из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]T_k = Q^T_k A Q[/math] размерности [math]k \times k[/math].

[math] T_k = \begin{pmatrix} \alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \beta_2 \\ & \beta_2 & \ddots & \ddots \\ & & \ddots & \ddots & \beta_{k-1} \\ & & & \beta_{k-1} & \alpha_k \end{pmatrix} [/math]

Нахождение собственных значений и собственных векторов такой матрицы значительно проще чем их вычисление для исходной матрицы. Например, они могут быть вычислены с помощью метода «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы. В этом методе эта процедура занимает [math]O(k^3)[/math] операций, где константа оказывается на практике довольно мала.

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

Исходные данные: симметрическая матрица [math]A[/math], начальный вектор [math]q_1[/math].

Вычисляемые данные: собственные вектора матрицы [math]T_k[/math] являющиеся столбцами матрицы [math]Q_k V[/math], и матрица собственных значений [math]\Lambda[/math], где [math]V, \Lambda[/math] из спектрального разложения [math]T_k = V\Lambda V^T[/math].

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

[math] \begin{align} l_{11} & = \sqrt{a_{11}}, \\ l_{j1} & = \frac{a_{j1}}{l_{11}}, \quad j \in [2, n], \\ l_{ii} & = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}, \quad i \in [2, n], \\ l_{ji} & = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}, \quad i \in [2, n - 1], j \in [i + 1, n]. \end{align} [/math]

Существует также блочная версия метода, однако в данном описании разобран только точечный метод.

В ряде реализаций деление на диагональный элемент выполняется в два этапа: вычисление [math]\frac{1}{l_{ii}}[/math] и затем умножение на него всех (видоизменённых) [math]a_{ji}[/math] . Здесь мы этот вариант алгоритма не рассматриваем. Заметим только, что он имеет худшие параллельные характеристики, чем представленный.