Участник:ArtyomKhakimov/Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 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>
 
 
 
= Литература =
 

Версия 12:00, 19 января 2017

Шаблон: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]

  1. Дж. Деммель «Вычислительная линейная алгебра», c. 232, алгоритм 5.2