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

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Konshin и ASA.

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

Основные авторы описания: Заспа А.Ю. (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7), Фролов А.А. (1.6, 1.7, 1.8, 1.9, 1.10, 2.4, 2.7)

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

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

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

$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}$

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

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

$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}$

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

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

Исходные данные: симметрическая матрица $A$, начальный вектор $b$.

Вычисляемые данные: собственные вектора матрицы $T_k$ являющиеся столбцами матрицы $Q_k V$, и матрица собственных значений $\Lambda$, где $V, \Lambda$ из спектрального разложения $T_k = V\Lambda V^T$. До начала первой итерации нормируется вектор начального приближения $q_1 = \frac{b}{\|b\|_2}$. Затем задаются константы $\beta_0 = 0,\; q_0 = 0$.

После выполнения данных действий начинаются итерации алгоритма. Которых будет не больше чем $k$. В ходе итераций постепенно формируется матрица $T_k$ из $\alpha_i$-x и $\beta_i$-х.

На $i$-й итерации происходят следующие вычисления. Вначале вычисляется произведение матрицы на вектор $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$. Если $\beta_i$ оказалась равной нулю, то больше никаких итераций в алгоритме не производят, сразу переходят к вычислению собственных значений и собственных векторов. Иначе, в конце итерации считают $q_{i+1} = \frac{z}{\beta_i}$ и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы $T_k$.

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

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

$z = Aq_i$

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

Для удобства понимания алгоритма можно выделить следующие макрооперации:

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

На каждой итерации будет вычисляться умножение матрицы $A$ размера $n \times n$ на вектор $q_i$. Эта операция и будет самой ресурсозатратной на каждой итерации. Результаты этой макрооперации будут использованы затем в скалярном произведении векторов и в линейной комбинации из 3 векторов. После этого результат линейной комбинации векторов необходимо нормировать. Для этого можно разбить эту процедуру на вычисление нормы вектора и деление вектора на получившееся число. Кроме того, после завершения итераций необходимо вычислить собственные вектора и собственные значения полученной матрицы. Данная операция не конкретизируется алгоритмом, и может быть выполнена любым способом. Поэтому удобнее её представить в виде макрооперации. В разделе 1.7 можно наглядно увидеть как передаются данные между этими макрооперациями.

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

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

$1.\, \beta_0 = 0,\; q_0 = 0$ #Инициализируются константы

$2.\, \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}$ #Вычисляем норму вектора начального приближения.

$3.\, q_{1_{j}} = \frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots\,, n$ #Нормализуем вектор начального приближения.

Начинаются итерации цикла, которых не больше чем $k$. На $i$-й итерации производятся следующие вычисления:

$i.1\, z_j = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\,\dots\,, n$ #Считаем результат применения линейного оператора $A$ к вектору $q_i$.

$i.2\, \alpha_i = \sum\limits_{j=1}^{n}q_{i_j} z_j$ #Получаем результат скалярного произведения векторов $q_i$ и $z$.

$i.3\, z_j = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\,\dots\,, n$ #Вычисляем линейную комбинацию векторов.

$i.4\, \beta_i = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_j^2}$ #Считаем норму вектора $z$.

$i.5$ Проверка равенства $\beta_i == 0$ # Если норма оказалась равной нулю, то завершаем итерации и переходим к вычислению собственных векторов и собственных значений полученной матрицы. В обратном случае, продолжаем выполнения итераций.

$i.6\, q_{i+1_j} = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n$ #Нормируем вектор $z$.

$i.7\,$ Если выполнили $k$ итераций, то завершаем выполнение итераций, переходим к следующему шагу. Иначе начинаем последующую итерацию цикла.

$4.$ Вычисляем собственные значения и собственные вектора полученной матрицы $T_k$.

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

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

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

• вычисление нормы:
  • $n$ операций умножения
  • $n-1$ операций сложений
  • извлечение квадратного корня
• $n$ операций делений (нормирование)

