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

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

2016-12-05T09:52:47Z

M.grigoriev: /* Замечания */

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

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

=== Замечания по тексту ===
* В разделе 1.4 кроме выделения макроопераций требуется показать и их взаимосвязь. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:44, 25 октября 2016 (MSK)
* Графы на рисунках в разделе 1.7 не являются графами алгоритма (граф алгоритма по определению ациклический, вершины соответствуют операторам, дуги - передаче информации между ними); рисунки должны быть пронумерованы. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:44, 25 октября 2016 (MSK)
** Если показываете на графе (рис.1) входные и выходные данные, лучше их описать не как "Input" и "Output", а подписать конкретные элементы данных. А на рис.2 почему-то входных и выходных данных нет. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 11:27, 1 ноября 2016 (MSK)
* В разделе 1.10 должна быть приведена оценка вычислительной сложности алгоритма. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 17:44, 25 октября 2016 (MSK)
** Ошибся в предыдущем замечании - имелась в виду оценка вычислительной мощности. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 11:24, 1 ноября 2016 (MSK)
* Рис.4 и 5 вставлены в статью крайне неудачно. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:02, 21 ноября 2016 (MSK)
* В разделе 2.4 не приведены все параметры запуска теста - какой компилятор, с какими опциями использовался, какие версии библиотек, на каких узлах проводился запуск и т.д. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:02, 21 ноября 2016 (MSK)
* В разделе 2.4 не объяснены множественные особенности получившихся графиков. [[Участник:ASA|Александр Сергеевич Антонов]] ([[Обсуждение участника:ASA|обсуждение]]) 14:02, 21 ноября 2016 (MSK)

=== Замечания ===
* убрать * в качестве умножения в математических формулах +
* 1.3 хотелось бы подробнее... +
* 1.6 "имеет сложность O(kn^2) и односится к квадратичным алгоритмам" - уточните, а как же "k"?
* 1.7 - хорошие графы, хотелось бы их словесного описания и, может быть, обсуждения (также чуть уменьшить поле отображения рисунков)
* 2.4 Сборка: информация "lanczos-mpi.cpp -o lanczos-mpi" совершенно лишняя +
* Рис.4-5 - странные значения на шкале "Cores" +
* хорошо смотрится самостоятельность реализации на С++

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

2016-12-05T09:52:20Z

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

{{Assignment|ASA}}

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.4, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.4, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительное ядро на каждой итерации алгоритма Ланцоша состоит в вычислении очередного промежуточного вектора <math>z</math> произведением исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации:
:<math>z = Aq_i</math>

При больших <math>k</math>, сопоставимых с <math>n</math>, к вычислительному ядру можно отнести вычисление собственных значений и собственных векторов трёхдиагональной матрицы. Однако на практике <math>k</math> много меньше <math>n</math>.

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

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

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

[[Файл:mgrigoriev3.png|thumb|center|600px|Рисунок 2. Информационный граф линейного оператора <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

[[Файл:mgrigoriev1.png|600px|thumb|center|Рисунок 3. Граф алгоритма с указанием входных и выходных данных. 
<math>\|x\|</math> — вычисление нормы, 
<math>/</math> — операция деления, 
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов, 
<math>*</math> — операция умножения, 
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2), 
▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>, 
<math>\mathsf{A, b, q_0 = 0, \beta_0 = 0}</math> — входные данные, 
<math>\mathsf{\alpha_i, \beta_i}</math> — выходной результат.]]

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

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

* Умножение <math>A^{n * n}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>\log(n)</math> [https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k\log n)</math>-ой сложностью по высоте ЯПФ и <math>O(n^2)</math>-ой по ширине ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

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

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

* Сложность последовательного алгоритма <math>O(kn^2)</math>;
* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k\log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;
** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;
* Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>. При <math>k</math> много меньше <math>n</math> вычислительная мощность ≈ <math> 2k</math>;
* Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);
* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
* Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.

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

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

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

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

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

