Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
Файл:Symbol wait.svgЭта работа ждет рассмотрения преподавателем
Дата последней правки страницы: 18.11.2016
Авторы этой статьи считают, что задание выполнено.



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

Основные авторы описания: Заспа А.Ю. (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)

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

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

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

[math]\displaystyle{ 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]\displaystyle{ Q_k = [q_1, q_2, \dots, q_k] }[/math] размерности [math]\displaystyle{ n \times k }[/math], состоящая из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_k = Q^T_k A Q }[/math] размерности [math]\displaystyle{ k \times k }[/math].

[math]\displaystyle{ 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]\displaystyle{ O(k^3) }[/math] операций, где константа оказывается на практике довольно мала.

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

Исходные данные: симметрическая матрица [math]\displaystyle{ A }[/math], начальный вектор [math]\displaystyle{ b }[/math].

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

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

На [math]\displaystyle{ i }[/math]-й итерации происходят следующие вычисления. Вначале вычисляется произведение матрицы на вектор [math]\displaystyle{ z = Aq_i }[/math]. Затем считается скалярное произведение [math]\displaystyle{ \alpha_i = q^T_i z }[/math] и значение сохраняется как элемент формируемой матрицы. После этого необходимо произвести вычисление линейной комбинации векторов [math]\displaystyle{ z = z - \alpha_i q_i - \beta_{i-1}q_{i-1} }[/math]. Норму полученного вектора запишем как элемент матрицы [math]\displaystyle{ \beta_i = \|z\|_2 }[/math]. Если [math]\displaystyle{ \beta_i }[/math] оказалась равной нулю, то больше никаких итераций в алгоритме не производят, сразу переходят к вычислению собственных значений и собственных векторов. Иначе, в конце итерации считают [math]\displaystyle{ q_{i+1} = \frac{z}{\beta_i} }[/math] и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы [math]\displaystyle{ T_k }[/math].[3]

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

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

[math]\displaystyle{ z = Aq_i }[/math]

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

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

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

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


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

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

[math]\displaystyle{ 1.\, \beta_0 = 0,\; q_0 = 0 }[/math] #Инициализируются константы

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

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

Начинаются итерации цикла, которых не больше чем [math]\displaystyle{ k }[/math]. На [math]\displaystyle{ i }[/math]-й итерации производятся следующие вычисления:

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

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

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

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

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

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

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

[math]\displaystyle{ 4. }[/math] Вычисляем собственные значения и собственные вектора полученной матрицы [math]\displaystyle{ T_k }[/math].

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

Вычисления [math]\displaystyle{ k }[/math] собственных значений матрицы порядка [math]\displaystyle{ n }[/math] и соответствующих им собственных векторов алгоритмом Ланцоша без переортогонализации состоит из инициализации, [math]\displaystyle{ k }[/math] итераций и нахождения собственных значений и векторов трехдиагональной симметричной матрицы порядка [math]\displaystyle{ k }[/math].

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

  • вычисление нормы:
    • [math]\displaystyle{ n }[/math] операций умножения
    • [math]\displaystyle{ n-1 }[/math] операций сложений
    • извлечение квадратного корня
  • [math]\displaystyle{ n }[/math] операций делений (нормирование)

Каждая из [math]\displaystyle{ k }[/math] итераций состоит из:

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

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

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

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

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

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

Рисунок 1. Граф алгоритма без отображения входных и выходных данных.
[math]\displaystyle{ \|x\| }[/math] — макрооперация взятия нормы,
[math]\displaystyle{ / }[/math] — операция деления,
[math]\displaystyle{ \cdot }[/math] — операция скалярного умножения векторов,
[math]\displaystyle{ \mathsf{F} }[/math] — вычисление линейной комбинации векторов,
[math]\displaystyle{ \mathsf{Ax} }[/math] — линейный оператор,
[math]\displaystyle{ \nabla }[/math] — условие преждевременного выхода.
Рисунок 2. Граф алгоритма с отображением входных и выходных данных.
[math]\displaystyle{ \|x\| }[/math] — макрооперация взятия нормы,
[math]\displaystyle{ / }[/math] — операция деления,
[math]\displaystyle{ \cdot }[/math] — операция скалярного умножения векторов,
[math]\displaystyle{ \mathsf{F} }[/math] — вычисление линейной комбинации векторов,
[math]\displaystyle{ \mathsf{Ax} }[/math] — линейный оператор,
[math]\displaystyle{ \nabla }[/math] — условие преждевременного выхода,
[math]\displaystyle{ \mathsf{In} }[/math] — входные данные,
[math]\displaystyle{ \mathsf{Out} }[/math] — выходные результат.
Ошибка создания миниатюры: Файл не найден
Рисунок 3. Граф макрооперации "умножение матрицы на вектор" без отображения входных и выходных данных.
[math]\displaystyle{ * }[/math] — операция умножения двух коэффициентов,
[math]\displaystyle{ + }[/math] — операция сложения.

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

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

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

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

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


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

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


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

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

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

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

Выходные данные: по [math]\displaystyle{ k }[/math] приближений собственных значений и собственных векторов.

Объём выходных данных: [math]\displaystyle{ k+kn }[/math].

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

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

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

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

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

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

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

Особенности реализации последовательного алгоритма

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

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

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

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

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

  • gcc-5.2.0
  • openmpi-1.8.4
  • mpicc main.c -o out

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

  • число процессоров [64 : 224] с шагом 16 / 32;
  • размер матрицы [25000 : 150000] с шагом 5000.

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

  • минимальная эффективность реализации 4,55%;
  • максимальная эффективность реализации 11,17%.

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

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

Как видно из графиков, производительность линейно зависит от количества процессоров, при этом минимальная эффективность наблюдается на матрицах размером 25000. Для всех последующих размеров минимальная эффективность составила 8.39% с попаданием большинства значений в интервал 10.5%-11%.

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

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

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

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

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

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

Ссылки

  1. ↑ Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).
  2. ↑ Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).
  3. ↑ Деммель Д. Вычислительная линейная алгебра
  4. ↑ http://www.comp-phys.org/software/ietl/lanczos.html
  5. ↑ http://www.caam.rice.edu/software/ARPACK/
Источник — https://xn--80aecvaldz.xn--80ahdrnmvk.xn--p1ai/w/ru/index.php?title=Участник:A.Freeman/Алгоритм_Ланцоша_для_точной_арифметики_(без_переортогонализации)&oldid=18384