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

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-24T09:56:03Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|'''Рисунок 1.''' Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|'''Рисунок 2.''' Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan] на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]
Набор и границы значений изменяемых параметров реализации алгоритма:
* размер матрицы от 5000 до 25000 с шагом 500.
* число процессоров 1, 2, 4, [8:144] с шагом 8.
Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:
{|align="center"
|-valign="top"
|[[file:Trlanlib_mpi1.png|thumb|center|700px|'''Рисунок 3.''' Производительность в зависимости от числа процессоров и размера матрицы.]]
|[[file:Trlanlib_mpi2.png|thumb|center|700px|'''Рисунок 4.''' Время выполнения в зависимости от числа процессоров и размера матрицы.]]
|}
Можно заметить, что данная реализация немасштабируема и не дает ускорения при увеличении числа процессов. Из рисунков 3, 4 хорошо видно, что наилучшая производительность достигается с использованием только одного-двух процессов. Дальнейшее увеличение числа процессов лишь приводит к снижению эффективности. После 8-ми процессов ситуация стабилизируется и дальнейшее распараллеливание мало сказывается на результат
{| class="wikitable mw-collapsible mw-collapsed"
! Использованные модули и компиляторы
|-
|
Модули:
mkl/4.0.2.146
impi/5.0.1-ofa
компилятор mpif90, параметры компиляции:
mpif90 -I/opt/intel/composer_xe_2013.2.146/mkl/include/ia32 -I/opt/intel/composer_xe_2013.2.146/mkl/include -O0 /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_lapack95_lp64.a -L/opt/intel/composer_xe_2013.2.146/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl simple.f90 -o simple lapack_aux.o ../../libtrlan_mpi.a

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

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*библиотека [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-23T12:30:04Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan] на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]
Набор и границы значений изменяемых параметров реализации алгоритма:
* размер матрицы от 5000 до 25000 с шагом 500.
* число процессоров 1, 2, 4, [8:144] с шагом 8.
Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:
{|align="center"
|-valign="top"
|[[file:Trlanlib_mpi1.png|thumb|center|700px|Рисунок 3: Производительность в зависимости от числа процессоров и размера матрицы.]]
|[[file:Trlanlib_mpi2.png|thumb|center|700px|Рисунок 4: Время выполнения в зависимости от числа процессоров и размера матрицы.]]
|}
Можно заметить, что данная реализация немасштабируема и не дает ускорения при увеличении числа процессов. Из рисунков 3, 4 хорошо видно, что наилучшая производительность достигается с использованием только одного-двух процессов. Дальнейшее увеличение числа процессов лишь приводит к снижению эффективности. После 8-ми процессов ситуация стабилизируется и дальнейшее распараллеливание мало сказывается на результат
{| class="wikitable mw-collapsible mw-collapsed"
! Использованные модули и компиляторы
|-
|
Модули:
mkl/4.0.2.146
impi/5.0.1-ofa
компилятор mpif90, параметры компиляции:
mpif90 -I/opt/intel/composer_xe_2013.2.146/mkl/include/ia32 -I/opt/intel/composer_xe_2013.2.146/mkl/include -O0 /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_lapack95_lp64.a -L/opt/intel/composer_xe_2013.2.146/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl simple.f90 -o simple lapack_aux.o ../../libtrlan_mpi.a

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

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*библиотека [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-21T22:47:21Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan] на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]
Набор и границы значений изменяемых параметров реализации алгоритма:
* размер матрицы от 5000 до 25000 с шагом 500.
* число процессоров 1, 2, 4, [8:144] с шагом 8.
Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:
{|align="center"
|-valign="top"
|[[file:Trlanlib_mpi1.png|thumb|center|700px|Рисунок 3: Производительность в зависимости от числа процессоров и размера матрицы.]]
|[[file:Trlanlib_mpi2.png|thumb|center|700px|Рисунок 4: Время выполнения в зависимости от числа процессоров и размера матрицы.]]
|}
Можно заметить, что данная реализация немасштабируема и не дает ускорения при увеличении числа процессов. Из рисунков 3, 4 хорошо видно, что наилучшая производительность достигается с использованием только одного-двух процессов. Дальнейшее увеличение числа процессов лишь приводит к снижению эффективности. После 8-ми процессов ситуация стабилизируется и дальнейшее распараллеливание мало сказывается на результат
{| class="wikitable mw-collapsible mw-collapsed"
! Использованные модули и компиляторы
|-
|
Модули:
mkl/11.2.0
impi/4.1.0
компилятор mpif90, параметры компиляции:
mpif90 -I/opt/intel/composer_xe_2013.2.146/mkl/include/ia32 -I/opt/intel/composer_xe_2013.2.146/mkl/include -O0 /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_lapack95_lp64.a -L/opt/intel/composer_xe_2013.2.146/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl simple.f90 -o simple lapack_aux.o ../../libtrlan_mpi.a

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

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*библиотека [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

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

2016-11-21T22:46:43Z

Sfedotov:

== Статья [[Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией]] ==

=== Замечания по тексту ===
* В разделе 1.4 должно быть не перечисление макроопераций, а нужно показать структуру алгоритма на макроуровне. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:31, 1 ноября 2016 (MSK)
Исправил
* В разделе 1.5 кроме блок-схемы нужны также текстовые пояснения. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:31, 1 ноября 2016 (MSK)
Исправил
* В разделе 1.8 должна быть оценка параллельной сложности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:31, 1 ноября 2016 (MSK)
Исправил
* В разделе 1.10 должна быть оценка вычислительной мощности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:31, 1 ноября 2016 (MSK)
Исправил

=== Отсутствующие части ===
* Не заполнен раздел 2.4 описания. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 15:43, 16 ноября 2016 (MSK)
Добавил

* Рисунки в статье не пронумерованы. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:17, 21 ноября 2016 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:17, 21 ноября 2016 (MSK)
* В разделе 2.4 не объяснены множественные особенности графиков. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:17, 21 ноября 2016 (MSK)
* Из графиков в разделе 2.4 следует сделать выводы о масштабируемости реализации. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:17, 21 ноября 2016 (MSK)
Исправлено

Файл:Trlanlib mpi1.png

2016-11-21T22:43:27Z

Sfedotov: Sfedotov загрузил новую версию Файл:Trlanlib mpi1.png

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-21T22:41:37Z

Sfedotov:

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan] на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]
Набор и границы значений изменяемых параметров реализации алгоритма:
* размер матрицы от 5000 до 25000 с шагом 500.
* число процессоров 1, 2, 4, [8:144] с шагом 8.
Для всех запусков требовалось найти 10 минимальных собственных значений диагональной матрицы. Величины производительности и времени выполнения брались из вывода встроенной функции состояния trl_print_info. На следующих рисунках приведены графики производительности и времени выполнения данной реализации в зависимости от изменяемых параметров запуска:
{|align="center"
|-valign="top"
|[[file:Trlanlib_mpi1.png|thumb|center|700px|Рисунок 3: Производительность в зависимости от числа процессоров и размера матрицы.]]
|[[file:Trlanlib_mpi2.png|thumb|center|700px|Рисунок 4: Время выполнения в зависимости от числа процессоров и размера матрицы.]]
|}
Можно заметить, что данная реализация немасштабируема и не дает ускорения при увеличении числа процессов. Из рисунков 3, 4 хорошо видно, что наилучшая производительность достигается с использованием только одного-двух процессов. Дальнейшее увеличение числа процессов лишь приводит к снижению эффективности. После 8-ми процессов ситуация стабилизируется и дальнейшее распараллеливание мало сказывается на результат
{| class="wikitable mw-collapsible mw-collapsed"
! Использованные модули и компиляторы
|-
|
Модули:
mkl/11.2.0
impi/4.1.0
компилятор mpif90, параметры компиляции:
mpif90 -I/opt/intel/composer_xe_2013.2.146/mkl/include/ia32 -I/opt/intel/composer_xe_2013.2.146/mkl/include -O0 /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/composer_xe_2013.2.146/mkl/lib/intel64/libmkl_lapack95_lp64.a -L/opt/intel/composer_xe_2013.2.146/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl simple.f90 -o simple lapack_aux.o ../../libtrlan_mpi.a

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

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*библиотека [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

Файл:Trlanlib mpi2.png

2016-11-21T22:15:30Z

Sfedotov:

Файл:Trlanlib mpi1.png

2016-11-21T22:15:13Z

Sfedotov: Sfedotov загрузил новую версию Файл:Trlanlib mpi1.png

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-20T15:38:43Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
Для исследования алгоритма была протестирована параллельная MPI-реализация библиотеки [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan] на языке Fortran90. Исследование проводилось на суперкомпьютере Ломоносов[http://parallel.ru/cluster/]
Параметры запуска алгоритма:
* размер матрицы от 5000 до 25000 с шагом 500.
* число процессоров от 8 до 128 с шагом 8.
В результате проведённых экспериментов максимальная эффективность реализации составила 9.03%, а минимальная - 0.31%.
На следующих рисунках приведены графики производительности и эффективности данной реализации в зависимости от изменяемых параметров запуска.
{|align="center"
|-valign="top"
|[[file:Trlanlib_mpi1.png|thumb|center|700px|График производительности в зависимости от числа процессоров и размера матрицы.]]
|[[file:Trlanlib_mpi2.png|thumb|center|700px|График эффективности в зависимости от числа процессоров и размера матрицы.]]
|}
{| class="wikitable mw-collapsible mw-collapsed"
! Использованные модули и компиляторы
|-
|
module add slurm

module add impi

module add mkl

Для тестовой программы был использован компилятор mpif90, флаги компиляции:
-I${MKLROOT}/include/intel64/lp64 -I${MKLROOT}/include -O0
линковки:
${MKLROOT}/lib/intel64/libmkl_blas95_lp64.a ${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a -L${MKLROOT}/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl

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

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*библиотека [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

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

2016-11-20T15:38:31Z

Sfedotov:

Файл:Trlanlib mpi1.png

2016-11-20T15:25:34Z

Sfedotov:

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

2016-11-18T16:02:24Z

Sfedotov: /* Статья Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией */

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-11-18T16:02:12Z

Sfedotov:

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Данный алгоритм разбивает исходную задачу поиска собственных значений и собственных векторов на два этапа: построение трехдиагональной матрицы и, непостредственно, сам поиск собственных значений и собственных векторов, аппроксимирующих исходные. Первый этап состоит из последовательного применения метода Ланцоша построения последовательности подпространств Крылова. На втором этапе одним из известных методов строятся оптимальные приближения. Из затратных операций можно выделить умножение матрицы на вектор, вычисление нормы вектора, процесс переортогонализации и поиск всех собственных значений и всех обственных векторов симметричной трехдиагональной матрицы.

== Схема реализации последовательного алгоритма ==
В разделе [[#Математическое описание алгоритма]] уже представлено описание алгоритма.

Пусть задана некая симметрическая матрица <math>A</math>

Шаг 1: Вычислить <math>q_1, \beta_0, q_0</math>

Шаг 2:В цикле последовательно вычислять матрицу <math>T_j,\alpha_j, \beta_j</math>

*Найти собственные значения матрицы <math>T_j</math>
*Если <math>\beta_j=0</math> Выход
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией является итерационным. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Внутри одной итерации ресурсами параллелизма могут быть процесс ортогонализации ''j''-го вектора Ритца, умножение исходной матрицы на вектор, вычисление нормы вектора и другие векторные операции.
Умножение вещественной симметричной матрицы ''A'' размера <math>n \times n</math> на вектор ''b'' длины ''n'' требует последовательного выполнения ''n'' ярусов умножений и сложений.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов с учётом ''k'' итераций оценивается в
<math>\frac{kn^2}{k \log{n}}=\frac{n^2}{\log{n}}</math>
*Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, оценивается в <math>\frac{k\,n^2}{0.5\,n^2+1.5\,n+k\,n+k+1}\approx 2\,k</math> так как в прикладных задачах часто требуется небольшое число собственных значений и справедливо <math>k<<n</math>
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
== Динамические характеристики и эффективность реализации алгоритма ==
== Выводы для классов архитектур ==

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*программный пакет [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

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

2016-10-28T12:54:54Z

Sfedotov:

== Статья [[Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией]] ==

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-28T12:31:37Z

Sfedotov:

{{Assignment}}
{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы.
== Схема реализации последовательного алгоритма ==
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Нетрудно заметить, что процесс ортогонализации ''j''-го вектора Ритца может быть вычислен параллельно, как и умножение исходной матрицы на вектор.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
== Динамические характеристики и эффективность реализации алгоритма ==
== Выводы для классов архитектур ==

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*программный пакет [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-18T17:18:22Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы.
== Схема реализации последовательного алгоритма ==
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализацией требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти.
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Нетрудно заметить, что процесс ортогонализации ''j''-го вектора Ритца может быть вычислен параллельно, как и умножение исходной матрицы на вектор.

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>.
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Алгоритм Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому алгоритм Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого алгоритма Ланцоша имеет смысл использовать блочный (ленточный) вариант. В блочном алгоритме Ланцоша в отличие от простого алгоритма вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого алгоритма в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный алгоритм с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного алгоритма все копии кратного собственного значения будут найдены на одном шаге, в простом же алгоритме Ланцоша они определяются последовательно. Блочный алгоритм, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
== Динамические характеристики и эффективность реализации алгоритма ==
== Выводы для классов архитектур ==

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*программный пакет [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GraphLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-18T16:51:10Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\;\beta_0=0,\;q_0=0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& if \; \beta_j =0\; break\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы
== Схема реализации последовательного алгоритма ==
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализации требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Нетрудно заметить, что процесс переортогонализации ''j''-го приближения может быть вычислен параллельно

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Метод Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому метод Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого метода Ланцоша имеет смысл использовать блочный (ленточный) метод Ланцоша. В блочном методе Ланцоша в отличие от простого метода вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого метода в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный метод с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного метода все копии кратного собственного значения будут найдены на одном шаге, в простом же методе Ланцоша они определяются последовательно. Блочный метод, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
== Динамические характеристики и эффективность реализации алгоритма ==
== Выводы для классов архитектур ==

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*программный пакет [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GrapLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-18T12:26:06Z

Sfedotov:

{{Assignment}}

{{algorithm
| name = Алгоритм Ланцоша с полной переортогонализацией
| serial_complexity = <math>O(n^2)</math>
| pf_height = <math>O(m)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>\frac{n(n + 1)}{2}</math>
| output_data = <math>m(n + 1)</math>
}}

Автор описания: [[Участник:Sfedotov|Федотов С.А.]]

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша''' - итерационный метод вычисления собственных значений и собственных векторов вещественной симметричной матрицы A, основными преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Применяется для сильно больших матриц. Прямой алгоритм, разработанный [https://ru.wikipedia.org/wiki/%D0%9B%D0%B0%D0%BD%D1%86%D0%BE%D1%88,_%D0%9A%D0%BE%D1%80%D0%BD%D0%B5%D0%BB%D0%B8%D0%B9 Корнелием Ланцошом], адаптирует степенной метод поиска наибольших собственных значений и собственных векторов линейной системы размерности ''n'' с ограниченным числом итераций, ''m'', где ''m'' много меньше ''n''. Метод достаточно эффективный с вычислительной точки зрения, однако в первоначальном виде неустойчив. В 1970, Ojalvo и Newman<ref>Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).</ref> показали как добиться устойчивости и применили метод для расчета очень крупных инженерных конструкций, подверженных динамической нагрузке. В своей работе авторы также предложили способ выбора начального приближения и эмпирически оценили число итераций ''m''. Оригинальный метод не учитывает влияние ошибок округления, поэтому на практике для арифметики с плавающей точкой используется вариант метода Ланцоша с полной переортогонализацией. Существует более быстая модификация при почти такой же точности - [[Алгоритм Ланцоша с выборочной ортогонализацией]]. В настоящее время метод широко используется в различных технических областях и имеет множество своих вариантов.

== Математическое описание алгоритма ==
Алгоритм Ланцоша соединяет в себе метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. Входными данными алгоритма служат квадратная матрица ''A=AT''и вектор начального приближения ''b''. Мы пытаемся найти трехдиагональную симметричную матрицу ''Tk=QkTAQk'', собственные значения которой приближают собственные значения матрицы ''A''. Иными словами на ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk,]'' и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''<ref>Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)</ref>:
<math>
\begin{align}
q_1 = & \frac{b}{\|b\|_2},\; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = Aq_j\\
& \alpha_j = q^T_j z\\
& z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}\\
& z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i,\;z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i \\
& \beta_j = \|z\|_2\\
& q_{j+1} = \frac{z}{\beta_j}\\
end \; & for
\end{align}
</math>

Диагональные элементы обозначены как <math>\alpha_j=t_{jj}</math>, а элементы побочной диагонали <math>\beta_j=t_{j-1,j}=t_{j,j-1}</math>. После каждой итерации мы вычисляем <math>\alpha_j, \beta_j</math>, из которых строится матрица
<math>T_j = \begin{pmatrix}
\alpha_1 & \beta_2 & & & & 0 \\
\beta_2 & \alpha_2 & \beta_3 & & & \\
& \beta_3 & \alpha_3 & \ddots & & \\
& & \ddots & \ddots & \beta_{j-1} & \\
& & & \beta_{j-1} & \alpha_{j-1} & \beta_j \\
0 & & & & \beta_j & \alpha_j \\
\end{pmatrix}</math>

Полная переортогонализация соответствует ''повторному'' проведению процесса ортогонализации Грама-Шмидта для того, чтобы почти гарантировать, что <math>z</math> будет ортогонален векторам ''q1,q2,...,qj-1''. Потеря ортогонализации не заставляет алгоритм вести себя непредсказуемым образом. В качестве расплаты мы приобретаем ''повторные копии'' сошедшихся чисел Ритца. То есть для больших ''k'' матрица ''Tk'' может иметь не одно собственное значение, очень близкое к <math>\lambda_i(A)</math>, но много таких собственных значений. Нахождение собственных значений и собственных векторов трехдиагональной матрицы значительно проще чем их вычисление для исходной матрицы. Например при помощи '''MRRR''' алгоритма они могут быть вычислены за ''O''(''k''2) флопов.

== Вычислительное ядро алгоритма ==
За один шаг алгоритма можно выделить следующие вычислительные ядра:
*умножение исходной матрицы на вектор <math>z=Aq_j</math>
*процесс ортогонализации Грама-Шмидта <math>z = z - \sum^{j-1}_{i=1}\left(z^Tq_i\right)q_i</math>

== Макроструктура алгоритма ==
Основные состовляющие алгоритма это умножение матрицы на вектор, процесс переортогонализации и поиск всех собственных значений и всех собственных векторов симметричной трехдиагональной матрицы
== Схема реализации последовательного алгоритма ==
[[file:Lanczos1.jpg|thumb|center|450px|Блок-схема алгоритма]]

== Последовательная сложность алгоритма ==
Для вычисления ''k''-й итерации алгоритма Ланцоша с полной переортогонализации требуется выполнить:
*''n+2k'' умножений вектора на вектор
*''2k'' умножение вектора на константу
*''2k'' сложения векторов
*''O''(''k''2) операций на разложение матрицы ''Tk''
Легко проверить, что ''m'' шагов этого алгоритма потребуют ''O''(''m''2''n'') флопов и ''O''(''mn'') слов памяти
== Информационный граф ==
[[file:Lanczos2.jpg|thumb|center|300px|Вычислительная схема ''j''-й итерации. Синим цветом помечены промежуточные данные]]

== Ресурс параллелизма алгоритма ==
Алгоритм Ланцоша с полной переортогонализацией в целом последовательный. На ''j''-м шаге вычисляются числа <math>\alpha_j, \beta_j</math> и ортонормированный вектор <math>q_{j+1}</math>, которые служат входными данными для ''j+1''-го шага. Нетрудно заметить, что процесс переортогонализации ''j''-го приближения может быть вычислен параллельно

== Входные и выходные данные алгоритма ==
'''Входные данные''': матрица <math>A\in\mathbb{R}^{n\times n}</math>.
Дополнительные ограничения:
* <math>A</math> – симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>.

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

'''Выходные данные''':
*вектор приближенных собственных значений матрицы <math>\Lambda</math> (элементы <math>\lambda_{ii}</math>).
*приближения к собственным векторам ''<math>v_i,i=1,\ldots, m </math>
'''Объём выходных данных''': <math> m(n+1) </math>.

== Свойства алгоритма ==
*Метод Ланцоша не требует хранения исходной матрицы полным массивом: метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать разреженность матрицы или регулярность структуры матрицы. Поэтому метод Ланцоша может быть использован для матриц таких больших порядков, когда методы вращений и отражений не могут быть применены.
*Часто уже при небольшом <math>( j \sim 2 \sqrt{n})</math> числе шагов метода крайние пары Ритца <math>(\lambda_i,v_i|[v_1,v_2,...,v_j]=Q_jV)</math> оказываются хорошими приближениями к крайним собственным парам исходной матрицы (имеются теоретические оценки скорости сходимости, обосновывающие эту возможность), а именно крайние собственные пары часто и нужны.
*Если известно, что искомые (крайние) собственные значения являются кратными, то вместо простого метода Ланцоша имеет смысл использовать блочный (ленточный) метод Ланцоша. В блочном методе Ланцоша в отличие от простого метода вместо одного начального вектора берется ортонормированная система из ''l, l > 1'', векторов. Аналогом векторов <math>z,q_i</math> простого метода в блочном являются прямоугольные матрицы размера ''n * l'', аналогом скаляров <math>\beta_i, \alpha_i</math> - квадратные матрицы порядка ''l'', аналогом трехдиагональной матрицы ''Ti'' - ленточная матрица с полушириной ленты ''l''. Такие ленточные матрицы могут иметь собственные значения кратности вплоть до ''l''. Поэтому блочный метод с ''l'' начальными вектрами используется, когда известно, что среди искомых собственных значений многие имеют кратность ''l''. При использовании блочного метода все копии кратного собственного значения будут найдены на одном шаге, в простом же методе Ланцоша они определяются последовательно. Блочный метод, как и [[Алгоритм Ланцоша для точной арифметики]] может быть дополнен выборочной ортогонализацией.

= Программная реализация алгоритма =
== Особенности реализации последовательного алгоритма ==
== Локальность данных и вычислений ==
== Возможные способы и особенности параллельной реализации алгоритма ==
== Масштабируемость алгоритма и его реализации ==
== Динамические характеристики и эффективность реализации алгоритма ==
== Выводы для классов архитектур ==

== Существующие реализации алгоритма ==
Алгоритм Ланцоша реализован во множестве библиотек и математических программных пакетах.
*библиотека [https://en.wikipedia.org/wiki/ARPACK ARPACK] на языке FORTRAN77
*[http://www.nag.co.uk/numeric/fl/nagdoc_fl23/html/INDEXES/KWIC/lanczos.html NAG Library]
*программный пакет [http://crd-legacy.lbl.gov/~kewu/trlan.html TRLan]
*[https://turi.com/products/create/open_source.html GrapLab] параллельная реализация на C++

= Литература =
<references \>

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-16T20:16:18Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T19:03:17Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T19:00:48Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T18:44:03Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T18:42:19Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T18:16:08Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T18:01:51Z

Sfedotov:

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-15T17:11:20Z

Sfedotov:

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

2016-10-12T07:49:20Z

Sfedotov: /* Перенос */

== Перенос ==

Статья:[[Алгоритм_Ланцоша_с_полной_переортогонализацией]];
Выполняет:[[Участник:Sfedotov| Федотов С.А.]]

Статья: algowiki-project.org/ru/Решение_задачи_Коши_для_системы_ОДУ_методом_Рунге-Кутты_4_порядка;
Выполняют: Алимов Дамир и Антон Лесничий

Статья: Алгоритм кластеризации, основанный на минимальном остовном дереве;
Выполняет: Стандрик Антон Сергеевич (Astandrik)

Статья: [[Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(3)|Нахождение собственных чисел квадратной матрицы методом QR разложения (3)]];
Выполняют: [[Участник:Смирнова_Александра|Смирнова А.С.]], [[Участник:Александра_Киямова|Киямова А.]]

Статья: [[Метод_Якоби_вычисления_сингулярных_чисел_и_векторов]];
Выполняют: [[Участница:Почернина_Елена| Почернина Елена]] и [[Участница:Светлана|Костюкова Светлана]]

Статья: [[Метод_Ньютона_для_решения_систем_нелинейных_уравнений(2)]];
Выполняют: [[Участник:Oleggium|Гирняк О.Р.]], [[Участник:Dimx19|Васильков Д.А.]]

Статья: [[Алгоритм_кластеризации_с_использованием_представлений]];
Выполняют: [[Участник:Cookie|Шатков С.А.]], [[Антохина Ю.А.]]

Статья: [[Быстрое_преобразование_Фурье_(FFT)|Быстрое преобразование Фурье (FFT)]];
Выполняют: [[Участник:Victoria_Stepanischeva|Степанищева Виктория]] и [[Участник:Igor_Munkin|Мунькин Игорь]]

Статья: [[Быстрое_дискретное_преобразование_Фурье|Быстрое_дискретное_преобразование_Фурье]];
Выполняют: [[Участник:ivanov.kir.m | К.М.Иванов]], [[Участник:Pumpkin | А.C.Сапин]]

Участник:Sfedotov/Алгоритм Ланцоша с полной переортогонализацией

2016-10-05T17:50:58Z

Sfedotov: Новая страница: «= Свойства и структура алгоритмов = == Общее описание алгоритма == '''Алгоритм Ланцоша с по…»

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

== Общее описание алгоритма ==
'''Алгоритм Ланцоша с полной переортогонализацией''' является инструментом поиска собственных значений и собственных векторов симметричной матрицы A. Это итерационный многошаговый метод для сильно больших матриц, преимуществами которого являются малые затраты памяти и вычислительных ресурсов. Алгоритм Ланцоша соединяет метод Ланцоша для построения крыловского подпространства с процедурой Рэлея-Ритца. На ''k''-м шаге из ортонормированных векторов Ланцоша строится матрица ''Qk = [q1,q2,...,qk]'', и в качестве приближенных собственных значений матрицы ''A'' принимаются числа Ритца, т. е. собственные значения симметричной трехдиагональной матрицы ''Tk=QkTAQk''. Пусть ''Tk=V<math>\Lambda</math>VT'' есть спектральное разложение матрицы ''Tk'', столбцы матрицы ''QkV'' рассматриваются как приближения к соответсвующим собственным векторам матрицы ''A''.

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

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

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

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

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

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

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

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

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

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

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

== Существующие реализации алгоритма ==

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