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

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

Автор — Козлов Владимир, 617 гр.

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]q_1, \dots, q_k[/math]: [math]T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k][/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] ([math]O(n^2)[/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] ([math]O(nk^2)[/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]

Следует заметить, что очень часто в макрооперациях возникает операция вычисления суммы всех элементов массива. Её необходимо выполнять и при умножении матрицы на вектор, и при вычислении скалярного произведения, и при вычислении нормы вектора. Разные методы суммирования дают разный ресурс параллелизма. В этой работе предполагается, что везде для суммирования используется последовательный способ.

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

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

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

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

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

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

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

Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица [math]A[/math] и вектор [math]b[/math]), Out — выходные данные ([math]\alpha_j, \beta_j[/math] и [math]q_j[/math]), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.

Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица [math]A[/math]), Out — выходные данные ([math]\alpha_j, \beta_j[/math] и [math]q_j[/math]), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.

Информационный граф для алгоритма умножения матрицы на вектор можно видеть по этой ссылке. Информационный граф для алгоритма вычисления скалярного произведения можно видеть по этой ссылке. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из [math]n[/math] изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из [math]n[/math] вершин.

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

Для реализации [math]j[/math]-ой итерации алгоритма нужно выполнить следующие ярусы:

1. [math]n[/math] ярусов сложений с [math]n[/math] операциями умножения в каждом (вычисление [math]z[/math]);
2. [math]n[/math] ярусов сложений с [math]1[/math] операцией умножения в каждом (вычисление [math]\alpha_j[/math]);
3. [math]j[/math] ярусов сложений с [math]n[/math] операциями умножения в каждом, [math]n[/math] ярусов сложений с [math]j[/math] операциями умножения в каждом, [math]1[/math] ярус вычитаний размера [math]n[/math] (первая переортогонализация);
4. [math]j[/math] ярусов сложений с [math]n[/math] операциями умножения в каждом, [math]n[/math] ярусов сложений с [math]j[/math] операциями умножения в каждом, [math]1[/math] ярус вычитаний размера [math]n[/math] (вторая переортогонализация);
5. [math]n[/math] ярусов сложений с [math]1[/math] операцией умножения в каждом (вычисление [math]\beta_j[/math]);
6. [math]1[/math] ярус делений размера [math]n[/math] (вычисление [math]q_{j + 1}[/math]).

Итого ширина ЯПФ одной итерации составляет [math]O(n)[/math], высота — [math]4n + 2j + 3[/math]. Ширина ЯПФ всего алгоритма составляет [math]O(n)[/math], высота — [math]O(nk + k^2)[/math].

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

Входные данные: вещественная симметричная матрица [math]A[/math], вектор [math]b[/math], число [math]k[/math].

Размер входных данных: [math]n^2 + n[/math].

Выходные данные: набор чисел [math]\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}[/math], набор векторов [math]q_1, \dots, q_k[/math].

Размер выходных данных: [math]2k - 1 + nk[/math].

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

• Алгоритм недетерминирован, так как он может досрочно завершить работу.
• Алгоритм устойчив благодаря операции переортогонализации.
• Соотношение последовательной и параллельной сложности алгоритма — [math]n[/math]. Распараллелирование позволяет сделать алгоритм из квадратичного по [math]n[/math] линейным по [math]n[/math].
• Вычислительная мощность — порядка [math]2k[/math].
• В алгоритме есть длинные дуги — на протяжение всех итераций необходимо хранить набор векторов [math]Q[/math].

2 Программная реализация алгоритма

2.7 Существующие реализации алгоритма

• The NAG Library от Numerical Algorithms Group;
• Библиотека ARPACK.

3 Литература

[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.