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

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

Основные авторы описания: Д.Р.Слюсарь (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.4, 2.7), М.А.Григорьев (1.2, 1.7, 1.8, 1.9, 2.4, 2.7, 3)

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

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

Алгоритм Ланцоша поиска собственных значений был опубликован Корнелием Ланцошем в 1950 году . Этот итерационный алгоритм применим только к эрмитовым матрицам [math]A[/math]. Метод позволяет за [math]k[/math] итераций вычислять [math]k[/math]-ое приближение собственных значений и собственных векторов исходной матрицы [math]A[/math].

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

Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный алгоритм Ланцоша с полной переортогонализацией предложенный в 1970 Ojalvo и Newman.

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

На вход алгоритма подается эрмитова матрица [math]A = A^\dagger[/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]

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

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

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы [math]T_k[/math], с помощью которых и будет найдены искомые собственные векторы исходной матрицы [math]A[/math].

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы [math]A[/math] на вектор [math]q_i[/math] с предыдущей итерации

[math]z = Aq_i[/math]

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

• Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
  • умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
  • скалярное произведение векторов;
  • линейная комбинация векторов (сложение/умножение на вещественные числа);
  • вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
• Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

Вычисление собственных значений полученной матрицы на практике выполняется с помощью QR-алгоритма и выходит за рамки описанного алгоритма.

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

Исходные данные: симметричная матрица [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 = & 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 \\ & \; \; \; \; exit\\ & else \\ & \; \; \; \; q_{i+1} = z / \beta_i \\ end \; & for \end{align} [/math]

После этого вычисляются собственные значения и собственные вектора симметричной трехдиагональной матрицы [math]T_k[/math] наиболее удобным образом.

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:

• Умножение квадратной матрицы [math]n * n[/math] на вектор длины [math]n[/math]. Требует [math]n * n[/math] умножений и сложений;
• Скалярное произведение векторов длины [math]n[/math]. Требует [math]n[/math] умножений и сложений;
• Сложение векторов длины [math]n[/math]. Требует [math]n[/math] сложений;
• Умножение вектора длины [math]n[/math] на скаляр. Требует [math]n[/math] умножений;
• Нахождение квадратичной нормы вектора длины [math]n[/math]. Требует [math]n[/math] умножений и сложений + извлечения квадратного корня;
• Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math] k \times k [/math]. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует [math]~O(k^{2.3})[/math] операций.

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

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

Таким образом, в худшем случае, алгоритм Ланцоша имеет сложность [math] O(k * n^2) [/math] и односится к квадратичным алгоритмам.

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

Информационный граф алгоритма можно разбить на две части:

• Граф алгоритма с отображением входных и выходных данных.
• Граф линейного оператора.

Рисунок 1. Информационный граф алгоритма.
[math]\|x\|[/math] — вычисление нормы,
[math]/[/math] — операция деления,
[math]\mathsf{F}[/math] — вычисление линейной комбинации векторов,
[math]*[/math] — операция умножения,
[math]\mathsf{Ax}[/math] — линейный оператор (Рисунок 2),
▽ - условие выхода [math]\mathsf{(\beta_i = 0)}[/math]

Рисунок 2. Информационный граф линейного оператора [math]\mathsf{Ax}[/math].
[math]+[/math], [math]*[/math] - операции сложения и умножения.

Рисунок 3. Граф алгоритма с указанием входных и выходных данных.
[math]\|x\|[/math] — вычисление нормы,
[math]/[/math] — операция деления,
[math]\mathsf{F}[/math] — вычисление линейной комбинации векторов,
[math]*[/math] — операция умножения,
[math]\mathsf{Ax}[/math] — линейный оператор (Рисунок 2),
▽ - условие выхода [math]\mathsf{(\beta_i = 0)}[/math],
[math]\mathsf{A, b, q_0 = 0, \beta_0 = 0}[/math] — входные данные,
[math]\mathsf{\alpha_i, \beta_i}[/math] — выходной результат.

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

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

• Умножение [math]A^{n * n}[/math] на вектор длины [math]n[/math] требует [math]n[/math] ярусов умножений и сложенийж
  • При этом сложение элементов вектора длины [math]n[/math] можно выполнить за [math]\log(n)[/math] [2];
• Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
• Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

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

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

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

Объём входных данных: [math]n * (n + 1) + 1 [/math].

Выходные данные: вектор собственных значений [math]\Lambda[/math], матрица собственных векторов [math]E[/math].

Объём выходных данных: [math]k * (n + 1)[/math].

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

• Сложность последовательного алгоритма [math]O(k*n^2)[/math];
• По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью [math]O(k * \log n)[/math], при классификации по ширине ЯПФ его сложность [math]O(n^2)[/math];
  • Таким образом, отношение последовательной сложности к параллельной [math]\frac{kn^2}{k \log{n}}[/math];
• Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма [math]\frac{k*(2*n^2+8*n-1)+3*n}{n^2+2*k}[/math]. При [math]k[/math] много меньше [math]n[/math] вычислительная мощность ≈ [math] 2*k[/math];
• Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
• Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы [math]A[/math], находящихся на границе ее спектра (в [math]T_{j}[/math] в первую очередь появляются максимальные по модулю собственные значения);
• Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
• Нестабильность алгоритма (эффект ложной сходимости) присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.

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

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

Для исследования масшабируемости алгоритма была написана реализация[3] на языке C++ с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[4].

Были исследованы:

• время выполнения программы в зависимости от размера входных данных (матрицы);
• время выполнения в зависимости от размера параллельных нод.

Дополнительно было измерено время работы исходя из оптимизационных флагов компилятора GCC [5].

Параметры запуска алгоритма:

• размер матрицы от 20000 до 175000 с шагом 2500;
• количество процессоров от 8 до 128 с шагом 8.

Программа запускалась на суперкомьютере "Ломоносов" со следующими характеристиками:

• Компилятор GCC 5.2.1;
• Версия MPI 1.8.4;
• Сборка проводилась командой:

  mpic++ -std=c++0x -O2 lanczos-mpi.cpp -lm -static-libstdc++ -o lanczos-mpi

Полученная эффективность реализации варьируется в пределах от 2% (на маленьких входных данных) до 45% (на больших входных данных и максимальном числе нодов).

На рисунке 4 представлена эффективность программы при отсутствии флага оптимизации O2.

Рисунок 4. Эффективность параллельной реализации алгоритма Ланцоща на MPI.

На рисунке 5 представлена эффективность алгоритма с оптимизацией компилятора GCC (-O2).

Рисунок 5. Эффективность параллельной реализации алгоритма Ланцоща на MPI с использованием оптимизации компилятора -O2.

Как видно из графиков выше, средняя эффективность минимальна при малых входных данных (на размере матрицы в 10000). В среднем данный показатель держится между 9% и 14%.

При этом, оптимизация компилятора MPIC++ позволяет получить выигрыш в эффективности более чем в 3 раза (45,5% против 14,7%). Однако, как можно заметить из графиков, при оптимизации алгоритм ведет себя более нестабильно, скорее всего из-за значительной траты времени на работу с памятью (более планые участки свидетельствует о том, что необходимые данные нашлись в кеше).

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

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

1. The IETL Project

2. MatLab

3. ARPACK

4. Julia Math

Также существуют официальные реализации на других языках (например R). Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

Для проверки масшабируемости, алгоритм был реализован на языке C++ с применением MPI [6].