|
|
Строка 1: |
Строка 1: |
| + | {{Assignmenta}} |
| Авторы: Хакимов А. С. | | Авторы: Хакимов А. С. |
| | | |
Строка 39: |
Строка 40: |
| <math>\alpha_j=q_j^Tz,</math> | | <math>\alpha_j=q_j^Tz,</math> |
| <math>z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},</math> | | <math>z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},</math> |
− | <math>\text{/* Провести выборочную ортогонализацию по отношению}, </math> | + | /* Провести выборочную ортогонализацию по отношению |
− | <math>\text{ к сошедшимся векторам Ритца */}</math> | + | к сошедшимся векторам Ритца */ |
− | <math>for\, i \leqslant k, \text{таких, что} \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math> | + | <math>for\, i \leqslant k, \,</math>таких, что<math \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math> |
| <math>z = z-(y^T_{i,k},z)y_{i,k}</math> | | <math>z = z-(y^T_{i,k},z)y_{i,k}</math> |
− | <math>end \, for</math> | + | <math>\text{end for}</math> |
| <math>\beta_{j}=\|z\|_2</math> | | <math>\beta_{j}=\|z\|_2</math> |
| <math>q_{j+1}=z/\beta_{j}, </math> | | <math>q_{j+1}=z/\beta_{j}, </math> |
− | <math>\text{Вычислить собственные значения и собственные векторы}</math> | + | Вычислить собственные значения и собственные векторы |
− | <math>\text{матрицы} \, \, T_{j} \, \text{и оценки погрешности в них}</math> | + | матрицы<math>\, \, T_{j} \, \,</math>и оценки погрешности в них |
− | <math>end \, for</math>
| + | <math>\text{end for}</math> |
− | | |
− | == Вычислительное ядро алгоритма ==
| |
− | | |
− | * Умножение матрицы на вектор, <math>z=Aq_j, </math>.
| |
− | | |
− | * Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>
| |
− | | |
− | == Макроструктура алгоритма ==
| |
− | | |
− | Для алгоритма Ланцоша можно выделить следующие макрооперации:
| |
− | | |
− | * Умножение матрицы на вектор. Состоит из операций умножения вектора на число и сложения векторов.
| |
− | | |
− | <math>z=Aq_j,</math>.
| |
− | | |
− | * Вычисление вектора <math>q_{j+1}</math> как скалярное произведение двух других векторов.
| |
− | | |
− | <math>\alpha_j=q_j^Tz,</math> | |
− | | |
− | <math>z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},</math>
| |
− | | |
− | <math>q_{j+1}=z/\|z\|_2</math>
| |
− | | |
− | * Ортогонализация вектора Ланцоша с помощью скалярного произведения:
| |
− | | |
− | <math>z = z-(y^T_{i,t},z)y_{i,t}</math>
| |
− | | |
− | * Алгоритм "разделяй и властвуй" для вычисления собственных значений матрицы.
| |
− | | |
− | == Схема реализации последовательного алгоритма ==
| |
− | | |
− | Алгоритм выполняется в следующей последовательности:
| |
− |
| |
− | <math>1.\, \, \beta_0 = 0,\; q_0 = 0</math> #Инициализируем векторы
| |
− |
| |
− | <math>2.\, \, \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}</math> #Вычисляем норму вектора начального приближения.
| |
− |
| |
− | <math>3.\, \, q_{1_{j}} = \frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots\,, n</math> #Нормализуем вектор начального приближения.
| |
− |
| |
− | <math>4.\, \,</math> Начинается цикл, в котором не более чем <math>k</math>. На <math>i</math>-й итерации производятся следующие вычисления:
| |
− |
| |
− | <math>4.1.\, \, z_j = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\,\dots\,, n</math> #Считаем результат применения линейного оператора <math>A</math> к вектору <math>q_i</math>.
| |
− |
| |
− | <math>4.2.\, \, \alpha_i = \sum\limits_{j=1}^{n}q_{i_j} z_j</math> #Получаем результат скалярного произведения векторов <math>q_i</math> и <math>z</math>.
| |
− |
| |
− | <math>4.3.\, \, z_j = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\,\dots\,, n</math> #Вычисляем линейную комбинацию векторов.
| |
− |
| |
− | <math>4.4.\, \,</math> Начинается цикл, в котором производится выборочная ортогонализация векторов Ланцоша
| |
− |
| |
− | <math>4.4.1.\, \,</math> Вычисляется вектор Ритца <math>y_{i,k_d} = \sum\limits_{x=1}^{k} q_{x_d} v_{i_x}, \; d = 1,\,\dots\,, n,</math>, где <math>y_{i,k_d}</math> - <math>d</math>-я компонента вектора Ритца <math>y,</math> | |
− | <math>4.4.2.\, \,</math> Производится скалярное умножение <math>y_{i,k}^T</math> и <math> z: \; \; \gamma = \sum\limits_{d=1}^{n}y_{i,k_d} z_d,</math>
| |
− | <math>4.4.3.\, \,</math> Ищется новое собственное значение <math> \theta_j </math> и собственный вектор <math> v_j </math> для полученной матрицы <math> T_j,</math>
| |
− |
| |
− | <math>4.5.\, \, \beta_i = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_j^2}</math> #Считаем норму вектора <math>z</math>.
| |
− |
| |
− | <math>4.6.</math> Проверка равенства <math>\beta_i == 0</math> # Если норма оказалась равной нулю, то завершаем итерации и переходим к вычислению собственных векторов и собственных значений полученной матрицы. В обратном случае, продолжаем выполнения итераций.
| |
− |
| |
− | <math>4.7.\, \, q_{i+1_j} = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n</math> #Нормируем вектор <math>z</math>.
| |
− |
| |
− | <math>4.8.\,</math> Если выполнили <math>k</math> итераций, то завершаем выполнение итераций, переходим к следующему шагу. Иначе начинаем последующую итерацию цикла.
| |
− |
| |
− | <math>5.</math> Вычисляем собственные значения и собственные вектора полученной матрицы <math>T_k</math>.
| |
− |
| |
− | == Последовательная сложность алгоритма ==
| |
− | | |
− | * Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
| |
− | * Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
| |
− | * Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
| |
− | * Оценка погрешностей - <math>kn</math> операций умножения
| |
− | * Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
| |
− | * Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
| |
− | * Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения
| |
− | | |
− | Итого на каждой итерации:
| |
− | | |
− | * Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
| |
− | * Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>
| |
− | | |
− | Итого за все время выполнения алгоритма:
| |
− | | |
− | * Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
| |
− | * Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>
| |
− | | |
− | == Информационный граф ==
| |
− | | |
− | [[File:anczos_graph.png|1024px|center|рис. 1: Информационный граф]]
| |
− | Здесь:
| |
− | <math>Ax</math> — перемножение матрицы и вектора,
| |
− | <math>\|x\|</math> — макрооперация взятия нормы,
| |
− | <math>\alpha x</math> — умножение скаляра на вектор,
| |
− | <math>\cdot </math> — операция скалярного умножения векторов,
| |
− | <math>\mathsf{F}</math> — операция <math>z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i</math>,
| |
− | <math>\mathsf{f}</math> — вычисление линейной комбинации векторов
| |
− | | |
− | == Ресурс параллелизма алгоритма ==
| |
− | | |
− | На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:
| |
− | | |
− | * Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
| |
− | * Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
| |
− | * Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
| |
− | * Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
| |
− | * Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
| |
− | * Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>
| |
− | | |
− | Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.
| |
− | | |
− | == Входные и выходные данные алгоритма ==
| |
− | | |
− | Входные данные: вещественная матрица <math>A</math> размера <math>n \times n</math>, вектор <math>b</math> длины <math>n</math>, число <math>k</math>
| |
− | | |
− | Объем входных данных: <math>n^2 + n + 1</math>
| |
− | | |
− | Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>
| |
− | | |
− | Объем выходных данных: <math>nk</math>
| |
− | | |
− | == Свойства алгоритма ==
| |
− | | |
− | = Программная реализация =
| |
− | | |
− | == Особенности реализации последовательного алгоритма ==
| |
− | | |
− | == Локальность данных и вычислений ==
| |
− | | |
− | == Возможные способы и особенности параллельной реализации алгоритма ==
| |
− | | |
− | == Масштабируемость алгоритма и его реализации ==
| |
− | | |
− | Исследование алгоритма на масштабируемость производилось на суперкомпьютере [http://parallel.ru/cluster "Ломоносов"] с использованием технологии MPI. Для тестирования была выбрана [http://pastebin.com/SbNrh7nz вот эта] реализация алгоритма.
| |
− | | |
− | Набор параметров, которые изменялись при тестировании:
| |
− | *Число процессоров от 16 до 128.
| |
− | *Размерность матрицы от 5000 до 50000 с шагом 5000.
| |
− | | |
− | {| class="wikitable" style="border: none; background: none;"
| |
− | |+ align="bottom" style="caption-side: bottom" | Табл. 1. Результаты тестирования. В ячейках таблицы -- время работы алгоритма в секундах.
| |
− | ! colspan="2" rowspan="2" style="border: none; background: none;"|
| |
− | ! colspan="4"| Число процессоров
| |
− | |-
| |
− | ! 16 !! 32 !! 64 !! 128
| |
− | |-
| |
− | ! rowspan="10"| Размер матрицы
| |
− | ! 5000
| |
− | | 1,07 || 0,71 || 0,62 || 0,54
| |
− | |-
| |
− | ! 10000
| |
− | | 2,39 || 1,94 || 1,17 || 0,79
| |
− | |-
| |
− | ! 15000
| |
− | | 4,35 || 3,38 || 1,89 || 1,59
| |
− | |-
| |
− | ! 20000
| |
− | | 7,63 || 5,19 || 3,78 || 2,35
| |
− | |-
| |
− | ! 15000
| |
− | | 12,41 || 7,33 || 4,32 || 3,44
| |
− | |-
| |
− | ! 30000
| |
− | | 16,68 || 9,36 || 5,02 || 4,68
| |
− | |-
| |
− | ! 35000
| |
− | | 22,58 || 11,92 || 8,04 || 6,20
| |
− | |-
| |
− | ! 40000
| |
− | | 29,31 || 15,90 || 8,47 || 7,43
| |
− | |-
| |
− | ! 45000
| |
− | | 37,04 || 19,42 || 10,55 || 9,04
| |
− | |-
| |
− | ! 50000
| |
− | | 45,64 || 23,39 || 16,09 || 10,25
| |
− | |}
| |
− | | |
− | [[File:grph.png|1024px|center|рис. 1: Информационный граф]]
| |
− | | |
− | == Динамические характеристики и эффективность реализации алгоритма ==
| |
− | | |
− | == Выводы для классов архитектур ==
| |
− | | |
− | == Существующие реализации ==
| |
− | | |
− | Реализация в проекте Vienna VL<ref>https://github.com/viennacl/viennacl-dev</ref>
| |
− | | |
− | Реализация в проекте IETL<ref>http://www.comp-phys.org/software/ietl/lanczos.html</ref>
| |
− | | |
− | Реализация в проекте ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>
| |
− | | |
− | Lanczos Plus Plus <ref>https://github.com/g1257/LanczosPlusPlus</ref>
| |
− | | |
− | ED Lanczos <ref>https://github.com/henhans/ED_lanczos</ref>
| |
− | | |
− | Paralleles Rechnen 2 <ref>https://github.com/soneyworld/ParallelesRechnen2</ref>
| |
− | | |
− | Cuda-Arnoldi <ref>https://github.com/trantalaiho/Cuda-Arnoldi/</ref>
| |
− | | |
− | = Литература =
| |
Шаблон:Assignmenta
Авторы: Хакимов А. С.
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 1950 году. Его эффективность обусловлена экономией памяти для хранения матриц и экономией вычислительных ресурсов. Алгоритм итерационный и использует степенной метод для поиска наибольших собственных значений и векторов матриц. Основной недостаток алгоритма заключается в накоплении ошибок округления, для решения которых появились методы поддержания ортогонализации т.н. векторов Ланцоша. Здесь мы рассмотрим выборочный метод поддержания ортогонализации, который существенно экономит процессорное время.
На вход алгоритма подаётся [math]A = A^T[/math],
- [math]
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] ,\, \;[/math]
случайный вектор [math]b[/math], как первое приближение собственного вектора матрицы и [math]k [/math] - количество собственных значений и собственных векторов, которые требуется найти.
Матрица [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n \times j[/math] строится на каждой итерации и состоит из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца [math]\theta_i [/math], - собственные значения симметричной трехдиагональной матрицы [math]T_j = Q^T_j A Q_j[/math] размерности [math]j \times j[/math].
- [math]
T_j = \begin{pmatrix}
\alpha_1 & \beta_1 \\
\beta_1 & \alpha_2 & \beta_2 \\
& \beta_2 & \ddots & \ddots \\
& & \ddots & \ddots & \beta_{j-1} \\
& & & \beta_{j-1} & \alpha_j
\end{pmatrix}\; (2).
[/math]
Однако, векторы [math]q_j [/math] теряют ортогональность вследствие приобретения больших компонент в направлениях векторов Ритца [math]y_{i,j} = Q_j v_i [/math], отвечающих сошедшимся числам Ритца [math] \theta_i [/math]. Поэтому чтобы построить [math]q_j [/math], предлагается на каждом шаге следить за оценками погрешностей [math]\beta_{t}|v_i(t)|, i = 1 \dots t, t = j - 1 [/math], где [math]v_i(t) [/math] - [math]t[/math]-я компонента собственного вектора [math]v_i [/math]. И когда какая-то оценка становится слишком малой, проводить ортогонализацию вектора Ланцоша [math]z [/math]. Величина [math]\beta_{t}|v_i(t)| [/math] считается малой, если она меньше, чем [math]\sqrt{\varepsilon}||T_{t}|| [/math], где [math]\varepsilon[/math] - доступная машинная точность чисел.
После следует вычисление собственных значений [math] \theta_j [/math] и собственных векторов [math]v_j [/math] полученной трехдиагональной матрицы [math]T_j[/math], например, с помощью метода "разделяй и властвуй"[1]
1.2 Математическое описание алгоритма
[math] q_{1} = b_{j}/\|b\|_2, \beta_0=0, q_0=0[/math]
[math] for\, j=1\,\, to\, \, k[/math]
[math]z=Aq_j,[/math]
[math]\alpha_j=q_j^Tz,[/math]
[math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},[/math]
/* Провести выборочную ортогонализацию по отношению
к сошедшимся векторам Ритца */
[math]for\, i \leqslant k, \,[/math]таких, что[math]
\lt math\gt z = z-(y^T_{i,k},z)y_{i,k}[/math]
[math]\text{end for}[/math]
[math]\beta_{j}=\|z\|_2[/math]
[math]q_{j+1}=z/\beta_{j}, [/math]
Вычислить собственные значения и собственные векторы
матрицы[math]\, \, T_{j} \, \,[/math]и оценки погрешности в них
[math]\text{end for}[/math]
- ↑ Дж. Деммель «Вычислительная линейная алгебра», c. 232, алгоритм 5.2