Для исследования масшабируемости алгоритма была написана реализация[https://github.com/danyanya/ss16-task1] на языке C++ с использованием MPI. Реализация была протестирована на суперкомпьютере Ломоносов[http://parallel.ru/cluster/].

Были исследованы:
* время выполнения программы в зависимости от размера входных данных (матрицы);
* время выполнения в зависимости от размера параллельных нод.

Дополнительно было измерено время работы исходя из оптимизационных флагов компилятора GCC [https://gcc.gnu.org].

Параметры запуска алгоритма:
* размер матрицы от 20000 до 175000 с шагом 2500;
* количество процессоров от 8 до 128 с шагом 8.

Программа запускалась на суперкомьютере "Ломоносов" со следующими характеристиками:
* Компилятор GCC 5.2.1;
* Версия MPI 1.8.4;
* Сборка проводилась командой: <pre>mpic++ -std=c++0x -O2 <CPP_FILE> -lm -static-libstdc++</pre>

Полученная эффективность реализации варьируется в пределах от 2% (на маленьких входных данных) до 45% (на больших входных данных и максимальном числе нодов).

На рисунке 4 представлена эффективность программы при отсутствии флага оптимизации O2.

[[Файл:Grigorev_lanczos_eff_no_02.png|thumb|center|600px|Рисунок 4. Эффективность параллельной реализации алгоритма Ланцоща на MPI.]]

На рисунке 5 представлена эффективность алгоритма с оптимизацией компилятора GCC (-O2).

[[Файл:Grigorev_lanczos_eff_o2.png|thumb|center|600px|Рисунок 5. Эффективность параллельной реализации алгоритма Ланцоща на MPI с использованием оптимизации компилятора -O2.]]

Как видно из графиков выше, средняя эффективность минимальна при малых входных данных (на размере матрицы в 10000).
В среднем данный показатель держится между 9% и 14%.

При этом, оптимизация компилятора MPIC++ позволяет получить выигрыш в эффективности более чем в 3 раза (45,5% против 14,7%). Однако, как можно заметить из графиков, при оптимизации алгоритм ведет себя более нестабильно, скорее всего из-за значительной траты времени на работу с памятью (более планые участки свидетельствует о том, что необходимые данные нашлись в кеше).

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R). Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

Для проверки масшабируемости, алгоритм был реализован на языке C++ с применением MPI [https://github.com/danyanya/ss16-task1].

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

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

2016-11-03T12:32:36Z

M.grigoriev:

{{Assignment}}

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

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

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

[[Файл:mgrigoriev3.png|800px|thumb|center|Рисунок 2. Информационный граф линейного оператора <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

[[Файл:mgrigoriev1.png|800px|thumb|center|Рисунок 3. Граф алгоритма с указанием входных и выходных данных. 
<math>\|x\|</math> — вычисление нормы, 
<math>/</math> — операция деления, 
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов, 
<math>*</math> — операция умножения, 
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2), 
▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>, 
<math>\mathsf{A, b, q_0 = 0, \beta_0 = 0}</math> — входные данные, 
<math>\mathsf{\alpha_i, \beta_i}</math> — выходной результат.]]

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

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

* Умножение <math>A^{n * n}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>\log(n)</math> [https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k * \log n)</math>-ой сложностью по высоте ЯПФ и <math>O(n^2)</math>-ой по ширине ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

* Сложность последовательного алгоритма <math>O(k*n^2)</math>;
* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k * \log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;
** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;
* Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма <math>\frac{k*(2*n^2+8*n-1)+3*n}{n^2+2*k}</math>. При <math>k</math> много меньше <math>n</math> вычислительная мощность ≈ <math> 2*k</math>;
* Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);
* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
* Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

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

2016-11-03T12:29:43Z

M.grigoriev:

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

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

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

[[Файл:mgrigoriev3.png|800px|thumb|center|Рисунок 2. Информационный граф линейного оператора <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

[[Файл:mgrigoriev1.png|800px|thumb|center|Рисунок 3. Граф алгоритма с указанием входных и выходных данных. 
<math>\|x\|</math> — вычисление нормы, 
<math>/</math> — операция деления, 
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов, 
<math>*</math> — операция умножения, 
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2), 
▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>, 
<math>\mathsf{A, b, q_0 = 0, \beta_0 = 0}</math> — входные данные, 
<math>\mathsf{\alpha_i, \beta_i}</math> — выходной результат.]]

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

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

* Умножение <math>A^{n * n}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>\log(n)</math> [https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k * \log n)</math>-ой сложностью по высоте ЯПФ и <math>O(n^2)</math>-ой по ширине ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

* Сложность последовательного алгоритма <math>O(k*n^2)</math>;
* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k * \log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;
** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;
* Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма <math>\frac{k*(2*n^2+8*n-1)+3*n}{n^2+2*k}</math>. При <math>k</math> много меньше <math>n</math> вычислительная мощность ≈ <math> 2*k</math>;
* Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);
* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
* Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

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

2016-11-03T12:26:29Z

M.grigoriev: /* Информационный граф */

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

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

Информационный граф алгоритма можно разбить на две части:
* Граф алгоритма с отображением входных и выходных данных.
* Граф линейного оператора.
[[Файл:mgrigoriev2.png|800px|thumb|left|Рисунок 1. Информационный граф алгоритма.
<math>\|x\|</math> — вычисление нормы,
<math>/</math> — операция деления,
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов,
<math>*</math> — операция умножения,
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),
▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>]]

[[Файл:mgrigoriev3.png|800px|thumb|right|Рисунок 2. Информационный граф линейного оператора <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

[[Файл:mgrigoriev1.png|800px|thumb|left|Рисунок 3. Граф алгоритма с указанием входных и выходных данных.
<math>\|x\|</math> — вычисление нормы,
<math>/</math> — операция деления,
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов,
<math>*</math> — операция умножения,
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),
▽ - условие выхода <math>\mathsf{(\beta_i = 0)}</math>,
<math>\mathsf{A, b, q_0 = 0, \beta_0 = 0}</math> — входные данные,
<math>\mathsf{\alpha_i, \beta_i}</math> — выходной результат.]]

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

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

* Умножение <math>A^{n * n}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>\log(n)</math> [https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k * \log n)</math>-ой сложностью по высоте ЯПФ и <math>O(n^2)</math>-ой по ширине ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

* Сложность последовательного алгоритма <math>O(k*n^2)</math>;
* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k * \log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;
** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;
* Вычислительная мощность алгоритма Ланцоша без переортогонализации из последовательной сложности алгоритма <math>\frac{k*(2*n^2+8*n-1)+3*n}{n^2+2*k}</math>. При <math>k</math> много меньше <math>n</math> вычислительная мощность ≈ <math> 2*k</math>;
* Алгоритм Ланцоша без переортогонализации не является детерминированным из-за того, что возможно выполнение меньшего числа итераций алгоритма, из-за того, что все собственные значения уже вычислены;
* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);
* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;
* Нестабильность алгоритма (эффект ''ложной сходимости'') присуще плавающей арифметике, из-за ошибок округления в которой на очередном этапе может быть построен линейно зависимый от исходных новый вектор, что повлечет невозможность дальнейшего приближения собственных значений.

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

Файл:Mgrigoriev3.png

2016-11-03T12:17:12Z

M.grigoriev:

Файл:Mgrigoriev2.png

2016-11-03T12:16:58Z

M.grigoriev:

Файл:Mgrigoriev1.png

2016-11-03T12:16:25Z

M.grigoriev:

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

2016-10-31T14:49:33Z

M.grigoriev: /* Информационный граф */

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Граф алгоритма с отображением входных и выходных данных.
[[Файл:Algorithm graph.png|800px|thumb|center| Рисунок 1. Граф алгоритма с отображением входных и выходных данных.
<math>\|x\|</math> — вычисление нормы,
<math>/</math> — операция деления,
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов,
<math>\cdot</math> — операция умножения,
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),
▽ - условие выхода,
<math>\mathsf{Input}</math> — входные данные,
<math>\mathsf{Output}</math> — выходные результат.]]

