Алговики - Вклад [ru]

Обсуждение участника:ArtyomKhakimov

2017-02-06T20:11:54Z

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

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил описание структуры алгоритма. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:42, 6 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил пояснения а графу. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 23:11, 6 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил оценки вычислительной мощности. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:48, 6 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Указал параметры компилцяии и запуска. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:56, 6 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Попробовал сделать выводы. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 23:06, 6 февраля 2017 (MSK)

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

2017-02-06T20:11:13Z

ArtyomKhakimov: /* Информационный граф */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

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

* Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
* Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math>
* При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Вариант с частичной переортогонализацией является промежуточным.
* Соотношение последовательной и параллельной сложности алгоритма равно <math>n</math>, что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш

= Программная реализация =

== Особенности реализации последовательного алгоритма ==

== Локальность данных и вычислений ==

== Возможные способы и особенности параллельной реализации алгоритма ==

== Масштабируемость алгоритма и его реализации ==

Исследование алгоритма на масштабируемость производилось на суперкомпьютере [http://parallel.ru/cluster "Ломоносов"] с использованием технологии MPI. Для тестирования была выбрана [http://pastebin.com/SbNrh7nz вот эта] реализация алгоритма.

Набор параметров, которые изменялись при тестировании:
*Число процессоров от 16 до 128.
*Размерность матрицы от 5000 до 50000 с шагом 5000.

Использовался компилятор mpicxx в режиме openmpi/1.5.5-gcc.
Компиляция осуществлялась с параметрами по-умолчанию.
Запуск вычислений производился на узле regular4.

{| 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
|}

По наклону графика при максимальной загрузке можно заметить, что на производительность алгоритма количество вычислительных процессов влияет сильнее, чем размерность входной матрицы. Наблюдается также ускорение при количестве параллельных процессоров 32 и размерности матрицы около 40 000 из чего можно сделать вывод, что такое сочетание является оптимальным для алгоритма Ланцоша.

[[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>

= Литература =

Обсуждение участника:ArtyomKhakimov

2017-02-06T20:06:32Z

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил описание структуры алгоритма. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:42, 6 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил оценки вычислительной мощности. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:48, 6 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Указал параметры компилцяии и запуска. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:56, 6 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Попробовал сделать выводы. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 23:06, 6 февраля 2017 (MSK)

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

2017-02-06T20:06:04Z

ArtyomKhakimov: /* Масштабируемость алгоритма и его реализации */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

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

* Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
* Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math>
* При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Вариант с частичной переортогонализацией является промежуточным.
* Соотношение последовательной и параллельной сложности алгоритма равно <math>n</math>, что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш

= Программная реализация =

== Особенности реализации последовательного алгоритма ==

== Локальность данных и вычислений ==

== Возможные способы и особенности параллельной реализации алгоритма ==

== Масштабируемость алгоритма и его реализации ==

Исследование алгоритма на масштабируемость производилось на суперкомпьютере [http://parallel.ru/cluster "Ломоносов"] с использованием технологии MPI. Для тестирования была выбрана [http://pastebin.com/SbNrh7nz вот эта] реализация алгоритма.

Набор параметров, которые изменялись при тестировании:
*Число процессоров от 16 до 128.
*Размерность матрицы от 5000 до 50000 с шагом 5000.

Использовался компилятор mpicxx в режиме openmpi/1.5.5-gcc.
Компиляция осуществлялась с параметрами по-умолчанию.
Запуск вычислений производился на узле regular4.

{| 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
|}

По наклону графика при максимальной загрузке можно заметить, что на производительность алгоритма количество вычислительных процессов влияет сильнее, чем размерность входной матрицы. Наблюдается также ускорение при количестве параллельных процессоров 32 и размерности матрицы около 40 000 из чего можно сделать вывод, что такое сочетание является оптимальным для алгоритма Ланцоша.

[[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>

= Литература =

Обсуждение участника:ArtyomKhakimov

2017-02-06T19:56:52Z

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил описание структуры алгоритма. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:42, 6 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил оценки вычислительной мощности. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:48, 6 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Указал параметры компилцяии и запуска. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:56, 6 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)

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

2017-02-06T19:56:19Z

ArtyomKhakimov: /* Масштабируемость алгоритма и его реализации */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

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

* Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
* Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math>
* При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Вариант с частичной переортогонализацией является промежуточным.
* Соотношение последовательной и параллельной сложности алгоритма равно <math>n</math>, что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш

= Программная реализация =

== Особенности реализации последовательного алгоритма ==

== Локальность данных и вычислений ==

== Возможные способы и особенности параллельной реализации алгоритма ==

== Масштабируемость алгоритма и его реализации ==

Исследование алгоритма на масштабируемость производилось на суперкомпьютере [http://parallel.ru/cluster "Ломоносов"] с использованием технологии MPI. Для тестирования была выбрана [http://pastebin.com/SbNrh7nz вот эта] реализация алгоритма.

Набор параметров, которые изменялись при тестировании:
*Число процессоров от 16 до 128.
*Размерность матрицы от 5000 до 50000 с шагом 5000.

Использовался компилятор mpicxx в режиме openmpi/1.5.5-gcc.
Компиляция осуществлялась с параметрами по-умолчанию.
Запуск вычислений производился на узле regular4.

{| 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>

= Литература =

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

2017-02-06T19:55:57Z

ArtyomKhakimov: /* Масштабируемость алгоритма и его реализации */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

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

* Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
* Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math>
* При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Вариант с частичной переортогонализацией является промежуточным.
* Соотношение последовательной и параллельной сложности алгоритма равно <math>n</math>, что говорит о том, что параллельная реализация теоретически должна давать ощутимый выигрыш

= Программная реализация =

== Особенности реализации последовательного алгоритма ==

== Локальность данных и вычислений ==

== Возможные способы и особенности параллельной реализации алгоритма ==

== Масштабируемость алгоритма и его реализации ==

Исследование алгоритма на масштабируемость производилось на суперкомпьютере [http://parallel.ru/cluster "Ломоносов"] с использованием технологии MPI. Для тестирования была выбрана [http://pastebin.com/SbNrh7nz вот эта] реализация алгоритма.

Набор параметров, которые изменялись при тестировании:
*Число процессоров от 16 до 128.
*Размерность матрицы от 5000 до 50000 с шагом 5000.

Использовался компилятор mpicxx в режиме openmpi/1.5.5-gcc
Компиляция осуществлялась с параметрами по-умолчанию.
Запуск вычислений производился на узле regular4.

{| 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>

= Литература =

Обсуждение участника:ArtyomKhakimov

2017-02-06T19:48:17Z

ArtyomKhakimov: /* Замечания по тексту */

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил описание структуры алгоритма. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:42, 6 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил оценки вычислительной мощности. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:48, 6 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)

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

2017-02-06T19:46:58Z

ArtyomKhakimov: /* Свойства алгоритма */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

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

* Алгоритм устойчив, так как операция переортогонализации как раз направлена на исключение проблемы округления вычисляемых значений
* Вычислительная мощность алгоритма - отношение числа операций к суммарному объему входных и выходных данных - примерно равно <math>2k</math>
* При реализации классического алгоритма Ланцоша возникает большая погрешность при округлении. Вариант с полной переортогонализацией позволяет избегать больших погрешностей, однако является более ресурсоемким. Вариант с частичной переортогонализацией является промежуточным.
* Соотношение последовательной и параллельной сложности алгоритма равно <math>n</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>

= Литература =

Обсуждение участника:ArtyomKhakimov

2017-02-06T19:42:18Z

ArtyomKhakimov: /* Замечания по тексту */

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Добавил список макроопераций. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:42, 6 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)

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

2017-02-06T19:41:51Z

ArtyomKhakimov: /* Макроструктура алгоритма */

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <math>z=Aq_j, </math>.

* Ортогонализация по отношению к сошедшимся векторам <math>z = z-(y^T_{i,k},z)y_{i,k} </math> для <math>i = 1 \dots t</math>

== Макроструктура алгоритма ==

Алгоритма Ланцоша представляется в виде совокупности метода Ланцоша, построения крыловского подпространства и процедуры Рэлея-Ритца.

Каждая итерация первого этапа может быть разложена на следующие макроединицы:
* Умножение матрицы на вектор
* Вычисление скалярного произведения векторов
* Умножение векторов на числа и вычитание векторов
* При выполнении условия необходимости переортогонализации – перемножение и вычитание векторов
* Вычисление нормы вектора
* Умножение вектора на число

Второй этап, процедура Рэлея-Ритца, заключается в нахождении собственных значений и собственных векторов матрицы и интерпретации их как приближенных значений собственных чисел и собственных векторов исходной матрицы, он может быть выполнен любым способом, например, QR-итерацией.

== Схема реализации последовательного алгоритма ==

Алгоритм выполняется в следующей последовательности:

<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>

= Литература =

Обсуждение участника:ArtyomKhakimov

2017-02-06T19:38:28Z

ArtyomKhakimov: /* Замечания по тексту */

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

=== Замечания по тексту ===

* В разделе 1.2 должен быть не псевдокод реализации, а описание алгоритма на математическом языке (в виде формул). [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
Учёл это замечание. --[[Участник:ArtyomKhakimov|ArtyomKhakimov]] ([[Обсуждение участника:ArtyomKhakimov|обсуждение]]) 22:38, 6 февраля 2017 (MSK)
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* К графу в разделе 1.7 требуются пояснения, в частности, при каких значениях внешних параметров он получен. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:19, 3 февраля 2017 (MSK)

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

2017-02-06T19:37:08Z

ArtyomKhakimov:

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>A</math> и начальный вектор <math>b</math>.

Алгоритм вычисляет собственные вектора матрицы <math>T_k</math>, которые являются столбцами матрицы <math>Q_k V</math>. Матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. При первой итерации вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math> нормируется. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>.

Затем начинаются итерации алгоритма. Всего итераций не более <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х.

Вычисления на <math>i</math>-й итерации таковы:
Сначала вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого подизводится вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Пусть норма полученного вектора - это элемент матрицы: <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то итерации заканчиваются. Далее будет следовать этап вычисления собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации.

После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref>

== Вычислительное ядро алгоритма ==

* Умножение матрицы на вектор, <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>

= Литература =

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

2017-01-26T11:15:48Z

ArtyomKhakimov:

{{Assignment}}
Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>end \, for</math>
<math>\beta_{j}=\|z\|_2</math>
<math>q_{j+1}=z/\beta_{j}, </math>
Вычислить собственные значения и собственные векторы
<math>\, \, T_{j} \, \, </math>и оценки погрешности в них
<math>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>

= Литература =

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

2017-01-20T08:24:50Z

ArtyomKhakimov: Содержимое страницы заменено на «{{Assignment}} Авторы: Хакимов А. С. [https://ru.wikipedia.org/wiki/Участник:ArtyomKhakimov/Алгорит…»

{{Assignment}}
Авторы: Хакимов А. С.

[https://ru.wikipedia.org/wiki/Участник:ArtyomKhakimov/Алгоритм_Ланцоша_для_арифметики_с_плавающей_точкой_с_выборочной_ортогонализацией Статья временно находится на Википедии]

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

2017-01-19T17:39:15Z

ArtyomKhakimov:

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>\text{/* Провести выборочную ортогонализацию по отношению}, </math>
<math>\text{ к сошедшимся векторам Ритца */}</math>
<math>for\, i \leqslant k, \text{таких, что} \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math>
<math>z = z-(y^T_{i,k},z)y_{i,k}</math>
<math>end \, for</math>
<math>\beta_{j}=\|z\|_2</math>
<math>q_{j+1}=z/\beta_{j}, </math>
<math>\text{Вычислить собственные значения и собственные векторы}</math>
<math>\text{матрицы} \, \, T_{j} \, \text{и оценки погрешности в них}</math>
<math>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>

= Литература =

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

2017-01-19T12:26:01Z

ArtyomKhakimov:

{{Assignment}}
Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>\text{/* Провести выборочную ортогонализацию по отношению}, </math>
<math>\text{ к сошедшимся векторам Ритца */}</math>
<math>for\, i \leqslant k, \text{таких, что} \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math>
<math>z = z-(y^T_{i,k},z)y_{i,k}</math>
<math>end \, for</math>
<math>\beta_{j}=\|z\|_2</math>
<math>q_{j+1}=z/\beta_{j}, </math>
<math>\text{Вычислить собственные значения и собственные векторы}</math>
<math>\text{матрицы} \, \, T_{j} \, \text{и оценки погрешности в них}</math>
<math>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>

= Литература =

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

2017-01-19T11:06:01Z

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

2017-01-19T10:41:36Z

ArtyomKhakimov:

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>\text{/* Провести выборочную ортогонализацию по отношению}, </math>
<math>\text{ к сошедшимся векторам Ритца */}</math>
<math>for\, i \leqslant k, \text{таких, что} \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math>
<math>z = z-(y^T_{i,k},z)y_{i,k}</math>
<math>end \, for</math>
<math>\beta_{j}=\|z\|_2</math>
<math>q_{j+1}=z/\beta_{j}, </math>
<math>\text{Вычислить собственные значения и собственные векторы}</math>
<math>\text{матрицы} \, \, T_{j} \, \text{и оценки погрешности в них}</math>
<math>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>

= Литература =

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

2017-01-19T09:45:05Z

ArtyomKhakimov:

{{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 Статья временно находится на Википедии]

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

2017-01-19T09:41:24Z

ArtyomKhakimov:

{{Assignmenta}}
Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>

= Литература =

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

2017-01-19T09:17:18Z

ArtyomKhakimov:

{{Assignmenta}}
Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>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>

= Литература =

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

2017-01-19T09:00:46Z

ArtyomKhakimov:

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

2017-01-19T08:57:04Z

ArtyomKhakimov:

Авторы: Хакимов А. С.

= Свойства и структура алгоритма =
== Общее описание алгоритма ==

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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>\text{/* Провести выборочную ортогонализацию по отношению}, </math>
<math>\text{ к сошедшимся векторам Ритца */}</math>
<math>for\, i \leqslant k, \text{таких, что} \, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{k}\| \, </math>
<math>z = z-(y^T_{i,k},z)y_{i,k}</math>
<math>end \, for</math>
<math>\beta_{j}=\|z\|_2</math>
<math>q_{j+1}=z/\beta_{j}, </math>
<math>\text{Вычислить собственные значения и собственные векторы}</math>
<math>\text{матрицы} \, \, T_{j} \, \text{и оценки погрешности в них}</math>
<math>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>

= Литература =

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

2017-01-19T08:32:51Z

ArtyomKhakimov:

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

2017-01-19T08:28:30Z

ArtyomKhakimov:

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

2017-01-19T08:26:40Z

ArtyomKhakimov:

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

2017-01-19T08:26:21Z

ArtyomKhakimov:

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

2017-01-18T18:30:23Z

ArtyomKhakimov: Новая страница: «Авторы: Хакимов А. С. == Свойства и структура алгоритма == === Общее описание алгоритма === Ал…»

Авторы: Хакимов А. С.

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===

Алгоритм Ланцоша служит для нахождения собственных значений и собственных векторов для больших разреженных матриц, к которым нельзя применить прямые методы из-за больших требований к памяти и времени. Он был опубликован Корнелием Ланцошем в 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\, \, do:</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>/* Провести выборочную ортогонализацию по отношению, </math>
<math> к сошедшимся векторам Ритца */, </math>
<math>t = j - 1, </math>
<math>for\, i=1\,\, to\, \, t\, \, do: </math>
<math>if\, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, then:</math>
<math>z = z-(y^T_{i,t},z)y_{i,t} </math>, где <math>y_{i,t} = Q_{t}v_i</math>
<math>\beta_{j}=\|z\|_2 </math>
<math>q_{j+1}=z/\beta_{j}, </math>
Строим матрицу <math> T_j</math> (2), вычисляем собственные значения <math> \theta_j </math> и собственные векторы <math>v_j </math> полученной матрицы <math>T_j</math>.