Каждая из $k$ итераций состоит из:

• умножение матрицы на вектор:
  • $n^2$ умножений
  • $n(n-1)$ сложений
• скалярное произведение двух векторов
  • $n$ умножений
  • $n-1$ сложение
• вычисление линейной комбинации векторов:
  • $2n$ умножений
  • $2n$ вычитаний
• норма вектора:
  • $n$ операций умножения
  • $n-1$ операций сложений
  • извлечение квадратного корня
• нормирование вектора:
  • $n$ операций деления

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

• $k(n^2 + 4n) + n$ умножений,
• $k(n^2+3n-2) + n-1$ сложений/вычитаний,
• $kn + n$ делений,
• $k + 1$ вычислений квадратного корня.

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

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

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

Рисунок 1. Граф алгоритма без отображения входных и выходных данных. $\|x\|$ — макрооперация взятия нормы, $/$ — операция деления, $\cdot$ — операция скалярного умножения векторов, $\mathsf{F}$ — вычисление линейной комбинации векторов, $\mathsf{Ax}$ — линейный оператор, $\nabla$ — условие преждевременного выхода.

Рисунок 2. Граф алгоритма с отображением входных и выходных данных. $\|x\|$ — макрооперация взятия нормы, $/$ — операция деления, $\cdot$ — операция скалярного умножения векторов, $\mathsf{F}$ — вычисление линейной комбинации векторов, $\mathsf{Ax}$ — линейный оператор, $\nabla$ — условие преждевременного выхода, $\mathsf{In}$ — входные данные, $\mathsf{Out}$ — выходные результат.

Рисунок 3. Граф макрооперации "умножение матрицы на вектор" без отображения входных и выходных данных.
$*$ — операция умножения двух коэффициентов,
$+$ — операция сложения.

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

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

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

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

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

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

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

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

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

Входные данные: симметричная матрица $A$ (элементы $a_{ij}, i \geq j$), число $k$.

Объём входных данных: $\frac{n (n + 1)}{2} + 1$ (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы, единица относится к параметру $k$).

Выходные данные: по $k$ приближений собственных значений и собственных векторов.

Объём выходных данных: $k+kn$.

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является субквадратичным: $\frac{kn^2}{k \log{n}}$.

Вычислительная мощность алгоритма нахождения коэффициентов $\alpha_i,\beta_i$ без экономии памяти состовляет $\frac{k(2n^2+8n-1)+3n}{n^2+2k}$, при $k \ll n$, вычислительная мощность приблизительно равна $2k$.

В случае хранения половины исходной матрицы получаем соотношношение $\frac{k(2n^2+8n-1)+3n}{n(n+1)/2+2k}$, что приблизительно равно $4k$.

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

Важное свойство метода Ланцоша состоит в том, что первыми в матрице $T_{j}$ появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы $A$, находящихся на краях её спектра.

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

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

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

2.4.1 Реализация на языке C

Исследованная параллельная реализация на языке C

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

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

• gcc-5.2.0
• openmpi-1.8.4
• аргументы компилятора: -std=c11 -Ofast

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

• число процессоров [32 : 256] с шагом 16;
• размер матрицы [5000 : 200000] с шагом 5000.

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

• минимальная эффективность реализации 1.92%;
• максимальная эффективность реализации 59.47%.

График полученного распределения эффективности:

Рисунок 4. Параллельная реализация алгоритма Ланцоша. Распределение производительности.

На следующих рисунках приведены графики производительности и эффективности данной реализации в зависимости от изменяемых параметров запуска.

Рисунок 5. Параллельная реализация алгоритма Ланцоша. Изменение производительности в зависимости от числа процессоров и размера матрицы.

Рисунок 6. Параллельная реализация алгоритма Ланцоша. Изменение эффективности в зависимости от числа процессоров и размера матрицы.

Средняя эффективность для матриц с размером больше 25000 составила 52.7%.

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

Реализация в проекте IETL

Реализация в проекте ARPACK