|
|
| Строка 1: |
Строка 1: |
| {{Assignmenta}} | | {{Assignment}} |
| Авторы: Хакимов А. С. | | Авторы: Хакимов А. С. |
|
| |
|
| = Свойства и структура алгоритма =
| | [https://ru.wikipedia.org/wiki/%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:ArtyomKhakimov/%D0%90%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88%D0%B0_%D0%B4%D0%BB%D1%8F_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B8_%D1%81_%D0%BF%D0%BB%D0%B0%D0%B2%D0%B0%D1%8E%D1%89%D0%B5%D0%B9_%D1%82%D0%BE%D1%87%D0%BA%D0%BE%D0%B9_%D1%81_%D0%B2%D1%8B%D0%B1%D0%BE%D1%80%D0%BE%D1%87%D0%BD%D0%BE%D0%B9_%D0%BE%D1%80%D1%82%D0%BE%D0%B3%D0%BE%D0%BD%D0%B0%D0%BB%D0%B8%D0%B7%D0%B0%D1%86%D0%B8%D0%B5%D0%B9 Статья временно находится на Википедии] |
| == Общее описание алгоритма ==
| |
| | |
| Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 1950 году. Его эффективность обусловлена экономией памяти для хранения матриц и экономией вычислительных ресурсов. Алгоритм итерационный и использует степенной метод для поиска наибольших собственных значений и векторов матриц. Основной недостаток алгоритма заключается в накоплении ошибок округления, для решения которых появились методы поддержания ортогонализации т.н. векторов Ланцоша. Здесь мы рассмотрим выборочный метод поддержания ортогонализации, который существенно экономит процессорное время.
| |
| | |
| | |
| На вход алгоритма подаётся <math>A = A^T</math>,
| |
| | |
| | |
| {{Шаблон:ASymmetric}} <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>, например, с помощью метода "разделяй и властвуй"<ref>Дж. Деммель «Вычислительная линейная алгебра», c. 232, алгоритм 5.2</ref>
| |
| | |
| == Математическое описание алгоритма ==
| |
| | |
| <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 \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math>
| |
| <math>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>
| |
| | |
| == Вычислительное ядро алгоритма ==
| |
| | |
| * Умножение матрицы на вектор, <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|рис. 2: Время работы алгоритма]]
| |
| | |
| == Динамические характеристики и эффективность реализации алгоритма ==
| |
| | |
| == Выводы для классов архитектур ==
| |
| | |
| == Существующие реализации ==
| |
| | |
| Реализация в проекте 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>
| |
| | |
| = Литература =
| |
