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

Файл:Symbol wait.svg

Эта работа прошла предварительную проверку Дата последней правки страницы: 04.10.2017 Данная работа соответствует формальным критериям. Проверено ASA.

Авторы: А.В.Ерошкин (ссылкаКод), М.М.Маликов (ссылка)

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

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

Алгоритм Ланцоша - итерационный метод, используемый для вычисления части собственных значений и соответствующих им собственных векторов матрицы [math]\displaystyle{ A }[/math] размера [math]\displaystyle{ n*n }[/math], изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра [math]\displaystyle{ k \lt \lt n }[/math], контролирующего количество итераций. Несмотря на то, что алгоритм является вычислительно эффективным, первоначально сформулированный метод был плохо применим из-за численной неустойчивости - метод хорошо работал на целочисленных значениях, однако в арифметике с плавающей точкой ошибки округления давали большую погрешность. В 1970 году Ojalvo и Newman показали, как сделать метод численно стабильным и применили его для расчета крупных инженерных сооружений, подверженных динамическим нагрузкам. Кроме того, они показали способ выбора начального приближения (с использованием ГПСЧ), а также эмпирический способ для выбора числа [math]\displaystyle{ k }[/math] (примерно в полтора раза больше искомого числа собственных векторов). В данный момент существует две основных модификации метода (с полной и выборочной переортогонализацией), а также большое количество модификаций, использующихся в различных технических областях. Алгоритм используется для больших [math]\displaystyle{ n }[/math].

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