* Подграф линейного оператора.

[[Файл:Linear operator.png|800px|thumb|center| Рисунок 2. Подграф алгоритма, описывающий линейный оператор <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>log(n)</math>[https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k \log n)</math>-ой сложностью по высоте ЯПФ и по ширине ЯПФ <math>O(n^2)</math>.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

* Сложность последовательного алгоритма <math>O(kn^2)</math>;

* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k \log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;

** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

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

2016-10-31T14:39:41Z

M.grigoriev: /* Информационный граф */

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]] (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 2.7), [[Участник:m.grigoriev|М.А.Григорьев]] (1.2, 1.7, 1.8, 1.9, 2.7, 3)

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца[https://en.wikipedia.org/wiki/Rayleigh%E2%80%93Ritz_method]. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Исходя из предложенной последовательной реализации метода, макрооперациями в алгоритме являются:

* Процедура итеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор (состоит из умножения вектора на число и сложения векторов);
** скалярное произведение векторов;
** линейная комбинация векторов (сложение/умножение на вещественные числа);
** вычисление нормы (скалярное произведение векторов и вычисление квадратного корня);
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

При этом, первая макрооперация (построение трехдиагональной матрицы) выполняется строго последовательно. Единственное, что можно распараллелить - это умножение матрицы на вектор. Однако это операция достаточно легковесная, если оперировать, например, разреженной матрицей.

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

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Граф алгоритма с отображением входных и выходных данных.
[[Файл:Algorithm graph.png|800px|thumb|center| Рисунок 1. Граф алгоритма с отображением входных и выходных данных.
<math>\|x\|</math> — вычисление нормы,
<math>/</math> — операция деления,
<math>\mathsf{F}</math> — вычисление линейной комбинации векторов,
<math>\cdot</math> — операция умножения,
<math>\mathsf{Ax}</math> — линейный оператор (Рисунок 2),
<math>\mathsf{Input}</math> — входные данные,
<math>\mathsf{Output}</math> — выходные результат.]]

* Подграф линейного оператора.

[[Файл:Linear operator.png|800px|thumb|center| Рисунок 2. Подграф алгоритма, описывающий линейный оператор <math>\mathsf{Ax}</math>. <math>+</math>, <math>*</math> - операции сложения и умножения.]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложенийж
** При этом сложение элементов вектора длины <math>n</math> можно выполнить за <math>log(n)</math>[https://algowiki-project.org/ru/%D0%9D%D0%B0%D1%85%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D0%B5_%D1%81%D1%83%D0%BC%D0%BC%D1%8B_%D1%8D%D0%BB%D0%B5%D0%BC%D0%B5%D0%BD%D1%82%D0%BE%D0%B2_%D0%BC%D0%B0%D1%81%D1%81%D0%B8%D0%B2%D0%B0_%D1%81%D0%B4%D0%B2%D0%B0%D0%B8%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5%D0%BC];
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает <math>O(k \log n)</math>-ой сложностью по высоте ЯПФ и по ширине ЯПФ <math>O(n^2)</math>.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

* Сложность последовательного алгоритма <math>O(kn^2)</math>;

* По классификации по высоте ЯПФ Алгоритм Ланцоша является алгоритмом со сложностью <math>O(k \log n)</math>, при классификации по ширине ЯПФ его сложность <math>O(n^2)</math>;

** Таким образом, отношение последовательной сложности к параллельной <math>\frac{kn^2}{k \log{n}}</math>;

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

Файл:Algorithm graph.png

2016-10-31T14:12:54Z

M.grigoriev: Граф алгоритма

Граф алгоритма

Файл:Linear operator.png

2016-10-31T14:11:59Z

M.grigoriev: граф линейного оператора без отображения входных и выходных данных. * — операция умножения двух коэффициентов, + — операция сложения.

граф линейного оператора без отображения входных и выходных данных.
* — операция умножения двух коэффициентов, + — операция сложения.

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

2016-10-15T19:57:47Z

M.grigoriev:

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]], [[Участник:m.grigoriev|М.А.Григорьев]]

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Макрооперациями в алгоритме являются:

* Процедура иттеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор;
** скалярное произведение векторов;
** деление (умножение) вектора на вещественное число;
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Общий граф алгоритма, от входа ''IN'' до функции вычислений собственных значений полученной трехдиагональной матрицы (''EIGEN'') (Рис. 1);
[[Файл:Danyanya-lanczos-infopraph-1.png|800px|thumb|center|Общий граф алгоритма. IN - вход, Iter - одна итерация алгоритма, EIGEN - процедура вычисления собственных значений]]

* Подграф на каждой итерации внутри цикла (от ''IN'' до ''NEXT'' в рисунке 2);

[[Файл:Danyanya-lanczos-infopraph-2.png|800px|thumb|center|Подграф алгоритма, описывающий одну итерацию. IN - вход в очедедную итерацию, NEXT - переход к следующей итерации или выход цикла, между ними возможен параллелизм при операциях умножения/сложения, вычитания из q(i) q(i-1)-ого, а также деления нового вектора на норму (полученную в блоке Norm). Красными стрелками показаны переходы в непараллельные секции (переход в NEXT неизбежен, из-за чего помечен зеленым)]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложений.
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и по высоте ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

M.grigoriev

2016-10-15T19:57:17Z

M.grigoriev: Полностью удалено содержимое страницы

Участник:M.grigoriev

2016-10-15T19:57:09Z

M.grigoriev: Новая страница: «Григорьев Михаил Александрович, студент 2 курса магистратуры ВМК, кафедра ИБ, 619 группа»

Григорьев Михаил Александрович, студент 2 курса магистратуры ВМК, кафедра ИБ, 619 группа

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

2016-10-15T19:54:27Z

M.grigoriev: /* Общее описание алгоритма */

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]], [[m.grigoriev|М.А.Григорьев]]

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Макрооперациями в алгоритме являются:

