Участник:A.Freeman/Алгоритм Ланцоша для точной арифметики (без переортогонализации)
| Алгоритм Ланцоша без переортогонализации | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(kn^2) }[/math] |
| Объём входных данных | [math]\displaystyle{ \frac{n (n + 1)}{2}+1 }[/math] |
| Объём выходных данных | [math]\displaystyle{ k + kn }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(k \log n) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(n^2) }[/math] |
Основные авторы описания: Заспа А.Ю., Фролов А.А.
Свойства и структура алгоритма
Общее описание алгоритма
Алгоритм Ланцоша был опубликовн физиком и математиком Корнелием Ланцощем в 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{ \begin{align} q_1 = & \frac{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 \; break\\ & q_{i+1} = \frac{z}{\beta_i}\\ end \; & for \end{align} }[/math]
Затем вычислить собственные значения и собственные вектора матрицы [math]\displaystyle{ T_k }[/math]
Вычислительное ядро алгоритма
Вычислительным ядром на каждой итерации является вычисление произведения матрицы на вектор:
- [math]\displaystyle{ z = Aq_i }[/math]
Макроструктура алгоритма
Макрооперациями алгоритма являются: умножение матрицы на вектор, скалярное произведение векторов, умножение вектора на число, вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.
Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
[math]\displaystyle{ \begin{align} 1. \, \beta_0 = & 0,\; q_0 = 0,\\ \|b\|_2 = &\sqrt{\sum\limits_{j=1}^{n} b_j^2},\\ 2. \, q_{1_{j}} = &\frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots \,, n\\ for \; i = & 1 \; to \; k \\ 3. \, z_j & = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\, \dots \,, n\\ 4. \, \alpha_i & = \sum\limits_{j=1}^{n}q_j z_j\\ 5. \, z_j & = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\, \dots \,, n\\ 6. \, \beta_i & = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_m^2}\\ if \; \beta_i & == 0 \; then \; break\\ 7. \, q_{i+1_j} & = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n\\ end \; & for \end{align} }[/math]
Последовательная сложность алгоритма
Для вычисления [math]\displaystyle{ k }[/math] собственных значений матрицы порядка [math]\displaystyle{ n }[/math] и соответствующих им собственных векторов алгоритмом Ланцоша без переортогонализации в худшем случае требуется:
- [math]\displaystyle{ kn^2 + 4 kn + n }[/math] умножений,
- [math]\displaystyle{ k + 1 }[/math] вычислений квадратного корня,
- [math]\displaystyle{ n(k + 1) }[/math] делений,
- [math]\displaystyle{ n + k(n^2+3n-1) }[/math] сложений (вычитаний),
- нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера [math]\displaystyle{ k \times k }[/math].
Умножения и сложения (вычитания) составляют основную часть алгоритма.
Для нахождение собственных значений и векторов трехдиагональной симметричной матрицы эффективней всего использовать метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы, который требует [math]\displaystyle{ O(k^3) }[/math] операций. При классификации по последовательной сложности, таким образом, метод Ланцоша без переортогонализации относится к алгоритмам с квадратичной сложностью.
Информационный граф
TBD
Ресурс параллелизма алгоритма
Алгоритм Ланцоша в параллельной форме состоит из инициализации и [math]\displaystyle{ k }[/math] итераций. Инициализации состоит из следующих ярусов:
- ярус с [math]\displaystyle{ n }[/math] умножений
- ярус с [math]\displaystyle{ n-1 }[/math] сложений
- ярус с вычислением квадратного корня
- ярус [math]\displaystyle{ n }[/math] делений
Заметим, что вычисление суммы [math]\displaystyle{ n }[/math] элементов имеет высоту [math]\displaystyle{ \log n }[/math] и линейную ширину ярусов [math]\displaystyle{ n/2, n/4, ..., 1 }[/math].
Итерационная часть алгоритма состоит из следующих ярусов:
- 4 яруса умножений с максимальной шириной [math]\displaystyle{ n^2 }[/math]
- 3 блока вычислений сумм [math]\displaystyle{ n }[/math] элементов, состоящих из [math]\displaystyle{ \log n }[/math] ярусов и наибольшей шириной [math]\displaystyle{ \frac{n^2}{2} }[/math]
- ярус из двух последовательных вычитаний с шириной [math]\displaystyle{ n }[/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} }[/math].
Программная реализация
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации
Реализация в проекте IETL[3]
Реализация в проекте ARPACK [4]
Ссылки
- ↑ 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).
- ↑ http://www.comp-phys.org/software/ietl/lanczos.html
- ↑ http://www.caam.rice.edu/software/ARPACK/