Первый этап алгоритма - использование метода Ланцоша для построения крыловского подпространства: [math]\displaystyle{ K_k(A,x) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1] }[/math]. Входные данные алгоритма: квадратная симметричная матрица [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n*n }[/math], вектор начального приближения [math]\displaystyle{ b }[/math], а так же число итераций [math]\displaystyle{ k }[/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{ \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 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]\displaystyle{ T_k=Q_k^TAQ_k }[/math]. Ее собственные значения приближают собственные значения исходной матрицы. Пусть T_k=V[math]\displaystyle{ \Lambda }[/math]V^T - спектральное разложение матрицы T_k, тогда столбцы матрицы Q_kV рассматриваются как приближения к соответствующим собственным векторам матрицы A и называются векторами Ритца. Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы A.

Поиск собственных значений матрицы T намного легче, чем для исходной матрицы, так как предполагается, что [math]\displaystyle{ k \lt \lt n }[/math], и матрица T - трехдиагональная.

Полная переортогонализация необходима для того, чтобы гарантировать, что каждый полученный вектор q_j+1 ортогонален уже имеющимся векторам q_1..j. Без этого процесса будут накапливаться существенные вычислительный ошибки.

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

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

• [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] - умножение симметричной матрицы [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n*n }[/math] на вектор q размерности n.

• [math]\displaystyle{ z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i. }[/math] - процесс ортогонализации Грама-Шмидта.

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

Основные операции алгоритма:

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

2. Двойная ортогонализация методом Грмма-Шмидта. [math]\displaystyle{ z=z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i, }[/math] [math]\displaystyle{ z=z-\alpha_jq_j-\beta_{j}q_{j-1}. }[/math]

3. Вычисление обновленного базисного вектора. [math]\displaystyle{ q_{i+1} = z/ \beta }[/math].

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

В параграфе 1.2 приводится полная схема последовательного алгоритма.

Заполняем начальные значения алгоритма (b - начальное преближение).

[math]\displaystyle{ \begin{align} & q_1=b/||b||,\\ & \beta_1=0,\\ &q_0=0, \\ \end{align} }[/math]

Для всех [math]\displaystyle{ j=1..k }[/math]:

• 1. Вычисляется j-й диагональный элемент матрицы [math]\displaystyle{ T_k }[/math]: [math]\displaystyle{ z=Aq_j; \alpha_j=q_j^Tz; }[/math]

• 2. Проводится полная переортогонализация Грамма-Шмидта: [math]\displaystyle{ z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i; }[/math]

• 3. Вычисляются значения [math]\displaystyle{ \beta_{j+1} }[/math] матрицы [math]\displaystyle{ T_k }[/math]: [math]\displaystyle{ \beta_{j+1}=||z||; }[/math]

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

• 5. Сохраняем значения для следующей итерации [math]\displaystyle{ q_{j+1}=z/\beta_{j+1}. }[/math]

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

• 1. Основная часть операций в алгоритме Ланцоша производится во время умножения матрицы [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n*n }[/math] на вектор [math]\displaystyle{ q }[/math] размерности [math]\displaystyle{ n }[/math] - вычислительная сложность: [math]\displaystyle{ n^2 }[/math] умножений и [math]\displaystyle{ n^2-n }[/math] сложений. Остальные операции основного цикла производят меньше [math]\displaystyle{ n^2 }[/math] операций сложения или умножения. Так как умножение матрицы на вектор производится [math]\displaystyle{ k }[/math] раз, то сложность этой части алгоритма - [math]\displaystyle{ O(kn^2) }[/math]
• 2. Процесс ортогонализации Грама-Шмидта - вычислительная сложность: [math]\displaystyle{ k^2n+k(n+2) }[/math] умножений и [math]\displaystyle{ k^2n + k(n + 1) + 2 }[/math] сложений. Производится в цикле [math]\displaystyle{ k }[/math] раз. Сложность - [math]\displaystyle{ O(nk^2) }[/math]
• 3. Процесс разложения матрицы [math]\displaystyle{ T }[/math] размерности [math]\displaystyle{ k*k }[/math]. Сложность - [math]\displaystyle{ O(k^2) }[/math]

Так как число итераций много меньше размерности матрицы [math]\displaystyle{ A }[/math], [math]\displaystyle{ k \lt \lt n }[/math], то общая сложность алгоритма сокращается до [math]\displaystyle{ O(kn^2) }[/math].

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

Рис. 1. Граф j-й итерации.

Черные маленькие стрелочки - куда идут данные.

Большие стрелочки - обозначение прохода по массиву (изменение индекса)

В больших кружках преобразования (операции над данными или вычисления)

В маленьких кружках запись в переменную.

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

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

• 1. умножение матрицы размерности [math]\displaystyle{ n * n }[/math] на вектор длины [math]\displaystyle{ n }[/math] требует последовательного выполнения [math]\displaystyle{ n }[/math] ярусов умножений и сложений;
• 2. вычисление [math]\displaystyle{ \alpha_j }[/math] требует [math]\displaystyle{ n }[/math] ярусов сложений с 1 операцией умножения в каждом;
• 3. переортагонализация требует вычисления [math]\displaystyle{ j }[/math] ярусов сложений с [math]\displaystyle{ n }[/math] операциями умножения в каждом, [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ j }[/math] операциями умножения в каждом;
• 4. вычисление нормы вектора длины [math]\displaystyle{ n }[/math] требует [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ 1 }[/math] операцией умножения в каждом.

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

Входные данные: [math]\displaystyle{ A\in\mathbb{R}^{n\times n} }[/math] - симметричная матрица, т. е. [math]\displaystyle{ a_{ij}= a_{ji}, i, j = 1, \ldots, n }[/math]. Объем данных: [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math].

Выходные данные: [math]\displaystyle{ \Lambda }[/math] - вектор собственных значений матрицы и соответствующие им собственные вектора [math]\displaystyle{ v_i,i=1,\ldots, k }[/math]. Объем данных: [math]\displaystyle{ k(n+1) }[/math].

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

• В классическом алгоритме Ланцоша возникает большая погрешность при округлении чисел с плавающей точкой. Выбранный вариант с полной переортогонализацией устраняет этот недостаток, однако является более ресурсоемким. На практике наиболее популярен вариант с частичной переортогонализацией.
• Вычислительная мощность алгоритма (отношение числа операций к суммарному объему данных) оценивается как [math]\displaystyle{ \approx 2k }[/math].
• Преимуществом алгоритма является то, что он начинает поиск собственных значений матрицы начиная с максимального в абсолютном смысле значения.
• Нет необходимости хранения исходной матрицы на каждом вычислительном ядре, так как метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать регулярность структуры матрицы.
• Существует модифицированный блочный вариант алгоритма, применяемый в случае кратных собственных значений.

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

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

Проверка масштабируемости алгоритма проходила на IBM Blue Gene/P ВМК МГУ.

Используемые компиляторы: intel/15.0.090 и OpenMPI/1.8.4-icc.

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

• размеры входной матрицы: [5000:30000] c шагом 5000;
• число процессоров: от 1, 2, 4 .. 64.

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

• openmpi/1.5.5-icc
• intel/13.1.0

Рис. 2. Зависимость времени выполнения программы от входных данных

По графику видно, что время выполнения уменьшается с увеличением количества процессоров примерно до 16, дальше начинает медленно увеличиваться вплоть до 64 и далее. Это можно объяснить тем, что с увеличением количества вычислительных ядер растет количество передаваемых данных, и дальнейшее увеличение их числа только ухудшает положение.

код программы

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

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 (уже распараллелена)

Julia Math https://github.com/JuliaMath/IterativeSolvers.jl Julia

SciPy https://scipy.org/ Python

Литература

Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)

<references \>