Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)
| Файл:Symbol wait.svg | Эта работа прошла предварительную проверку Дата последней правки страницы: 03.02.2017 Данная работа соответствует формальным критериям. Проверено ASA. |
| Алгоритм Ланцоша для точной арифметики (без ортогонализации) | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(kn^2) }[/math] |
| Объём входных данных | [math]\displaystyle{ \frac{n(n+1)}{2} + 1 }[/math] |
| Объём выходных данных | [math]\displaystyle{ kn+k }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(k*log(n)) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(n^2) }[/math] |
Авторы: А.В. Шостик.
Свойства и структура алгоритма
Общее описание алгоритма
Алгоритм Ланцоша - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача [math]\displaystyle{ Ax = \lambda x }[/math]). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.
Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на понятии построения крыловского подпространства и процедуре Рэлея-Ритца: строит последовательность трехдиагональных матриц с тем свойством, что экстремальные собственные значения для каждой последующей трехдиагональной матрицы дают все более точные оценки экстремальных собственных значений для А.
В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)[1].
В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.
Математическое описание алгоритма
Исходные данные:
- симметрическая матрица [math]\displaystyle{ A=A^T }[/math] размерности [math]\displaystyle{ n }[/math], элементы матрицы [math]\displaystyle{ a_{ij}=a_{ji} }[/math]
- начальный вектор [math]\displaystyle{ v = {v_1, v_2, ..., v_n}, v \neq \theta }[/math]
Так как матрица является симметрической, достаточно хранить [math]\displaystyle{ \frac{n(n+1)}{2} }[/math] ее элементов.
Для исходной матрицы [math]\displaystyle{ A }[/math] и вектора [math]\displaystyle{ v }[/math] строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов [math]\displaystyle{ [v, Av, A^2v, ..., A^{k-1}v] }[/math], называемых векторами Ланцоша. Каждый из векторов Ланцоша ортонормируется: [math]\displaystyle{ q_k=\frac{A^kv}{\|A^kv\|} }[/math], и на шаге [math]\displaystyle{ k }[/math] алгоритма имеется Крыловская матрица [math]\displaystyle{ Q = {q_0, ..., q_{k-1}} }[/math] размерности [math]\displaystyle{ n \times k }[/math].
Соответствующий текущей итерации [math]\displaystyle{ k }[/math] вектор Ланцоша считается [math]\displaystyle{ k }[/math]-ым приближением собственного вектора исходной матрицы [math]\displaystyle{ A }[/math]. В качестве приближенных собственных значений матрицы [math]\displaystyle{ A }[/math] берутся собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_k = Q^T_k A Q }[/math] (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует [math]\displaystyle{ O(k^3) }[/math] операций.[2].
Выходные данные:
- совокупность собственных векторов матрицы [math]\displaystyle{ T_k }[/math] (т.е. [math]\displaystyle{ k }[/math]-ых приближений собственных векторов матрицы [math]\displaystyle{ A }[/math])
- совокупность собственных значений матрицы [math]\displaystyle{ T_k }[/math] (т.е. [math]\displaystyle{ k }[/math]-ых приближений собственных значений матрицы [math]\displaystyle{ A }[/math])
Вычислительное ядро алгоритма
Вычислительное ядро. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора [math]\displaystyle{ z=Aq_i }[/math]
Макроструктура алгоритма
Итерация алгоритма:
- вычисление нормы вектора,
- деление вектора на число,
- умножение матрицы на вектор,
- линейная комбинация векторов (умножение на число и сложение).
Финальный расчет:
- вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.
Схема реализации последовательного алгоритма
Псевдокод алгоритма[3]:
Вход : размерность матрицы [math]\displaystyle{ n }[/math], элементы симметричной матрицы [math]\displaystyle{ A }[/math], количество итераций [math]\displaystyle{ k }[/math] Выход : собственные вектора матрицы [math]\displaystyle{ T_k }[/math], матрица собственных значений
[math]\displaystyle{ \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]\displaystyle{ procedure(T_k) }[/math];
[math]\displaystyle{ procedure() }[/math] - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.
Последовательная сложность алгоритма
Глобально алгоритм состоит из инициализации, [math]\displaystyle{ k }[/math] итераций и финального расчета собственных векторов и собственных значений.
Рассмотрим сложность каждой макрооперации алгоритма:
- вычисление нормы вектора - [math]\displaystyle{ (n-1) }[/math] сложений, [math]\displaystyle{ n }[/math] умножений, операция извлечения корня
- деление вектора на число - [math]\displaystyle{ n }[/math] операций
- умножение матрицы на вектор - [math]\displaystyle{ n^2 }[/math] операций
- скалярное произведение векторов - [math]\displaystyle{ n }[/math] операций
- линейная комбинация векторов - [math]\displaystyle{ 2n }[/math] умножений и вычитаний
- вычисление собственных векторов и собственных значений матрицы [math]\displaystyle{ T_k }[/math] методом "Разделяй и властвуй" - [math]\displaystyle{ O(k^3) }[/math] операций.
Инициализация:
- вычисление нормы вектора,
- деление вектора на число.
Одна итерация:
- умножение матрицы на вектор,
- скалярное произведение векторов,
- линейная комбинация векторов,
- вычисление нормы вектора,
- деление вектора на число.
Финальный расчет:
- вычисление собственных векторов и собственных значений матрицы
Итого последовательная сложность алгоритма составляет [math]\displaystyle{ O(3n+k(n^2+n+2n+2n+3n))+O(k^3) }[/math].
Так как на практике порядок матрицы [math]\displaystyle{ n }[/math] намного превышает количество итераций [math]\displaystyle{ k }[/math], сложность алгоритма Ланцоша можно рассматривать [math]\displaystyle{ O(kn^2) }[/math].
Информационный граф
Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:
[math]\displaystyle{ \|.\| }[/math] — операция вычисления нормы,
[math]\displaystyle{ / }[/math] — операция деления,
[math]\displaystyle{ A }[/math] — операция умножения матрицы на вектор,
[math]\displaystyle{ * }[/math] — операция скалярного произведения,
[math]\displaystyle{ L }[/math] — операция вычисления линейной комбинации векторов,
▽ - проверка условия [math]\displaystyle{ \mathsf{(\beta_i = 0)} }[/math] и выход из цикла в случае выполнения условия
Ресурс параллелизма алгоритма
Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации [math]\displaystyle{ n }[/math] умножений и [math]\displaystyle{ n-1 }[/math] сложение. В таком виде сложение [math]\displaystyle{ n }[/math] элементов можно выполнить за [math]\displaystyle{ \log_2 n }[/math] действий.
Входные и выходные данные алгоритма
Входные данные: симметрическая матрица [math]\displaystyle{ A=A^T }[/math] размерности [math]\displaystyle{ n }[/math], количество итераций [math]\displaystyle{ k }[/math].
Объем входных данных: [math]\displaystyle{ \frac{n(n+1)}{2} + 1 }[/math].
Выходные данные: [math]\displaystyle{ k }[/math] собственных векторов, собственных значений матрицы на [math]\displaystyle{ k }[/math]-ом приближении.
Объем выходных данных: [math]\displaystyle{ kn+k }[/math].
Свойства алгоритма
- Сложность последовательной реализации алгоритма [math]\displaystyle{ O(kn^2) }[/math].
- Сложность параллельного алгоритма Ланцоша по высоте ЯПФ [math]\displaystyle{ O(k\log n) }[/math], по ширине ЯПФ [math]\displaystyle{ O(n^2) }[/math].
- В умножении матрицы на вектор сложение n элементов выполняется в [math]\displaystyle{ \log_2 n }[/math] ярусов шириной [math]\displaystyle{ \frac{n}{2}, \frac{n}{4}, ... , 1 }[/math], остальные операции внутри итерации выполняются последовательно
- Отношение последовательной сложности к параллельной [math]\displaystyle{ \frac{kn^2}{k \log{n}} }[/math].
- Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма [math]\displaystyle{ \frac{k(2n^2+8n-1)+3n}{n^2+2k} }[/math]. Учитывая [math]\displaystyle{ k }[/math] много меньше [math]\displaystyle{ n }[/math], вычислительная мощность ≈ [math]\displaystyle{ 2k }[/math].
- В алгоритме Ланцоша возможно выполнение числа итераций менее [math]\displaystyle{ k }[/math] если все собственные значения матрицы вычисляются раньше.
- Для найденных собственных значений справедливо: [math]\displaystyle{ \lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k}) }[/math], т.е. в первую очередь находятся максимальные по модулю собственные значения.
- Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Исследование проводилось на суперкомпьютере "Ломоносов"[1].
Исследуемая реализация на языке C [2] [4].
Компилятор gcc, аргумент -std=c11. Openmpi-1.8.4
- Число процессоров варьировалось от 32 до 256, step = 16.
- Размер матрицы изменялся от [math]\displaystyle{ 5 000 \times 5 000 }[/math] до [math]\displaystyle{ 200 000 \times 200 000 }[/math], step = 5 000.
Ниже приведена зависимость эффективности параллельной реализации от указанных выше параметров.
[math]\displaystyle{ Size }[/math] — размерность матрицы [math]\displaystyle{ N }[/math],
[math]\displaystyle{ Cores }[/math] — количество процессоров,
[math]\displaystyle{ Eff }[/math] — эффективность алгоритма в %.
Выводы:
- минимальная эффективность алгоритма при указанных параметрах составила 1.9% (при числе процессоров 32, размерности матрицы 5 000),
- максимальная эффективность алгоритма при указанных параметрах составила 59.5% (при числе процессоров 256, размерности матрицы 200 000),
- средняя эффективность алгоритма при числе процессоров от 32 до 256 и размерности матрицы от 25 000 до 200 000 составила 52,7%,
- при размерности матрицы менее 25 000 эффективность алгоритма менее 23%.
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе [5][6][7] проекты на базе данного алгоритма.
- UNDERWOOD (1977) [8]. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
- PROPACK (1998) [9]. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
- BLKLAN (2004) [10]. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
- IETL (2006) [11]. Только без реортогонализации. Язык C++.
- SLEPc (2009) [12]. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
- Lanczos Plus Plus (2009) [13]. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.
Литература
- ↑ Деммель Д. Вычислительная линейная алгебра
- ↑ https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm
- ↑ Деммель Д. Вычислительная линейная алгебра
- ↑ https://github.com/xrstalker/lanczos
- ↑ Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.
- ↑ Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.
- ↑ http://www.grycap.upv.es/slepc
- ↑ Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA
- ↑ http://www.daimi.au.dk/PB/537
- ↑ Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada
- ↑ http://www.comp-phys.org/software/ietl/
- ↑ Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.
- ↑ https://github.com/g1257/LanczosPlusPlus