Участник:Yastrebovks/Алгоритм Ланцоша с полной переортогонализацией

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

Авторы: Алексейчук Н.Н 616, Ястребов К.С. 609. Авторы внесли равнозначный вклад в каждый из разделов данной статьи.

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

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

Алгоритм Ланцоша ищет собвственные значения и собственные векторы для симетричной матрицы A вещественных чисел. Является итерацонным алгоритмом. Алгоритм Ланцоша использует степенной метод ([math]\displaystyle{ b_{k+1} = \frac{Ab_k}{||Ab_k||} }[/math]) для поиска наибольших собственных значений и векторов матриц.

В отличие от прямых алгоритмов требует мешьше памяти и мощности, что является несомненным плюсом для больших матриц.

Несмотря на свою скорость работы и экономию памяти, сначала не был популярным алгоритмом из – за недостаточной вычислительной устойчивости. В 1970 году Ojalvo и Newman [1] показали способ сделать алгоритм достаточно устойчивым. В этой же статье алгоритм был применен к расчету инженерной конструкции с большим количеством узлов, которые подвергались динамической нагрузке.

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

Памятка: Степенной метод нахождения наибольшего собственного числа матрицы можно сформулировать в предельном виде: если [math]\displaystyle{ b_0 }[/math] – случайный вектор, и [math]\displaystyle{ b_n+1 = Ab_n }[/math], тогда для больших чисел n предел [math]\displaystyle{ x_n/||x_n|| }[/math] стремится к нормированному наибольшему собственному вектору.

Алгоритм Ланцоша комбинирует метод Ланцоша для нахождения крыловского подпространства и метод Релэя – Ритца.

Подпространство Крылова для степенного метода: [math]\displaystyle{ K_m(v,A) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1] }[/math]

В качестве входных данных для алгоритма Ланцоша подаются квадратная матрица размерности [math]\displaystyle{ n }[/math]X[math]\displaystyle{ n }[/math]: [math]\displaystyle{ A=A^T }[/math]; а так же вектор начального приближения [math]\displaystyle{ b }[/math].

Метод осуществляет поиск трехдиагональной симметричной матрицы [math]\displaystyle{ T_k=Q_k^TAQ_k }[/math].

Таким образом матрица [math]\displaystyle{ T_k=\begin{bmatrix} \alpha_1 & \beta_2 \\ \beta_2 & \alpha_2 & \beta_3 &\\ &. & . & .\\ &&\beta_{k-1} & \alpha_{k-1} & \beta_k\\ &&&\beta_k & \alpha_k \end{bmatrix} }[/math]

Причем собственные значения [math]\displaystyle{ T_k }[/math] таковы, что приближают собственные значения исходной матрицы [math]\displaystyle{ A }[/math]. То есть на каждом [math]\displaystyle{ k }[/math]-м шаге из ортонормированных векторов Ланцоша строится матрица [math]\displaystyle{ Q_k = [q_1,q_2,...,q_k] }[/math] и в качестве приближенных собственных значений матрицы [math]\displaystyle{ A }[/math] принимаются числа Ритца.

Из-за ошибок округления вектора [math]\displaystyle{ A^{k-1}x_1 }[/math], формирующие подпространство Крылова, становятся неортогональными. Чтобы решить данную прорблему, проводят переортогонализацию методом Грамма-Шмидта.

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

Вычислительным ядром данного алгоритма являются следующие шаги:

• 1. [math]\displaystyle{ Aq=( \sum\nolimits_{i=^n}a_{1i}q_i, \sum\nolimits_{i=2}^na_{2i}q_i, ..., \sum\nolimits_{i=1}^na_{ni}q_i) }[/math].

• 2. [math]\displaystyle{ z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i. }[/math]

В данном случае первая операция выполняет умножение симметричной матрицы [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n }[/math]X[math]\displaystyle{ n }[/math] на [math]\displaystyle{ n }[/math]-мерный вектор [math]\displaystyle{ q }[/math], вследствие чего вычислительная сложность выполнения заключается в [math]\displaystyle{ n^2 }[/math] умножений и [math]\displaystyle{ n^2-n }[/math] сложений. Второе действие является процессом ортогонализации Грама-Шмидта. В последнем действии вычисляются [math]\displaystyle{ k^2n+k(n+2) }[/math] умножений и [math]\displaystyle{ k^2n + k(n + 1) + 2 }[/math] операций сложения.

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

Основные операции:

• 1. умножение матрицы на вектор [math]\displaystyle{ b=Ax }[/math].

