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

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

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

Автор описания: Федотов С.А.

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

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

Алгоритм Ланцоша - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный Корнелием Ланцошом, адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности n с ограниченным числом итераций, m, где m много меньше n. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций m. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - Алгоритм Ланцоша с выборочной ортогонализацией. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

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

Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица A=A^Tи вектор начального приближения b. Мы пытаемся найти трехдиагональную симметричную матрицу T_k=Q_k^TAQ_k, собственные значения которой приближают собственные значения матрицы A. Иными словами на k-м шаге из ортонормированных векторов Ланцоша строится матрица Q_k = [q₁,q₂,...,q_k,] и в качестве приближенных собственных значений матрицы A принимаются числа Ритца. Пусть T_k=V[math]\Lambda[/math]V^T есть спектральное разложение матрицы T_k, столбцы матрицы Q_kV рассматриваются как приближения к соответсвующим собственным векторам матрицы A: [math] \begin{align} q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\ for \; & j = 1 \; to \; k \\ & z = Aq_j\\ & \alpha_j = q^T_j z\\ & z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\ & z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\ & \beta_j = \|z\|_2\\ & if \; \beta_j =0\; break\\ & q_{j+1} = \frac{z}{\beta_j}\\ end \; & for \end{align} [/math]

Диагональные элементы обозначены как [math]\alpha_j=t_{jj}[/math], а элементы побочной диагонали [math]\beta_j=t_{j-1,j}=t_{j,j-1}[/math]. После каждой итерации мы вычисляем [math]\alpha_j, \beta_j[/math], из которых строится матрица [math]T_j = \begin{pmatrix} \alpha_1 & \beta_2 & & & & 0 \\ \beta_2 & \alpha_2 & \beta_3 & & & \\ & \beta_3 & \alpha_3 & \ddots & & \\ & & \ddots & \ddots & \beta_{j-1} & \\ & & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\ 0 & & & & \beta_j & \alpha_j \\ \end{pmatrix}[/math]

Полная переортогонализация соответствует повторному проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что [math]z[/math] будет ортогонален векторам q₁,q₂,...,q_j-1. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем повторные копии сошедшихся чисел Ритца. То есть для больших k матрица T_k может иметь не одно собственное значение, очень близкое к [math]\lambda_i(A)[/math], но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи MRRR алгоритма они могут быть вычислены за O(k²) флопов.

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

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

• умножение исходной матрицы на вектор [math]z=Aq_j[/math]
• процесс ортогонализации Грама-Шмидта [math]z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i[/math]

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

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

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

В разделе #Математическое описание алгоритма уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица [math]A[/math]

Шаг 1: Вычислить [math]q_1, \beta_0, q_0[/math]

Шаг 2:В цикле последовательно вычислять матрицу [math]T_j,\alpha_j, \beta_j[/math]

• Найти собственные значения матрицы [math]T_j[/math]
• Если [math]\beta_j=0[/math] Выход

Рисунок 1. Блок-схема алгоритма

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

Для вычисления k-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:

• n+2k умножений вектора на вектор
• 2k умножение вектора на константу
• 2k сложения векторов
• O(k²) операций на разложение матрицы T_k

Легко проверить, что m шагов этого алгоритма потребуют O(m²n) флопов и O(mn) слов памяти.

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

Рисунок 2. Вычислительная схема j-й итерации. Синим цветом помечены промежуточные данные

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

Алгоритм Ланцоша с полной переортогонализацией является итерационным. На j-м шаге вычисляются числа [math]\alpha_j, \beta_j[/math] и ортонормированный вектор [math]q_{j+1}[/math], которые служат входными данными для j+1-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации j-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции. Умножение вещественной симметричной матрицы A размера [math]n \times n[/math] на вектор b длины n требует последовательного выполнения n ярусов умножений и сложений.

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

Входные данные: матрица [math]A\in\mathbb{R}^{n\times n}[/math]. Дополнительные ограничения:

• [math]A[/math] – симметричная матрица, т. е. [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math].

Объём входных данных: [math]\frac{n (n + 1)}{2}[/math].

Выходные данные:

• вектор приближенных собственных значений матрицы [math]\Lambda[/math] (элементы [math]\lambda_{ii}[/math]).
• приближения к собственным векторам [math]v_i,i=1,\ldots, m [/math].

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

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

• Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом k итераций оценивается в

[math]\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}[/math]

• Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в [math]\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k[/math] так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо [math]k\lt \lt n[/math]
• Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
• Часто уже при небольшом [math]( j \sim 2 \sqrt{n})[/math] числе шагов метода крайние пары Ритца [math](\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)[/math] оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
• Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из l, l > 1, векторов. Аналогом векторов [math]z,q_i[/math] простого алгоритма в блочном являются прямоугольные матрицы размера n * l, аналогом скаляров [math]\beta_i, \alpha_i[/math] - квадратные матрицы порядка l, аналогом трехдиагональной матрицы T_i - ленточная матрица с полушириной ленты l. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до l. Поэтому блочный алгоритм с l начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность l. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и Алгоритм Ланцоша для точной арифметики может быть дополнен выборочной ортогонализацией.

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

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

Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки TRLan на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[1] Набор и границы значений изменяемых параметров реализации алгоритма:

• размер матрицы от 5000 до 25000 с шагом 500.
• число процессоров 1, 2, 4, [8:144] с шагом 8.

Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:

Рисунок 3. Производительность в зависимости от числа процессоров и размера матрицы.

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

Можно заметить, что данная реализация немасштабируема и не дает ускорения при увеличении числа процессов. Из рисунков 3, 4 хорошо видно, что наилучшая производительность достигается с использованием только одного-двух процессов. Дальнейшее увеличение числа процессов лишь приводит к снижению эффективности. После 8-ми процессов ситуация стабилизируется и дальнейшее распараллеливание мало сказывается на результат

Использованные модули и компиляторы

Модули: mkl/4.0.2.146 impi/5.0.1-ofa компилятор mpif90, параметры компиляции: mpif90 -I/opt/intel/composer_xe_2013.2.146/mkl/include/ia32 -I/opt/intel/composer_xe_2013.2.146/mkl/include -O0 /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_lapack95_lp64.a -L/opt/intel/composer_xe_2013.2.146/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl simple.f90 -o simple lapack_aux.o ../../libtrlan_mpi.a

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

Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.

• библиотека ARPACK на языке FORTRAN77
• NAG Library
• библиотека TRLan
• GraphLab параллельная реализация на C++

3 Литература

<references \>