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

Алгоритм Ланцоша без переортогонализации
Последовательный алгоритм
Последовательная сложность	[math]O(kn^2)[/math]
Объём входных данных	[math]\frac{n(n + 1)}{2}[/math]
Объём выходных данных	[math]k(n + 1)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(k)[/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]b[/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} q_1 = & \frac{b}{\|b\|_2},\; \beta_0 = 0,\; q_0 = 0\\ for \; & i = 1 \; to \; k \\ & z = Aq_i\\ & \alpha_i = q^T_i z\\ & z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\ & \beta_i = \|z\|_2\\ & if \; \beta_i == 0 \; then \; break\\ & q_{i+1} = \frac{z}{\beta_i}\\ end \; & for \end{align} [/math]

Затем вычислить собственные значения и собственные вектора матрицы [math]T_k[/math]

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

Вычислительным ядром на каждой итерации является вычисление произведения матрицы на вектор:

[math]z = Aq_i[/math]

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

Макрооперациями алгоритма являются:

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

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

Последовательность исполнения метода следующая:

[math] \begin{align} 1. \, \beta_0 = & 0,\; q_0 = 0,\\ \|b\|_2 = &\sqrt{\sum\limits_{j=1}^{n} b_j^2},\\ 2. \, q_{1_{j}} = &\frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots \,, n\\ for \; i = & 1 \; to \; k \\ 3. \, z_j & = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\, \dots \,, n\\ 4. \, \alpha_i & = \sum\limits_{j=1}^{n}q_j z_j\\ 5. \, z_j & = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\, \dots \,, n\\ 6. \, \beta_i & = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_m^2}\\ if \; \beta_i & == 0 \; then \; break\\ 7. \, q_{i+1_j} & = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n\\ end \; & for \end{align} [/math]

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

Для вычисления [math]k[/math] собственных значений матрицы порядка [math]n[/math] и соответствующих им собственных векторов алгоритмом Ланцоша без переортогонализации в худшем случае требуется:

• [math] kn^2 + 4 kn + n [/math] умножений,
• [math] k + 1 [/math] вычислений квадратного корня,
• [math] n(k + 1) [/math] делений,
• [math] n + k(n^2+3n-1)[/math] сложений (вычитаний),
• нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math] k \times k [/math].

Умножения и сложения (вычитания) составляют основную часть алгоритма.

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

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

TBD

1.8 Ресурс параллелизма алгоритма

Алгоритм Ланцоша в параллельной форме состоит из инициализации, [math]k[/math] итераций и вычисления собственных значений матрицы [math]T_k[/math], рассмотрение которого выходит за рамки данной статьи.

Заметим, что вычисление суммы [math]n[/math] элементов имеет высоту [math]\log n[/math] и линейную ширину ярусов [math]\frac{n}{2}, \frac{n}{4}, ... , 1[/math].

Инициализации состоит из следующих операций:

• вычисление нормы:
  • вычисление скалярного произведения вектора самого на себя:
    • ярус умножений шириной [math]n[/math]
    • [math]\log{n}[/math] ярусов сложений с наибольшей шириной [math]\frac{n}{2}[/math]
  • ярус извлечения квадратного корня с единичным вычислением
• деление вектора на число (нормирование вектора):
  • ярус [math]n[/math] операций деления

Итерационная часть алгоритма состоит из следующих операций:

• умножение матрицы на число:
  • ярус умножений с шириной [math]n^2[/math]
  • [math]\log n[/math] ярусов с наибольшей шириной [math]\frac{n^2}{2}[/math] (блок вычисления n сумм n элементов)
• вычисление скалярного произведения двух векторов:
  • ярус умножений шириной [math]n[/math]
  • [math]\log{n}[/math] ярусов сложений с наибольшей шириной [math]\frac{n}{2}[/math]
• вычисление линейной комбинации векторов:
  • ярус умножений шириной [math]2n[/math]
  • два яруса сложений шириной [math]n[/math]
• вычисление нормы:
  • вычисление скалярного произведения вектора самого на себя:
    • ярус умножений шириной [math]n[/math]
    • [math]\log{n}[/math] ярусов сложений с максимальной шириной [math]\frac{n}{2}[/math]
  • ярус извлечения квадратного корня с единичным вычислением
• проверка на выход из цикла
• деление вектора на число (нормирование вектора, может отсутствовать на последней итерации):
  • ярус [math]n[/math] операций деления

Таким образом, в параллельном варианте основную долю времени будут занимать операции сложения при умножении матрицы на вектор. При классификации по высоте ЯПФ Алгоритм Ланцоша относится к алгоритмам со сложностью [math]O(k \log n)[/math]. При классификации по ширине ЯПФ его сложность будет [math]O(n^2)[/math].

1.9 Входные и выходные данные алгоритма

TBD

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

TBD