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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-08T12:07:13Z

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

{{Assignment|ASA}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''': А.В. Шостик.

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследовалась реализация<ref>http://pastebin.com/3W4uUVGj</ref> на языке C++ с использованием MPI. Компилятор GCC 5.2.1; версия MPI 1.8.4;

Для компиляции и запуска использовались следующие команды:
<pre>$ mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++ -o <EXE_FILE>
$ mpirun -np <N_PROCESSORS> <EXE_FILE>
</pre>

* Число процессоров: 2, 4, 8, 16, 32, 64, 128.
* Размерность матрицы <math>N</math> от 5 000 до 50 000, шаг 5 000.
* Количество итераций алгоритма Ланцоша для всех экспериментов 200.

Результаты:

{| class="wikitable"
|+ Время выполнения алгоритма
|-
!
! 2
! 4
! 8
! 16
! 32
! 48
! 64
! 96
! 128
|-
! 5000
| 1,5243
| 1,0951
| 0,9334
| 0,6631
| 0,7456
| 0,7038
| 0,7989
| 0,9985
| 0,985
|-
! 10 000
| 5,667
| 4,2625
| 2,7502
| 3,6762
| 3,4435
| 3,6192
| 3,5194
| 3,8486
| 3,2956
|-
! 15 000
| 12,5207
| 6,5995
| 7,2897
| 8,5086
| 8,539
| 9,0018
| 8,9427
| 15,0132
| 7,7088
|-
! 20 000
| 23,8143
| 15,851
| 10,6313
| 13,7241
| 12,8051
| 14,3183
| 13,2946
| 14,0864
| 13,7499
|-
! 25 000
| 37,5419
| 17,9436
| 15,5451
| 18,3944
| 19,6076
| 18,5444
| 19,1347
| 21,8366
| 21,4454
|-
! 30 000
| 55,014
| 26,2702
| 22,3943
| 26,4935
| 25,5618
| 29,6784
| 28,0853
| 27,8872
| 31,2415
|-
! 35 000
| 71,3969
| 37,4811
| 30,2274
| 34,0041
| 38,1028
| 38,4912
| 38,0174
| 37,1407
| 40,4463
|-
! 40 000
| 100,326
| 47,4248
| 40,2295
| 44,0332
| 50,0015
| 52,9428
| 46,2961
| 51,4372
| 49,7504
|-
! 45 000
| 142,93
| 66,4829
| 52,1382
| 55,4683
| 62,3042
| 65,0885
| 60,0512
| 57,8891
| 64,3764
|}

Ниже приведен график зависимости времени выполнения программы от числа процессоров и размера исходной матрицы.

[[File:Lanzcos.png|600px|thumb|center|Рисунок 2. Время работы программы в зависимости от количества процессоров и размера матрицы. 

<math>Procs</math> — количество процессоров, 
<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Time</math> — время работы алгоритма в сек. ]]

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-08T03:08:21Z

Shostix:

{{Assignment|ASA}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''': А.В. Шостик.

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследовалась реализация<ref>http://pastebin.com/3W4uUVGj</ref> на языке C++ с использованием MPI. Компилятор GCC 5.2.1; версия MPI 1.8.4;

Для компиляции и запуска использовались следующие команды:
<pre>$ mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++ -o <EXE_FILE>
$ mpirun -np <N_PROCESSORS> <EXE_FILE>
</pre>

* Число процессоров: 2, 4, 8, 16, 32, 64, 128.
* Размерность матрицы <math>N</math> от 5 000 до 50 000, шаг 5 000.
* Количество итераций алгоритма Ланцоша для всех экспериментов 200.

Ниже приведена зависимость времени выполнения программы от числа процессоров и размера исходной матрицы.

[[|600px|thumb|center|Рисунок 2. Эффективность параллельной реализации алгоритма Ланцоша. 

<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Cores</math> — количество процессоров, 
<math>Eff</math> — эффективность алгоритма в %. ]]

'''Выводы''':

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-08T03:04:23Z

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

{{Assignment|ASA}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''': А.В. Шостик.

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследовалась реализация<ref>http://pastebin.com/3W4uUVGj</ref> на языке C++ с использованием MPI. Компилятор GCC 5.2.1; версия MPI 1.8.4;

Для компиляции и запуска использовались следующие команды:
<pre>$ mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++ -o <EXE_FILE>
$ mpirun -np <N_PROCESSORS> <EXE_FILE>
</pre>

* Число процессоров: 2, 4, 8, 16, 32, 64, 128.
* Размерность матрицы <math>N</math> от 5 000 до 50 000, шаг 5 000.
* Количество итераций алгоритма Ланцоша для всех экспериментов 200.

Ниже приведена зависимость времени выполнения программы от числа процессоров и размера исходной матрицы.

[[File:Lanczos Without Orthogonalization Eff2.png|600px|thumb|center|Рисунок 2. Эффективность параллельной реализации алгоритма Ланцоша. 

<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Cores</math> — количество процессоров, 
<math>Eff</math> — эффективность алгоритма в %. ]]

'''Выводы''':
* минимальная эффективность алгоритма при указанных параметрах составила 1.9% (при числе процессоров 32, размерности матрицы 5 000),
* максимальная эффективность алгоритма при указанных параметрах составила 59.5% (при числе процессоров 256, размерности матрицы 200 000),
* средняя эффективность алгоритма при числе процессоров от 32 до 256 и размерности матрицы от 25 000 до 200 000 составила 52,7%,
* при размерности матрицы менее 25 000 эффективность алгоритма менее 23%.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-02T21:05:44Z

Shostix:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''': А.В. Шостик.

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследуемая реализация на языке C [https://algowiki-project.org/ru/%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:A.Freeman/%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_%D1%82%D0%BE%D1%87%D0%BD%D0%BE%D0%B9_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B8_%28%D0%B1%D0%B5%D0%B7_%D0%BF%D0%B5%D1%80%D0%B5%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%B8%29#.D0.A1.D1.83.D1.89.D0.B5.D1.81.D1.82.D0.B2.D1.83.D1.8E.D1.89.D0.B8.D0.B5_.D1.80.D0.B5.D0.B0.D0.BB.D0.B8.D0.B7.D0.B0.D1.86.D0.B8.D0.B8] <ref>https://github.com/xrstalker/lanczos</ref>.

Компилятор gcc, аргумент -std=c11. Openmpi-1.8.4

* Число процессоров варьировалось от 32 до 256, step = 16.
* Размер матрицы изменялся от <math>5 000 \times 5 000</math> до <math>200 000 \times 200 000</math>, step = 5 000.

Ниже приведена зависимость эффективности параллельной реализации от указанных выше параметров.

[[File:Lanczos Without Orthogonalization Eff2.png|600px|thumb|center|Рисунок 2. Эффективность параллельной реализации алгоритма Ланцоша. 

<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Cores</math> — количество процессоров, 
<math>Eff</math> — эффективность алгоритма в %. ]]

'''Выводы''':
* минимальная эффективность алгоритма при указанных параметрах составила 1.9% (при числе процессоров 32, размерности матрицы 5 000),
* максимальная эффективность алгоритма при указанных параметрах составила 59.5% (при числе процессоров 256, размерности матрицы 200 000),
* средняя эффективность алгоритма при числе процессоров от 32 до 256 и размерности матрицы от 25 000 до 200 000 составила 52,7%,
* при размерности матрицы менее 25 000 эффективность алгоритма менее 23%.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-02T20:54:19Z

Shostix:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследуемая реализация на языке C [https://algowiki-project.org/ru/%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:A.Freeman/%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_%D1%82%D0%BE%D1%87%D0%BD%D0%BE%D0%B9_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B8_%28%D0%B1%D0%B5%D0%B7_%D0%BF%D0%B5%D1%80%D0%B5%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%B8%29#.D0.A1.D1.83.D1.89.D0.B5.D1.81.D1.82.D0.B2.D1.83.D1.8E.D1.89.D0.B8.D0.B5_.D1.80.D0.B5.D0.B0.D0.BB.D0.B8.D0.B7.D0.B0.D1.86.D0.B8.D0.B8] <ref>https://github.com/xrstalker/lanczos</ref>.

Компилятор gcc, аргумент -std=c11. Openmpi-1.8.4

* Число процессоров варьировалось от 32 до 256, step = 16.
* Размер матрицы изменялся от <math>5 000 \times 5 000</math> до <math>200 000 \times 200 000</math>, step = 5 000.

Ниже приведена зависимость эффективности параллельной реализации от указанных выше параметров.

[[File:Lanczos Without Orthogonalization Eff2.png|600px|thumb|center|Рисунок 2. Эффективность параллельной реализации алгоритма Ланцоша. 

<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Cores</math> — количество процессоров, 
<math>Eff</math> — эффективность алгоритма в %. ]]

'''Выводы''':
* минимальная эффективность алгоритма при указанных параметрах составила 1.9% (при числе процессоров 32, размерности матрицы 5 000),
* максимальная эффективность алгоритма при указанных параметрах составила 59.5% (при числе процессоров 256, размерности матрицы 200 000),
* средняя эффективность алгоритма при числе процессоров от 32 до 256 и размерности матрицы от 25 000 до 200 000 составила 52,7%,
* при размерности матрицы менее 25 000 эффективность алгоритма менее 23%.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализации алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-02T20:53:47Z

Shostix:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследуемая реализация на языке C [https://algowiki-project.org/ru/%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:A.Freeman/%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_%D1%82%D0%BE%D1%87%D0%BD%D0%BE%D0%B9_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B8_%28%D0%B1%D0%B5%D0%B7_%D0%BF%D0%B5%D1%80%D0%B5%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%B8%29#.D0.A1.D1.83.D1.89.D0.B5.D1.81.D1.82.D0.B2.D1.83.D1.8E.D1.89.D0.B8.D0.B5_.D1.80.D0.B5.D0.B0.D0.BB.D0.B8.D0.B7.D0.B0.D1.86.D0.B8.D0.B8] <ref>https://github.com/xrstalker/lanczos</ref>.

Компилятор gcc, аргумент -std=c11. Openmpi-1.8.4

* Число процессоров варьировалось от 32 до 256, step = 16.
* Размер матрицы изменялся от <math>5 000 \times 5 000</math> до <math>200 000 \times 200 000</math>, step = 5 000.

Ниже приведена зависимость эффективности параллельной реализации от указанных выше параметров.

[[File:Lanczos Without Orthogonalization Eff2.png|600px|thumb|center|Рисунок 2. Эффективность параллельной реализации алгоритма Ланцоша. 

<math>Size</math> — размерность матрицы <math>N</math>, 
<math>Cores</math> — количество процессоров, 
<math>Eff</math> — эффективность алгоритма в %. ]]

'''Выводы''':
* минимальная эффективность алгоритма при указанных параметрах составила 1.9% (при числе процессоров 32, размерности матрицы 5 000),
* максимальная эффективность алгоритма при указанных параметрах составила 59.5% (при числе процессоров 256, размерности матрицы 200 000),
* средняя эффективность алгоритма при числе процессоров от 32 до 256 и размерности матрицы от 25 000 до 200 000 составила 52,7%,
* при размерности матрицы менее 25 000 эффективность алгоритма менее 23%.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.
* Lanczos Plus Plus (2009) <ref>https://github.com/g1257/LanczosPlusPlus</ref>. Пример реализация алгоритма Ланцоша в применении к задаче расчета взаимодействия сильно кореллированных электронов. С частичной переортогонализацией, с полной переортогонализацией. Язык C++.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-02T20:14:00Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

Исследование проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster/].

Исследуемая реализация на языке C [https://algowiki-project.org/ru/%D0%A3%D1%87%D0%B0%D1%81%D1%82%D0%BD%D0%B8%D0%BA:A.Freeman/%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_%D1%82%D0%BE%D1%87%D0%BD%D0%BE%D0%B9_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B8_%28%D0%B1%D0%B5%D0%B7_%D0%BF%D0%B5%D1%80%D0%B5%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%B8%29#.D0.A1.D1.83.D1.89.D0.B5.D1.81.D1.82.D0.B2.D1.83.D1.8E.D1.89.D0.B8.D0.B5_.D1.80.D0.B5.D0.B0.D0.BB.D0.B8.D0.B7.D0.B0.D1.86.D0.B8.D0.B8] <ref>https://github.com/xrstalker/lanczos</ref>.

Компилятор gcc, аргумент -std=c11. Openmpi-1.8.4

* Число процессоров варьировалось от 32 до 256, step = 16.
* Размер матрицы изменялся от <math>5 000 \times 5 000</math> до <math>200 000 \times 200 000</math>, step = 5 000.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref><ref>http://www.grycap.upv.es/slepc</ref> проекты на базе данного алгоритма в хронологическом порядке.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует реализация с распараллеливанием на базе OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>Hernandez V., Roman J. E., Vidal V. SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems //ACM Transactions on Mathematical Software (TOMS).– Т. 31. – №. 3. – С. 351-362.</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C. Присутствует реализация с распараллеливанием на базе MPI.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-02-02T19:18:25Z

Shostix: /* Существующие реализации алгоритма */

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
|=== Существующие реализации алгоритма ===

Алгоритм Ланцоша для точной арифметики без переортогонализации на практике используется гораздо реже, чем алгоритм Ланцоша с частичной или полной ортогонализацией, но, тем не менее, существует множество реализаций данного алгоритма, как последовательных, так и параллельных. Ниже перечислены наиболее известные и обозреваемые в литературе <ref>Hernandez V. et al. A survey of software for sparse eigenvalue problems //Universidad Politecnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-6, http://www. grycap. upv. es/slepc. – 2009.</ref><ref>Hernandez V. et al. Lanczos methods in SLEPc //Universidad Politécnica de Valencia, Valencia, Spain, SLEPc Technical Report STR-5. – 2006. – С. 136.</ref> проекты на базе данного алгоритма в хронологическом порядке.

* UNDERWOOD (1977) <ref>Golub, G. H. and R. Underwood. The Block Lanczos Method for Computing Eigenvalues. In Mathematical Software III (edited by J. Rice), pp. 364{377. Academic Press, New York, NY, USA</ref>. Без переортогонализации и с частичной реортогонализацией. Язык Fortran 77.
* PROPACK (1998) <ref>http://www.daimi.au.dk/PB/537</ref>. Частичная реортогонализация. Язык Fortran 77, Matlab. Присутствует параллельная реализация OpenMP.
* BLKLAN (2004) <ref>Liu, G., W. Xu, and S. Qiao. Block Lanczos Tridiagonalization of Complex Symmetric Matrices. Technical Report CAS-04-07-SQ, Department of Computing and Software, McMaster University, Hamilton, Ontario, Canada</ref>. Без реортогонализации, с частичной реортогонализацией. Язык C, Matlab.
* IETL (2006) <ref>http://www.comp-phys.org/software/ietl/</ref>. Только без реортогонализации. Язык C++.
* SLEPc (2009) <ref>http://www.comp-phys.org/software/ietl/</ref>. Без реортогонализации, с выборочной реортогонализацией, с частичной реортогонализацией, с полной реортогонализацией. Язык C, Fortran 77. Присутствует реализация с распараллеливанием на базе MPI.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-24T08:34:49Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[File:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-24T08:29:28Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>\frac{n(n+1)}{2} + 1</math>
| output_data = <math>kn+k</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений и собственных векторов симметричных матриц (задача <math>Ax = \lambda x</math>). Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений неприменимы из-за высокой вычислительной сложности. На практике чаще всего используется для работы с большими разреженными матрицами.

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

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). К сожалению, на практике использование метода Ланцоша без ортогонализации затруднено ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. Такая проблема решается усовершенствованием алгоритма Ланшоца (алгоритмами Ланшоца с выборочной и полной ортогонализацией)<ref>Деммель Д. Вычислительная линейная алгебра</ref>.

В данной статье будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

* симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, элементы матрицы <math>a_{ij}=a_{ji}</math>
* начальный вектор <math>v = {v_1, v_2, ..., v_n}, v \neq \theta</math>
Так как матрица является симметрической, достаточно хранить <math>\frac{n(n+1)}{2}</math> ее элементов.

Для исходной матрицы <math>A</math> и вектора <math>v</math> строится крыловское подпространство. Крыловское подпространство - это линейная оболочка векторов <math>[v, Av, A^2v, ..., A^{k-1}v]</math>, называемых векторами Ланцоша.
Каждый из векторов Ланцоша ортонормируется: <math>q_k=\frac{A^kv}{\|A^kv\|}</math>, и на шаге <math>k</math> алгоритма имеется Крыловская матрица <math>Q = {q_0, ..., q_{k-1}}</math> размерности <math>n \times k</math>.

Соответствующий текущей итерации <math>k</math> вектор Ланцоша считается <math>k</math>-ым приближением собственного вектора исходной матрицы <math>A</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> (числа Ритца). В качестве алгоритма нахождения собственных векторов и собственных значений симметрической трехдиагональной матрицы будем использовать метод "разделяй и властвуй", который требует <math>O(k^3)</math> операций.<ref>https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm</ref>.

''Выходные данные'':

* совокупность собственных векторов матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных векторов матрицы <math>A</math>)
* совокупность собственных значений матрицы <math>T_k</math> (т.е. <math>k</math>-ых приближений собственных значений матрицы <math>A</math>)

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

'''Вычислительное ядро'''. т.е. часть алгоритма, требующая наибольших вычислительных затрат на каждой итерации состоит из вычисления промежуточного вектора <math>z=Aq_i</math>

=== Макроструктура алгоритма ===
''' Итерация алгоритма: '''
* вычисление нормы вектора,
* деление вектора на число,
* умножение матрицы на вектор,
* линейная комбинация векторов (умножение на число и сложение).

''' Финальный расчет: '''
* вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.

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

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

'''Вход''' : размерность матрицы <math>n</math>, элементы симметричной матрицы <math>A</math>, количество итераций <math>k</math>
'''Выход''' : собственные вектора матрицы <math>T_k</math>, матрица собственных значений
<math>
\begin{align}
q_1 = & b/ \|b\|_2,\;
\beta_0 = 0,\;
q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align};
</math>

<math>procedure(T_k)</math>;

<math>procedure()</math> - процедура вычисления собственных векторов и собственных значений симметрической трехдиагональной матрицы.

=== Последовательная сложность алгоритма ===

Глобально алгоритм состоит из инициализации, <math>k</math> итераций и финального расчета собственных векторов и собственных значений.

Рассмотрим сложность каждой макрооперации алгоритма:
* вычисление нормы вектора - <math>(n-1)</math> сложений, <math>n</math> умножений, операция извлечения корня
* деление вектора на число - <math>n</math> операций
* умножение матрицы на вектор - <math>n^2</math> операций
* скалярное произведение векторов - <math>n</math> операций
* линейная комбинация векторов - <math>2n</math> умножений и вычитаний
* вычисление собственных векторов и собственных значений матрицы <math>T_k</math> методом "Разделяй и властвуй" - <math>O(k^3)</math> операций.

'''Инициализация:'''
* вычисление нормы вектора,
* деление вектора на число.
'''Одна итерация:'''
* умножение матрицы на вектор,
* скалярное произведение векторов,
* линейная комбинация векторов,
* вычисление нормы вектора,
* деление вектора на число.
'''Финальный расчет:'''
* вычисление собственных векторов и собственных значений матрицы

Итого '''последовательная сложность алгоритма''' составляет <math>O(3n+k(n^2+n+2n+2n+3n))+O(k^3)</math>.

Так как на практике порядок матрицы <math>n</math> намного превышает количество итераций <math>k</math>, сложность алгоритма Ланцоша можно рассматривать <math>O(kn^2)</math>.

=== Информационный граф ===

Информационный граф алгоритма с входными и выходными данными можно описать в виде рисунка:

[[Файл:Ланцош.png|600px|thumb|center|Рисунок 1. Информационный граф алгоритма Ланцоша. 
<math>\|.\|</math> — операция вычисления нормы, 
<math>/</math> — операция деления, 
<math>A</math> — операция умножения матрицы на вектор, 
<math>*</math> — операция скалярного произведения, 
<math>L</math> — операция вычисления линейной комбинации векторов, 
▽ - проверка условия <math>\mathsf{(\beta_i = 0)}</math> и выход из цикла в случае выполнения условия]]

=== Ресурс параллелизма алгоритма ===

Внутри каждой итерации алгоритма можно задействовать параллельные вычисления для суммирования векторов (использовать суммирование методом сдваивания элементов). Эта операция будет использована в макрооперации умножения матрицы на вектор, которая требует в последовательной реализации <math>n</math> умножений и <math>n-1</math> сложение. В таком виде сложение <math>n</math> элементов можно выполнить за <math>\log_2 n</math> действий.

=== Входные и выходные данные алгоритма ===

''Входные данные'': симметрическая матрица <math>A=A^T</math> размерности <math>n</math>, количество итераций <math>k</math>.

''Объем входных данных'': <math>\frac{n(n+1)}{2} + 1</math>.

''Выходные данные'': <math>k</math> собственных векторов, собственных значений матрицы на <math>k</math>-ом приближении.

''Объем выходных данных'': <math>kn+k</math>.

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

* Сложность последовательной реализации алгоритма <math>O(kn^2)</math>.
* Сложность параллельного алгоритма Ланцоша по высоте ЯПФ <math>O(k\log n)</math>, по ширине ЯПФ <math>O(n^2)</math>.
** В умножении матрицы на вектор сложение n элементов выполняется в <math>\log_2 n</math> ярусов шириной <math>\frac{n}{2}, \frac{n}{4}, ... , 1</math>, остальные операции внутри итерации выполняются последовательно
* Отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>.
* Вычислительная мощность алгоритма Ланцоша из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. Учитывая <math>k</math> много меньше <math>n</math>, вычислительная мощность ≈ <math> 2k</math>.
* В алгоритме Ланцоша возможно выполнение числа итераций менее <math>k</math> если все собственные значения матрицы вычисляются раньше.
* Для найденных собственных значений справедливо: <math>\lambda_i(T_{k+1})\geq \lambda_i(T_{k})\geq \lambda_{i+1}(T_{k+1})\geq \lambda_{i+1}(T_{k})</math>, т.е. в первую очередь находятся максимальные по модулю собственные значения.
* Как следствие предыдущего пункта алгоритм Ланцоша особенно удобен для вычисления собственных значений матрицы, находящихся на границе ее спектра.

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===

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

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

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-19T00:52:52Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>n*(n-1)+2</math>
| output_data = <math>k*(n+1)</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений симметричных матриц. Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений наприменимы из-за высокой вычислительной сложности.

Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на методе Ланцоша для построения крыловского подпространства и процедуре Рэлея-Ритца.

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). Он не является устойчивым к влиянию ошибок округления, в отличие от своих модификаций - алгоритма Ланцоша с выборочной ортогонализацией и полной ортогонализацией. На практике часто применяется алгоритм Ланцоша с полной ортогонализацией, предложенный в 1970 году. Он является оптимальным по точности вычислений, хоть и наиболее дорогостоящим алгоритмом из указанных.

Для наименее дорогостоящего в вычислениях алгоритма Ланцоша без ортогонализации будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-19T00:48:45Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>n*(n+1)/2</math>
| output_data = <math>k*(n+1)</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений симметричных матриц. Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений наприменимы из-за высокой вычислительной сложности.

Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на методе Ланцоша для построения крыловского подпространства и процедуре Рэлея-Ритца.

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). Он не является устойчивым к влиянию ошибок округления, в отличие от своих модификаций - алгоритма Ланцоша с выборочной ортогонализацией и полной ортогонализацией. На практике часто применяется алгоритм Ланцоша с полной ортогонализацией, предложенный в 1970 году. Он является оптимальным по точности вычислений, хоть и наиболее дорогостоящим алгоритмом из указанных.

Для наименее дорогостоящего в вычислениях алгоритма Ланцоша без ортогонализации будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

''Исходные данные'':

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-19T00:08:43Z

Shostix:

{{algorithm
| name = Алгоритм Ланцоша для точной арифметики (без ортогонализации)
| serial_complexity = <math>O(kn^2)</math>
| input_data = <math>n</math>
| output_data = <math>k*(n+1)</math>
| pf_height = <math>O(k*log(n))</math>
| pf_width = <math>O(n^2)</math>
}}

'''Авторы''':
* [[Участник:Shostix|А.В. Шостик]].

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

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

'''Алгоритм Ланцоша''' - один из итерационных методов вычисления собственных значений симметричных матриц. Используется для работы с матрицами столь большими, что классические численные методы нахождения собственных значений наприменимы из-за высокой вычислительной сложности.

Алгоритм был предложен в 1950 году венгерским физиком и математиком Корнелием Ланцошем. Основан на методе Ланцоша для построения крыловского подпространства и процедуре Рэлея-Ритца.

В данной статье будет рассмотрен простой метод Ланцоша (без ортогонализации). Он не является устойчивым к влиянию ошибок округления, в отличие от своих модификаций - алгоритма Ланцоша с выборочной ортогонализацией и полной ортогонализацией. На практике часто применяется алгоритм Ланцоша с полной ортогонализацией, предложенный в 1970 году. Он является оптимальным по точности вычислений, хоть и наиболее дорогостоящим алгоритмом из указанных.

Для наименее дорогостоящего в вычислениях алгоритма Ланцоша без ортогонализации будем подразумевать отсутствие влияния ошибок округления на вычислительный процесс.

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-18T23:50:17Z

Shostix:

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-18T23:44:15Z

Shostix:

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

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

Участник:Shostix/Алгоритм Ланцоша для точной арифметики (без переортогонализации)

2017-01-18T23:38:13Z

Shostix: Новая страница: «{{algorithm | name = Алгоритм Ланцоша для точной арифметики (без ортогонализации) | serial_complexity = <math>…»