* Процедура иттеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор;
** скалярное произведение векторов;
** деление (умножение) вектора на вещественное число;
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Общий граф алгоритма, от входа ''IN'' до функции вычислений собственных значений полученной трехдиагональной матрицы (''EIGEN'') (Рис. 1);
[[Файл:Danyanya-lanczos-infopraph-1.png|800px|thumb|center|Общий граф алгоритма. IN - вход, Iter - одна итерация алгоритма, EIGEN - процедура вычисления собственных значений]]

* Подграф на каждой итерации внутри цикла (от ''IN'' до ''NEXT'' в рисунке 2);

[[Файл:Danyanya-lanczos-infopraph-2.png|800px|thumb|center|Подграф алгоритма, описывающий одну итерацию. IN - вход в очедедную итерацию, NEXT - переход к следующей итерации или выход цикла, между ними возможен параллелизм при операциях умножения/сложения, вычитания из q(i) q(i-1)-ого, а также деления нового вектора на норму (полученную в блоке Norm). Красными стрелками показаны переходы в непараллельные секции (переход в NEXT неизбежен, из-за чего помечен зеленым)]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложений.
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и по высоте ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

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

2016-10-15T19:45:01Z

M.grigoriev: /* Программная реализация алгоритма */

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]], [[m.grigoriev|М.А.Григорьев]]

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Макрооперациями в алгоритме являются:

