Участник: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>
 
= Литература =

Версия от 09:00, 19 января 2017

Шаблон:Assignmenta Авторы: Хакимов А. С.

Свойства и структура алгоритма

Общее описание алгоритма

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 1950 году. Его эффективность обусловлена экономией памяти для хранения матриц и экономией вычислительных ресурсов. Алгоритм итерационный и использует степенной метод для поиска наибольших собственных значений и векторов матриц. Основной недостаток алгоритма заключается в накоплении ошибок округления, для решения которых появились методы поддержания ортогонализации т.н. векторов Ланцоша. Здесь мы рассмотрим выборочный метод поддержания ортогонализации, который существенно экономит процессорное время.


На вход алгоритма подаётся [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{ ,\, \; }[/math]


случайный вектор [math]\displaystyle{ b }[/math], как первое приближение собственного вектора матрицы и [math]\displaystyle{ k }[/math] - количество собственных значений и собственных векторов, которые требуется найти.

Матрица [math]\displaystyle{ Q_j = [q_1, q_2, \dots, q_j] }[/math] размерности [math]\displaystyle{ n \times j }[/math] строится на каждой итерации и состоит из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца [math]\displaystyle{ \theta_i }[/math], - собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_j = Q^T_j A Q_j }[/math] размерности [math]\displaystyle{ j \times j }[/math].

[math]\displaystyle{ 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]\displaystyle{ q_j }[/math] теряют ортогональность вследствие приобретения больших компонент в направлениях векторов Ритца [math]\displaystyle{ y_{i,j} = Q_j v_i }[/math], отвечающих сошедшимся числам Ритца [math]\displaystyle{ \theta_i }[/math]. Поэтому чтобы построить [math]\displaystyle{ q_j }[/math], предлагается на каждом шаге следить за оценками погрешностей [math]\displaystyle{ \beta_{t}|v_i(t)|, i = 1 \dots t, t = j - 1 }[/math], где [math]\displaystyle{ v_i(t) }[/math] - [math]\displaystyle{ t }[/math]-я компонента собственного вектора [math]\displaystyle{ v_i }[/math]. И когда какая-то оценка становится слишком малой, проводить ортогонализацию вектора Ланцоша [math]\displaystyle{ z }[/math]. Величина [math]\displaystyle{ \beta_{t}|v_i(t)| }[/math] считается малой, если она меньше, чем [math]\displaystyle{ \sqrt{\varepsilon}||T_{t}|| }[/math], где [math]\displaystyle{ \varepsilon }[/math] - доступная машинная точность чисел.

После следует вычисление собственных значений [math]\displaystyle{ \theta_j }[/math] и собственных векторов [math]\displaystyle{ v_j }[/math] полученной трехдиагональной матрицы [math]\displaystyle{ T_j }[/math], например, с помощью метода "разделяй и властвуй"[1]

Математическое описание алгоритма

[math]\displaystyle{  q_{1} = b_{j}/\|b\|_2, \beta_0=0, q_0=0 }[/math]
[math]\displaystyle{  for\, j=1\,\, to\, \, k }[/math]
    [math]\displaystyle{ z=Aq_j, }[/math]
    [math]\displaystyle{ \alpha_j=q_j^Tz, }[/math]
    [math]\displaystyle{ z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, }[/math]
    /* Провести выборочную ортогонализацию по отношению
    к сошедшимся векторам Ритца */
    [math]\displaystyle{ for\, i \leqslant k, \, }[/math]таких, что[math]\displaystyle{ 
     \lt math\gt z = z-(y^T_{i,k},z)y_{i,k} }[/math]
    [math]\displaystyle{ \text{end for} }[/math]
    [math]\displaystyle{ \beta_{j}=\|z\|_2 }[/math]
    [math]\displaystyle{ q_{j+1}=z/\beta_{j},  }[/math]
    Вычислить собственные значения и собственные векторы
    матрицы[math]\displaystyle{ \, \, T_{j} \, \, }[/math]и оценки погрешности в них
[math]\displaystyle{ \text{end for} }[/math]
  1. ↑ Дж. Деммель «Вычислительная линейная алгебра», c. 232, алгоритм 5.2