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

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

2016-10-30T12:30:44Z

Kozlov Vladimir:

{{Assignment}}

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

Автор — Козлов Владимир, 617 гр.
== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм комбинирует '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk + k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-27T10:01:18Z

Kozlov Vladimir: /* Общее описание алгоритма */

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

Автор — Козлов Владимир, 617 гр.
== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм комбинирует '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk + k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-27T01:08:52Z

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

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

* Нет краткой сводки свойств алгоритма в начале статьи
** Добавил.

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

2016-10-27T01:06:47Z

Kozlov Vladimir:

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

Автор — Козлов Владимир, 617 гр.
== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk + k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-16T19:18:05Z

Kozlov Vladimir:

{{Assignment}}
Автор — Козлов Владимир, 617 гр.
== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-16T19:16:50Z

Kozlov Vladimir:

{{Assignment}}
== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-16T19:15:28Z

Kozlov Vladimir: /* Литература */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [https://www.nag.co.uk/content/nag-library The NAG Library] от Numerical Algorithms Group;
* Библиотека [http://www.caam.rice.edu/software/ARPACK/ ARPACK].

== Литература ==
[1] Деммель Дж. Вычислительная линейна алгебра. Теория и приложения. -М.: Мир, 2001. - 430 c.

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

2016-10-16T19:13:28Z

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

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

2016-10-16T19:12:45Z

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

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
* [[https://www.nag.co.uk/content/nag-library|The NAG Library]] от Numerical Algorithms Group;
* Библиотека [[http://www.caam.rice.edu/software/ARPACK/|ARPACK]].

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

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

2016-10-16T19:01:07Z

Kozlov Vladimir: /* Программная реализация алгоритма */

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

2016-10-16T18:59:40Z

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

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

Информационный граф для алгоритма умножения матрицы на вектор можно видеть [[Умножение плотной неособенной матрицы на вектор (последовательный вещественный вариант)#Информационный граф|по этой ссылке]]. Информационный граф для алгоритма вычисления скалярного произведения можно видеть [[Скалярное произведение векторов, вещественная версия, последовательный вариант#Информационный граф|по этой ссылке]]. Информационные графы для алгоритма вычитания двух векторов и для алгоритма деления вектора на число состоят из <math>n</math> изолированных вершин (каждый элемент вектора результата считается независимо). Информационный граф для алгоритма вычисления нормы вектора представляет собой цепь из <math>n</math> вершин.

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

Файл:Lanczoc-full-orth-macro.png

2016-10-16T18:45:03Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-macro.png

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Файл:Lanczoc-full-orth-iter.png

2016-10-16T18:41:34Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-iter.png

Информационный граф одной итерации алгоритма Ланцоша с полной переортогонализацией

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

2016-10-16T18:31:47Z

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

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T18:10:24Z

Kozlov Vladimir: /* Последовательная сложность алгоритма */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k) + O(nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T18:09:42Z

Kozlov Vladimir: /* Ресурс параллелизма алгоритма */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

=== Ресурс параллелизма алгоритма ===
Для реализации <math>j</math>-ой итерации алгоритма нужно выполнить следующие ярусы:
# <math>n</math> ярусов сложений с <math>n</math> операциями умножения в каждом (вычисление <math>z</math>);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\alpha_j</math>);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (первая переортогонализация);
# <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом, <math>1</math> ярус вычитаний размера <math>n</math> (вторая переортогонализация);
# <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом (вычисление <math>\beta_j</math>);
# <math>1</math> ярус делений размера <math>n</math> (вычисление <math>q_{j + 1}</math>).

Итого ширина ЯПФ одной итерации составляет <math>O(n)</math>, высота — <math>4n + 2j + 3</math>. Ширина ЯПФ всего алгоритма составляет <math>O(n)</math>, высота — <math>O(nk) + O(k^2)</math>.

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T17:47:50Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T17:39:42Z

Kozlov Vladimir: /* Входные и выходные данные алгоритма */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

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

Размер входных данных: <math>n^2 + n</math>.

Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.

Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T17:39:27Z

Kozlov Vladimir: /* Входные и выходные данные алгоритма */

== Свойства и структура алгоритма ==
=== Общее описание алгоритма ===
'''Алгоритм Ланцоша''' — это итерационный алгоритм поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Алгоритм применяется, когда матрица <math>A</math> слишком велика, чтобы к ней можно было применять точные прямые методы вычисления собственных значений. Алгоритм '''метод Рэлея — Ритца''' поиска приближённых собственных значений и '''метод Ланцоша''' построения крыловского подпространства.

Метод Рэлея — Ритца является методом поиска <math>k</math> приближённых собственных значений симметричной вещественной матрицы <math>A</math> размера <math>n \times n</math>. Если <math>Q = [Q_k, Q_u]</math> — ортонормированная матрица размера <math>n \times n</math>, <math>Q_k</math> имеет размер <math>n \times k</math>, <math>Q_u</math> имеет размер <math>n \times n - k</math>, то можно записать равенство

<math>T = Q^T A Q = [Q_k, Q_u]^T A [Q_k, Q_u] =
\left[ \begin{array}{cc}
Q_k^T A Q_k & Q_k^T A Q_u\\
Q_u^T A Q_k & Q_u^T A Q_u
\end{array} \right]
= \left[ \begin{array}{cc}
T_{k} & T_{ku}^T\\
T_{ku} & T_{u}
\end{array} \right].</math>

Метод Рэлея — Ритца заключается в том, что собственные значения матрицы <math>T_k = Q_k^T A Q_k</math> объявляются приближёнными собственными значениями матрицы <math>A</math>. Такое приближение является в некотором смысле «наилучшим»: можно показать, что если <math>T_k = V \Lambda V^{-1}</math> — спектральное разложение <math>T_k</math>, то пара <math>(Q_k V, \Lambda)</math> минимизирует функционал <math>L(P_k, D) = \Vert A P_k - P_k D \Vert_2</math>, причём <math>L(Q_k V, \Lambda) = \Vert T_{ku} \Vert_2</math>, то есть <math>A \approx (Q_k V) \Lambda (Q_k V)^{-1}</math>. Из этого также видно, что метод Рэлея — Ритца позволяет получать приближения для собственных векторов матрицы <math>A</math>. Более того, можно показать, что собственные значения <math>T</math> отличаются от некоторых собственных значений <math>A</math> не более чем на <math>\Vert T_{ku} \Vert_2</math>.

Метод Ланцоша — это метод построения матрицы <math>Q</math>, при использовании которого, во-первых, матрица <math>T</math> оказывается симметричной трёхдиагональной, во-вторых, столбцы <math>Q</math> и <math>T</math> вычисляются последовательно. Трёхдиагональность <math>T</math> приводит к следующим явлениям:
# матрица <math>T_k</math> является трёхдиагональной матрицей меньшей размерности, а для трёхдиагональных матриц существуют высокоэффективные методы поиска собственных значений;
# матрица <math>T_{ku}</math> имеет только один ненулевой (возможно) элемент — правый верхний, а значит, для оценки погрешности полученных собственных значений достаточно знать только этот элемент.

В теории в методе Ланцоша для вычисления каждого следующего столбца <math>q_{j + 1}</math> матрицы <math>Q</math> достаточно знать только <math>q_{j - 1}</math> и <math>q_{j}</math> в силу трёхдиагональности матрицы <math>T</math>. На практике из-за ошибок округления, если не предпринимать специальных мер, набор векторов <math>q_{1}, \dots, q_{k}</math> перестаёт быть ортогональным. Для борьбы с этим явлением на каждом шаге метода Ланцоша приходится выполнять так называемую переортогонализацию — повторно запускать процесс ортогонализации Грама — Шмидта.

=== Математическое описание алгоритма ===
Исходные данные: симметричная вещественная матрица <math>A</math> — матрица, для которой будут вычисляться собственные значения, вектор <math>b</math> — начальное приближение для метода Ланцоша.

Выходные данные: трёхдиагональная симметричная вещественная матрица <math>T_k = \left[ \begin{array}{cccc}
\alpha_1 & \beta_1 & & \\
\beta_1 & \ddots & \ddots & \\
& \ddots & \ddots & \beta_{k - 1}\\
& & \beta_{k - 1} & \alpha_k
\end{array} \right]</math>, ортонормированный набор векторов <math>q_1, \dots, q_k</math>: <math>T_k = [q_1, \dots, q_k]^T A [q_1, \dots, q_k]</math>.

Формулы метода:

<math>
\begin{array}{l}
q_1 = b / \Vert b \Vert_2\\
j = \overline{1, k}:\\
\quad z_j = A q_j \\
\quad \alpha_j = q_j^T A q_j = q_j^T z_j \\
\quad z_j' = z_j - \sum_{i=1}^j (z_j^T q_i) q_i \\
\quad z_j'' = z_j' - \sum_{i=1}^j (z_j'^T q_i) q_i\\
\quad \beta_j = \Vert z_j'' \Vert_2\\
\quad q_{j+1} = z_j'' / \Vert z_j'' \Vert_2 = z_j''/\beta_j
\end{array}
</math>

Собственные значения и собственные векторы матрицы <math>T_k</math> ищутся любым прямым алгоритмом.

Вычисление <math>z_j', z_j''</math> — это полная переортогонализация <math>z_j</math> методом Грама — Шмидта. Двойной запуск практически гарантирует, что <math>z_j''</math> будет ортогонален <math>q_1, \dots, q_j</math>. Заметим, что <math>\sum_{i=1}^j (z^T q_i) q_i = Q_j Q_j^T z</math>, где <math>Q_j = [q_1, \dots, q_j]</math>, поэтому переортогонализацию можно переписать в виде <math>z'_j = z_j - Q_j Q_j^T z_j, z_j'' = z_j' - Q_j Q_j^T z_j'</math>.

Если <math>\beta_j = 0</math> для какого-либо <math>j < k</math>, то <math>\Vert T_{ju}\Vert_2 = 0</math>, а значит, собственные значения <math>T_j</math> в точности совпадают с каким-то собственными значениями <math>A</math>. В этом случае дальнейшие вычисления прекращаются и либо используется полученная <math>T_j</math> размерности меньшей, чем <math>k</math>, либо процедура запускается заново с другим начальным вектором <math>b</math>.

=== Вычислительное ядро алгоритма ===
У алгоритма построения матрицы <math>T_k</math> в том виде, как он описан выше, можно выделить два ядра:
# умножение матрицы <math>A</math> на вектор <math>q_j</math>: <math>z_j = A q_j</math> (<math>O(n^2)</math> умножений вещественных чисел);
# переортогонализация: <math>z'_j = z_j - Q_j Q_j^T z_j,\; z_j'' = z_j' - Q_j Q_j^T z_j'</math> (<math>O(nk^2)</math> умножений вещественных чисел).

=== Макроструктура алгоритма ===
Для реализации алгоритма по указанным выше формулам имеет смысл выделить следующие макрооперации:
* арифметические операции над матрицами и векторами (сложение и умножение на число);
* умножение матрицы на вектор и скалярное произведение векторов;
* вычисление нормы вектора;
* добавление к матрице столбца.

=== Схема реализации последовательного алгоритма ===
<math>
\begin{array}{l}
\mathtt{input:}\; A, b, k\\
q_1 = b / \Vert b \Vert_2\\
Q = q_1\\
\mathtt{for}\; j = \overline{1, k}:\\
\quad z = A q_j \\
\quad \alpha_j = q_j^T z \\
\quad z = z - Q \left( Q^T z \right) \\
\quad z = z - Q \left( Q^T z \right)\\
\quad \beta_j = \Vert z \Vert_2\\
\quad \mathtt{if }\; \beta_j = 0\\
\quad \quad \mathtt{exit}\\
\quad \mathtt{end \; if}\\
\quad q_{j+1} = z/\beta_j\\
\quad Q = [Q, q_{j + 1}]\\
\mathtt{end \; for}\\
\mathtt{output:}\; [\alpha_1, \dots, \alpha_k], [\beta_1, \dots, \beta_{k - 1}], [q_1, \dots, q_k].
\end{array}
</math>

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

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

Таким образом, алгоритм имеет последовательную сложность <math>O(n^2k + nk^2)</math>.

=== Информационный граф ===
[[Файл:Lanczoc-full-orth-macro.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 1. Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными. In — входные данные (матрица <math>A</math> и вектор <math>b</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), nrm — вычисление еклидовой нормы вектора, / — деление вектора на число, Iter — итерации алгоритма.]]

[[Файл:Lanczoc-full-orth-iter.png|thumb|1000px|center|alt=Информационный граф алгоритма Ланцоша на уровне отдельных итераций|Рисунок 2. Информационный граф итерации алгоритма Ланцоша на уровне макроопераций. In — входные данные (матрица <math>A</math>), Out — выходные данные (<math>\alpha_j, \beta_j</math> и <math>q_j</math>), dprod — скалярное произведение векторов, mprod — произведение матрицы на вектор, - — разность векторов, nrm — вычисление еклидовой нормы вектора, / — деление вектора на число.]]

=== Ресурс параллелизма алгоритма ===
=== Входные и выходные данные алгоритма ===
Входные данные: вещественная симметричная матрица <math>A</math>, вектор <math>b</math>, число <math>k</math>.
Размер входных данных: <math>n^2 + n</math>.
Выходные данные: набор чисел <math>\alpha_1, \dots, \alpha_k, \beta_1, \dots, \beta_{k - 1}</math>, набор векторов <math>q_1, \dots, q_k</math>.
Размер выходных данных: <math>2k - 1 + nk</math>.

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

== Программная реализация алгоритма ==
=== Масштабируемость алгоритма и его реализации ===
=== Существующие реализации алгоритма ===

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

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

2016-10-16T17:09:05Z

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

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

2016-10-16T17:06:54Z

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

Файл:Lanczoc-full-orth-macro.png

2016-10-16T17:03:58Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-macro.png

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Файл:Lanczoc-full-orth-macro.png

2016-10-16T17:03:18Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-macro.png

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Файл:Lanczoc-full-orth-iter.png

2016-10-16T17:02:00Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-iter.png

Информационный граф одной итерации алгоритма Ланцоша с полной переортогонализацией

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

2016-10-16T17:00:45Z

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

Файл:Lanczoc-full-orth-macro.png

2016-10-16T16:59:28Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-macro.png

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Файл:Lanczoc-full-orth-macro.png

2016-10-16T16:58:30Z

Kozlov Vladimir: Kozlov Vladimir загрузил новую версию Файл:Lanczoc-full-orth-macro.png

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Файл:Lanczoc-full-orth-iter.png

2016-10-16T16:55:45Z

Kozlov Vladimir: Информационный граф одной итерации алгоритма Ланцоша с полной переортогонализацией

Информационный граф одной итерации алгоритма Ланцоша с полной переортогонализацией

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

2016-10-16T15:31:18Z

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

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

2016-10-16T15:30:50Z

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

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

2016-10-16T15:27:43Z

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

Файл:Lanczoc-full-orth-macro.png

2016-10-16T15:27:16Z

Kozlov Vladimir: Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

Информационный граф алгоритма Ланцоша на уровне итераций с входными и выходными данными.

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

2016-10-16T15:17:51Z

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

Участник:Kozlov Vladimir

2016-10-16T15:03:21Z

Kozlov Vladimir: Новая страница: «Козлов Владимир Дмитриевич. Студент ВМК, 2 курс магистратуры, группа 617.»

Козлов Владимир Дмитриевич. Студент ВМК, 2 курс магистратуры, группа 617.

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

2016-10-15T23:23:04Z

Kozlov Vladimir: /* Математическое описание алгоритма */

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

2016-10-15T23:19:36Z

Kozlov Vladimir: /* Последовательная сложность алгоритма */

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

2016-10-15T23:19:12Z

Kozlov Vladimir: /* Последовательная сложность алгоритма */

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

2016-10-15T23:11:50Z

Kozlov Vladimir: /* Вычислительное ядро алгоритма */

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

2016-10-15T23:09:40Z

Kozlov Vladimir: /* Последовательная сложность алгоритма */

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

2016-10-15T22:40:09Z

Kozlov Vladimir: /* Последовательная сложность алгоритма */

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

2016-10-15T22:22:10Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T22:18:31Z

Kozlov Vladimir: /* Математическое описание алгоритма */

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

2016-10-15T22:17:33Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T22:17:01Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T22:16:30Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T22:16:18Z

Kozlov Vladimir: /* Вычислительное ядро алгоритма */

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

2016-10-15T22:15:41Z

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

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

2016-10-15T22:13:21Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T21:54:51Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */

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

2016-10-15T21:20:10Z

Kozlov Vladimir: /* Схема реализации последовательного алгоритма */