* Процедура иттеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор;
** скалярное произведение векторов;
** деление (умножение) вектора на вещественное число;
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Общий граф алгоритма, от входа ''IN'' до функции вычислений собственных значений полученной трехдиагональной матрицы (''EIGEN'') (Рис. 1);
[[Файл:Danyanya-lanczos-infopraph-1.png|800px|thumb|center|Общий граф алгоритма. IN - вход, Iter - одна итерация алгоритма, EIGEN - процедура вычисления собственных значений]]

* Подграф на каждой итерации внутри цикла (от ''IN'' до ''NEXT'' в рисунке 2);

[[Файл:Danyanya-lanczos-infopraph-2.png|800px|thumb|center|Подграф алгоритма, описывающий одну итерацию. IN - вход в очедедную итерацию, NEXT - переход к следующей итерации или выход цикла, между ними возможен параллелизм при операциях умножения/сложения, вычитания из q(i) q(i-1)-ого, а также деления нового вектора на норму (полученную в блоке Norm). Красными стрелками показаны переходы в непараллельные секции (переход в NEXT неизбежен, из-за чего помечен зеленым)]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложений.
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и по высоте ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается встроенной возможности параллелизма - самыми стабильными в этом плане являются ARPACK, а также IETL.

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

M.grigoriev

2016-10-15T19:43:31Z

Григорьев Михаил Александрович, студент 2 курса магистратуры ВМК, кафедра ИБ, 619 группа

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

2016-10-15T19:41:54Z

M.grigoriev: Содержимое страницы заменено на «Участник:Danyanya/Алгоритм_Ланцоша_для_точной_арифметики_(без_переортогон…»

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

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

2016-10-15T19:37:55Z

M.grigoriev:

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

Основные авторы описания: [[Участник:danyanya|Д.Р.Слюсарь]], [[m.grigoriev|М.А.Григорьев]]

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

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

''' Алгоритм Ланцоша поиска собственных значений ''' был опубликован Корнелием Ланцошем в 1950 году <ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref>. Этот итерационный алгоритм применим только к эрмитовым матрицам <math>A</math>.
Метод позволяет за <math>k</math> итераций вычислять <math>k</math>-ое приближение собственных значений и собственных векторов исходной матрицы <math>A</math>.

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

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

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

На вход алгоритма подается эрмитова матрица <math>A = A^\dagger</math> '' (в вещественном случае матрица симметрична)'' ,
{{Шаблон: ASymmetric}}

Алгоритм Ланцоша соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Иными словами, из оргонормированных векторов Ланцоша <ref>https://ru.wikipedia.org/wiki/Биортогонализация_Ланцоша</ref> на каждой итерации строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math>.
В качестве приближенных собственных значений матрицы <math>A</math> берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math>:

<math>
T_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & \dots & 0 \\
\beta_1 & \alpha_2 & \beta_2 & \dots & 0 \\
0 & \beta_2 & \ddots & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & \dots & \dots & \beta_{k-1} & \alpha_k
\end{pmatrix}
</math>

На выходе алгоритма получается собственные векторы и вектор собственных значений матрицы <math>T_k</math>, с помощью которых и будет найдены искомые собственные векторы исходной матрицы <math>A</math>.

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

Вычислительным ядром на каждой итерации является вычисление произведения исходной матрицы <math>A</math> на вектор <math>q_i</math> с предыдущей итерации
:<math>z = Aq_i</math>

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

Макрооперациями в алгоритме являются:

* Процедура иттеративного построения трехдиагональной симметричной матрицы, включающая:
** умножение матрицы на вектор;
** скалярное произведение векторов;
** деление (умножение) вектора на вещественное число;
* Вычисление собственных значений и собственных векторов полученной в ходе работы трехдиагональной симметричной матрицы.

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

'''Исходные данные''': симметричная матрица <math>A</math>, случайный вектор <math>b</math>.

'''Вычисляемые данные''': собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>.

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

<math>
\begin{align}
q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\
for \; & i = 1 \; to \; k \\
& z = Aq_i\\
& \alpha_i = q^T_i z\\
& z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}\\
& \beta_i = \|z\|_2\\
& If \; \beta_i == 0 \; then \\
& \; \; \; \; exit\\
& else \\
& \; \; \; \; q_{i+1} = z / \beta_i \\
end \; & for
\end{align}
</math>

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

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

Последовательная сложность алгоритма рассчитана на основе приведенной выше реализации алгоритма. Исходя из псевдокода последовательно выполняются следующие операции:
* Умножение квадратной матрицы <math>n * n</math> на вектор длины <math>n</math>. Требует <math>n * n</math> умножений и сложений;
* Скалярное произведение векторов длины <math>n</math>. Требует <math>n</math> умножений и сложений;
* Сложение векторов длины <math>n</math>. Требует <math>n</math> сложений;
* Умножение вектора длины <math>n</math> на скаляр. Требует <math>n</math> умножений;
* Нахождение квадратичной нормы вектора длины <math>n</math>. Требует <math>n</math> умножений и сложений + извлечения квадратного корня;
* Нахождение собственных значений и векторов трехдиагональной симметричной матрицы размера <math> k \times k </math>. Наиболее эффективный метод метода «Разделяй-и-властвуй» в среднем требует <math>~O(k^{2.3})</math> операций.

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

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

Информационный граф алгоритма можно разбить на две части:
* Общий граф алгоритма, от входа ''IN'' до функции вычислений собственных значений полученной трехдиагональной матрицы (''EIGEN'') (Рис. 1);
[[Файл:Danyanya-lanczos-infopraph-1.png|800px|thumb|center|Общий граф алгоритма. IN - вход, Iter - одна итерация алгоритма, EIGEN - процедура вычисления собственных значений]]

* Подграф на каждой итерации внутри цикла (от ''IN'' до ''NEXT'' в рисунке 2);

[[Файл:Danyanya-lanczos-infopraph-2.png|800px|thumb|center|Подграф алгоритма, описывающий одну итерацию. IN - вход в очедедную итерацию, NEXT - переход к следующей итерации или выход цикла, между ними возможен параллелизм при операциях умножения/сложения, вычитания из q(i) q(i-1)-ого, а также деления нового вектора на норму (полученную в блоке Norm). Красными стрелками показаны переходы в непараллельные секции (переход в NEXT неизбежен, из-за чего помечен зеленым)]]

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

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

* Умножение <math>A^{(n * n)}</math> на вектор длины <math>n</math> требует <math>n</math> ярусов умножений и сложений.
* Остальные операции в рамках итерации выполняются последовательно (вычисление значений векторов может быть выполнено за 1 ярус):
* Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма.

Исходя из вышеизложенного, алгоритм Ланцоша обладает линейной сложностью по ширине ЯПФ и по высоте ЯПФ.

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

'''Входные данные:''' симметричная вещественная матрица <math>A</math>, случайный вектор <math>b</math>, число итераций <math>k</math>.

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

'''Выходные данные:''' вектор собственных значений <math>\Lambda</math>, матрица собственных векторов <math>E</math>.

'''Объём выходных данных:''' <math>k * (n + 1)</math>.

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

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

* Также алгоритм Ланцоша быстро сходится при вычислении собственных значений матрицы <math>A</math>, находящихся на границе ее спектра (в <math>T_{j}</math> в первую очередь появляются максимальные по модулю собственные значения);

* Из-за использования точной арифметики алгоритм Ланшоца может найти кратные собственные значения, которые на деле оными не являются;

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

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

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

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

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

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

TDB next week

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

В настоящее время существует множество реализаций алгоритма Ланцоша итеративного поиска собственных значений, как входящих в официальные дистрибутивы для вычислений (''ARPACK''), так и неофициальных реализаций, выложенных на Github. Среди них:

1. The IETL Project <ref>http://www.comp-phys.org/software/ietl/</ref>

2. MatLab <ref>https://www.mathworks.com/matlabcentral/newsreader/view_thread/10554</ref>

3. ARPACK <ref>http://www.caam.rice.edu/software/ARPACK/</ref>

4. Julia Math <ref>https://github.com/JuliaMath/IterativeSolvers.jl</ref>

Также существуют официальные реализации на других языках (например R).

Что касается про встроенную возможность параллелизма - самым стабильным в этом плане является ARPACK, а также IETL.

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

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

2016-10-15T15:08:10Z

M.grigoriev: /* Общее описание алгоритма */

{{algorithm
| name = Алгоритм Ланцоша без переортогонализации
| serial_complexity =
| pf_height =
| pf_width =
| input_data =
| output_data =
}}

Основные авторы описания: [[Участник:m.grigoriev|Григорьев М.А.]]

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

=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' был опубликовн физиком и математиком Корнелием Ланцошем<ref>Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).</ref> в 1950 году. Этот метод является частным случаем алгоритма Арнольда в случае, если исходная матрица <math>A</math> - симметрична, и был представлен как итерационный метод вычисления собственных значений симметричной матрицы. Этот метод позволяет за <math>k</math> итераций вычислять <math>k</math> приближений собственных значений и собственных векторов исходной матрицы. Хотя алгоритм и был эффективным в вычислительном смысле, но он на некоторое время был предан забвению из-за численной неустойчивости. Только в 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> модифицировали алгоритм для использования в арифметике с плавающей точкой. Новый метод получил название [[Алгоритм Ланцоша с полной переортогонализацией | алгоритма Ланцоша с полной переортогонализацией]].

Алгоритм соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. Из ортонормированных векторов Ланцоша строится матрица <math>Q_k = [q_1, q_2, \dots, q_k]</math> размерности <math>n \times k</math> и в качестве приближенных собственных значений принимаются числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы <math>T_k = Q^T_k A Q</math> размерности <math>k \times k</math>.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

=== Динамические характеристики и эффективность реализации алгоритма ===

=== Выводы для классов архитектур ===

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

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

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

2016-10-15T14:37:47Z

M.grigoriev: Содержание

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

2016-10-14T11:49:23Z

M.grigoriev: Новая страница: «{{algorithm | name = Алгоритм Ланцоша без переортогонализации | serial_complexity = | pf_height = | pf…»

{{algorithm
| name = Алгоритм Ланцоша без переортогонализации
| serial_complexity =
| pf_height =
| pf_width =
| input_data =
| output_data =
}}