• 2. ортогонализация методом Грмма-Шмидта.

• 3. вычисление нового вектора [math]\displaystyle{ q_{i+1} = z/ \beta }[/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{ 1.z=Aq_j, }[/math]

Вычисляется диагональный элемент конечной матрицы [math]\displaystyle{ T_k }[/math],

[math]\displaystyle{ 2.\alpha_j=q_j^Tz, }[/math]

Далее дважды проводится полная переортогонализация по Грамму-Шмидту. Данные действия необходимы для гарантии ортогональности базиса.

[math]\displaystyle{ 3.z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i, }[/math]

[math]\displaystyle{ 4.z=z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i, }[/math]

[math]\displaystyle{ 5.z=z-\alpha_jq_j-\beta_{j}q_{j-1}, }[/math]

Далее вычисляются недостающие элементы матрицы [math]\displaystyle{ T_k }[/math].

[math]\displaystyle{ 6.\beta_{j+1}=||z||, }[/math]

Если [math]\displaystyle{ \beta_{j+1}=0 }[/math], то алгоритм завершается.

[math]\displaystyle{ 7.q_{j+1}=z/\beta_{j+1}. }[/math]

В конце каждой итерации получаем очередной вектор [math]\displaystyle{ q_{j+1} }[/math] базиса с единичной нормой, который ортогонален всем предыдущим векторам базиса,

и элементы [math]\displaystyle{ \alpha_j }[/math], [math]\displaystyle{ \beta_{j+1} }[/math] трехдиагональной матрицы [math]\displaystyle{ T_k }[/math] .

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

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

Количество операций складывается из количества операций для классического метода Ланцоша[2] и количества операций для переортогонализации методом грамма-Шмидта[3].

Итоговая сложность составляет [math]\displaystyle{ O(n^2)+O(nk^2) }[/math].

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

Троеточием в данном алгоритме обозначаются операции, выполняющиеся для всех элементов от [math]\displaystyle{ i=1..n }[/math].

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

Хотя алгоритм является итерационным, возможно распараллелить внутри каждой итерации умножение матрицы на вектор и процесс переортогонализации Грамма-Шмидта.

Процесс умножения матриц matvec можно распараллелить несколькими способами[4].

Как видно из графа, для выполнения j-й операции, необходимо выполнить следующие ярусы:

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

Если учесть весь ресурс параллелизма можно получить сложность [math]\displaystyle{ O(k) }[/math].

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

На вход принимается сама матрица [math]\displaystyle{ A \in R^{n \times n} }[/math], случайный вектор [math]\displaystyle{ b }[/math] (впрочем возможен вариант, при котором этот вектор генерируется случайным образом). Ввиду симметричности входной матрицы достаточно передать лишь ее верхнюю треугольную матрицу, что дает объем входных данных [math]\displaystyle{ \frac{n(n + 1)}{2} }[/math].

Выходными данными для [math]\displaystyle{ j }[/math]ой итерации являются вектор [math]\displaystyle{ Aq_j }[/math] и собственное число [math]\displaystyle{ \lambda }[/math]. Общий объем выходных данных [math]\displaystyle{ k(n+1) }[/math].

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

Если [math]\displaystyle{ A }[/math] эрмитова матрица, то алгоритм Ланцоша и Bi-Lanczos сходятся к одинаковым трехдиагональным матрицам ритца[5].

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

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

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

Масштабируемость алгоритма и его реализации

Исследование масштабируемости реализации алгоритма проводилось с использованием стандарта MPI на суперкомпьютере "Ломоносов".

Набор изменяемых параметров запуска реализации алгоритма и границы значений параметров алгоритма:

• Число процессоров [math]\displaystyle{ 1,2,4,8,16,32,64,96,128 }[/math]
• Размерность матрицы [math]\displaystyle{ 1000:15000,1700,21000,25000,30000 }[/math]
• Количество итераций [math]\displaystyle{ 500 }[/math]

Сборка осуществлялась с параметрами:

• openmpi/1.5.5-icc
• intel/13.1.0

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

The IETL Project http://www.comp-phys.org/software/ietl/ C++

NAG Library http://www.nag.com/content/nag-library C, C++, Fortran, C#, MATLAB, R

ARPACK https://people.sc.fsu.edu/~jburkardt/m_src/arpack/arpack.html MATLAB

GrapLab https://turi.com/products/create/open_source.html C++

С частичной переортаганализацией

LANSO/PLANSO http://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_PLANSO.html Fortran (уже распараллелена)

Литература

Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).