Участник:Danyanya/Алгоритм Ланцоша для точной арифметики (без переортогонализации)
| Алгоритм Ланцоша для точной арифметики (без переортогонализации) | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(kn^2) }[/math] |
| Объём входных данных | [math]\displaystyle{ \frac{n(n + 1)}{2} }[/math] |
| Объём выходных данных | [math]\displaystyle{ k(n + 1) }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(k) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(n^2) }[/math] |
Основные авторы описания: Д.Р.Слюсарь (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), М.А.Григорьев (1.2, 1.7, 1.8, 1.9, 2.7, 3)
Свойства и структура алгоритма
Общее описание алгоритма
Алгоритм Ланцоша поиска собственных значений был опубликован Корнелием Ланцошем в 1950 году [1]. Этот итерационный алгоритм применим только к эрмитовым матрицам [math]\displaystyle{ A }[/math]. Метод позволяет за [math]\displaystyle{ k }[/math] итераций вычислять [math]\displaystyle{ k }[/math]-ое приближение собственных значений и собственных векторов исходной матрицы [math]\displaystyle{ A }[/math].
В данной статье рассмотрен упрощенный вариант алгоритма Ланцоша, подразумевающий отсутствие влияния ошибок округления на вычислительный процесс.
Данный алгоритм является неустойчивым, вследствие чего на практике применяется модифицированный алгоритм Ланцоша с полной переортогонализацией предложенный в 1970 Ojalvo и Newman[2].
Математическое описание алгоритма
На вход алгоритма подается эрмитова матрица [math]\displaystyle{ A = A^\dagger }[/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]
Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[1]. Иными словами, из оргонормированных векторов Ланцоша [3] на каждой итерации строится матрица [math]\displaystyle{ Q_k = [q_1, q_2, \dots, q_k] }[/math] размерности [math]\displaystyle{ n \times k }[/math]. В качестве приближенных собственных значений матрицы [math]\displaystyle{ A }[/math] берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_k = Q^T_k A Q }[/math]:
[math]\displaystyle{ T_k = \begin{pmatrix} \alpha_1 & \beta_1 & 0 & \dots & 0 \\ \beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\ 0 & \beta_2 & \ddots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\ 0 & \dots & \dots & \beta_{k-1} & \alpha_k \end{pmatrix} }[/math]
На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы [math]\displaystyle{ T_k }[/math], с помощью которых и будет найдены искомые собственные векторы исходной матрицы [math]\displaystyle{ A }[/math].
Вычислительное ядро алгоритма
Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы [math]\displaystyle{ A }[/math] на вектор [math]\displaystyle{ q_i }[/math] с предыдущей итерации
- [math]\displaystyle{ z = Aq_i }[/math]
Макроструктура алгоритма
Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:
- Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
- умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
- скалярное произведение векторов;
- линейная комбинация векторов (сложение/умножение на вещественные числа);
- вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
- Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.
При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.
Вычисление собственных значений полученной матрицы на практике выполняется с помощью QR-алгоритма и выходит за рамки описанного алгоритма.
Схема реализации последовательного алгоритма
Исходные данные: симметричная матрица [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].
Алгоритм [4] на псевдокоде:
[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{ T_k }[/math] наиболее удобным образом.
Последовательная сложность алгоритма
Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
- Умножение квадратной матрицы [math]\displaystyle{ n * n }[/math] на вектор длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n * n }[/math] умножений и сложений;
- Скалярное произведение векторов длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] умножений и сложений;
- Сложение векторов длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] сложений;
- Умножение вектора длины [math]\displaystyle{ n }[/math] на скаляр. Требует [math]\displaystyle{ n }[/math] умножений;
- Нахождение квадратичной нормы вектора длины [math]\displaystyle{ n }[/math]. Требует [math]\displaystyle{ n }[/math] умножений и сложений + извлечения квадратного корня;
- Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math]\displaystyle{ k \times k }[/math]. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует [math]\displaystyle{ ~O(k^{2.3}) }[/math] операций.
Таким образом, в худшем случае, алгоритм Ланцоша имеет сложность [math]\displaystyle{ O(k * n^2) }[/math] и односится к квадратичным алгоритмам.
Информационный граф
Информационный граф алгоритма можно разбить на две части:
- Общий граф алгоритма, от входа IN до функции вычислений собственных значений полученной трехдиагональной матрицы (EIGEN) (Рис. 1);
- Подграф на каждой итерации внутри цикла (от IN до NEXT в рисунке 2);
Ресурс параллелизма алгоритма
Важно заметить, что итерации алгоритма выполняются строго последовательно, а распараллеливание возможно только внутри итераций.
- Умножение [math]\displaystyle{ A^{(n * n)} }[/math] на вектор длины [math]\displaystyle{ n }[/math] требует [math]\displaystyle{ n }[/math] ярусов умножений и сложений.
- Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
- Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.
Исходя из вышеизложенного, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и по высоте ЯПФ.
Входные и выходные данные алгоритма
Входные данные: симметричная вещественная матрица [math]\displaystyle{ A }[/math], случайный вектор [math]\displaystyle{ b }[/math], число итераций [math]\displaystyle{ k }[/math].
Объём входных данных: [math]\displaystyle{ n * (n + 1) + 1 }[/math].
Выходные данные: вектор собственных значений [math]\displaystyle{ \Lambda }[/math], матрица собственных векторов [math]\displaystyle{ E }[/math].
Объём выходных данных: [math]\displaystyle{ k * (n + 1) }[/math].
Свойства алгоритма
- Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
- Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы [math]\displaystyle{ A }[/math], находящихся на границе ее спектра (в [math]\displaystyle{ T_{j} }[/math] в первую очередь появляются максимальные по модулю собственные значения);
- Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
- Нестабильность алгоритма (эффект ложной сходимости) присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
TDB next week
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (ARPACK), так и неофициальных реализаций, выложенных на Github. Среди них:
1. The IETL Project [5]
2. MatLab [6]
3. ARPACK [7]
4. Julia Math [8]
Также существуют официальные реализации на других языках (например R).
Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.
Литература
- ↑ 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).
- ↑ Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).
- ↑ https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша
- ↑ Деммель Д. Вычислительная линейная алгебра
- ↑ http://www.comp-phys.org/software/ietl/
- ↑ https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554
- ↑ http://www.caam.rice.edu/software/ARPACK/
- ↑ https://github.com/JuliaMath/IterativeSolvers.jl