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

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

2017-01-14T15:44:20Z

Dlitvinov:

{{assignment}}
Авторы статьи: Литвинов Д.А., 611, Шер А.В., 628. Авторы внесли равнозначный вклад в написание статьи.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

Макроструктура алгоритма Ланцоша с выборочной переортогонализацией представляется в виде совокупности метода Ланцоша с выборочной переортогонализацией построения крыловского подпространства и процедуры Рэлея-Ритца.

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

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

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков.

[[Файл:Algorithm_scheme_(v._3).svg|1350px|thumb|center| Рисунок 1. Детальный граф алгоритма Ланцоша с выборочной переортогонализацией. 
Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой. 
'?' - проверка условия на необходимость переортогонализации. 
На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]]

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

Исследование масштабируемости алгоритма проводилось на суперкомпьютере
[//hpc.cmc.msu.ru/ IBM Blue Gene/P ВМК МГУ.]

Для тестирования была выбрана реализация
[//crd-legacy.lbl.gov/~osni/marques.html#BLZPACK BLZPACK]. BLZPACK написан
на FORTRAN 77 в 2001 году и использует MPI для параллелизма.

Для запуска BLZPACK необходима программа-драйвер. Использовался один из
встроенных драйверов DRVMPI (файл <tt>drv/double/drvmpi.f</tt>) c измененными
параметрами -- они меняются в зависимости от размера матрицы и количества
требуемых пар собственных значений/векторов. Параметры устанавливаются в начале
файла drvmpi.f.

Поскольку в мануале (<tt>doc/manual.ps</tt>) нету полного описания этих параметров,
прокомментируем некоторые из них. Например, для матрицы размера 7500х7500
использовались следующие значения:

<pre>
INTEGER LEIG
PARAMETER (LEIG = 30)

INTEGER LISTOR
PARAMETER (LISTOR = 15000)

INTEGER LN
PARAMETER (LN = 7500)

INTEGER LRSTOR
PARAMETER (LRSTOR = 10000000)

INTEGER NCUV
PARAMETER (NCUV = 1)

INTEGER NNZMAX
PARAMETER (NNZMAX = 28128750)

DOUBLE PRECISION ONE
PARAMETER (ONE = 1.0D0)

DOUBLE PRECISION ZERO
PARAMETER (ZERO = 0.0D0)
</pre>

*<tt>LEIG</tt> -- размер массива для хранения полученных собственных
векторов/значений.
*<tt>LISTOR</tt> -- размер массива <tt>ISTOR</tt>. Формула определения его минимального размера приведена в мануале.
*<tt>LN</tt> -- порядок матрицы <tt>A</tt>, т.е. <tt>n</tt>.
*<tt>LRSTOR</tt> -- размер массива <tt>RSTOR</tt>. Формула определения его минимального размера приведена в мануале.
*<tt>NNZMAX</tt> -- количество строк во входном файле с матрицей А. Драйвер <tt>DRVMPI</tt> прочитывает симметричную верхнетреугольную матрицу из файла в координатном формате (строка, столбец, значение); таким образом, для плотной матрицы <tt>NNZMAX</tt> будет равен <math>(n^2 + n)/2</math>

Сам запуск <tt>DRVMPI</tt> конфигурируется в файле <tt>drvsp1.dat</tt> -- там задаётся количество требуемых
собственных значений/векторов, имя файла со входной матрицей и некоторые другие параметры.

Входные матрицы генерировались следующим простым скриптом на Python 2.7:

<pre>
import random

if __name__ == "__main__":
n = 7500
num_rows = (n*n + n) / 2

with open('matr7500.data', 'w') as f:
for i in xrange(1, n + 1):
j = i # column number
while j <= n:
value = random.uniform(-1, 1)
f.write(" " + str(i) + " " + str(j) + " " + '{:.14e}'.format(value) + "\n")
j += 1

</pre>

Был использован компилятор GNU Fortran 77 (g77), точнее, его обертка
<tt>mpif77</tt> с поддержкой MPI. На Blue Gene/P используется собственная
реализация MPI от IBM, основанная на MPICH. Для сборки BLZPACK на Blue Gene/P
нужно заменить имя компилятора (ключ FC) в файле <tt>sys/MACROS/Linux</tt> на
<tt>mpif77</tt> и запустить <tt>./creator -f -mpi</tt> в корне проекта. После
этого можно собрать драйвер, зайдя в директорию <tt>drv/data</tt> и выполнив
<tt>make drvmpi.x</tt>.

Использовался SMP-режим работы Blue Gene/P, когда на каждом вычислительном
узле выполняется всего один MPI-процесс.

Пример запуска программы:
<pre>
mpisubmit.bg -n 64 -w 00:15:00 -m smp --stdout eigen.out --stderr eigen.err drvmpi.x
</pre>

Тестировались следующие конфигурации:
* Порядок матриц 1000, 2500, 5000, 7500
* Количество процессоров 4, 8, 16, 32, 64, 128

У матрицы размера 1000 производился поиск 100 пар собственных векторов/значений,
у остальных -- 30 пар. Результаты приведены в таблице 1:

{| class="wikitable" style="border: none; background: none;"
|+ align="bottom" style="caption-side: bottom" | Табл. 1. Результаты тестирования. В ячейках таблицы -- время работы алгоритма в секундах.
! colspan="2" rowspan="2" style="border: none; background: none;"|[[File:Pfeil_SO.svg|none|20px]]
! colspan="6"| Число процессоров
|-
! 4 !! 8 !! 16 !! 32 !! 64 !! 128
|-
! rowspan="4"| Размер матрицы
! 1000
| 24 || 17 || 14 || 15 || 16 || 21
|-
! 2500
| 48 || 27 || 18 || 13 || 12 || 13
|-
! 5000
| 207 || 92 || 61 || 37 || 24 || 23
|-
! 7500
| 612 || 298 || 153 || 90 || 50 || 42
|}

Заметно, при всех размерах матриц алгоритм перестает масштабироваться начиная с
некоторого числа процессоров: для матрицы размера 1000 это 32 узла, для матрицы
размера 2500 64 узла и т.д. Это связано с тем, что затраты на коммуникацию между
узлами становится больше, чем выгода от распараллеливания вычислений. При
меньшем же количестве процессоров алгоритм хорошо масштабируется, соответствуя
теоретически вычисленной высоте ярусно-параллельной формы <math>O(kn + k^2)</math>

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

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

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

Неплохой обзор существующих реализаций находится здесь:
[http://slepc.upv.es/documentation/reports/str6.pdf A Survey of Software for Sparse Eigenvalue Problems]

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

* [http://crd-legacy.lbl.gov/~osni/ BLZPACK] FORTRAN 77, MPI. 2000 год.
* [http://www.cas.mcmaster.ca/~qiao/BLKLAN BLKLAN] C, MATLAB. Нет поддержки параллелизма. 2003 год.
* [http://sun.stanford.edu/~rmunk/PROPACK/ PROPACK] FORTRAN 77, MATLAB. OpenMP. Предназначен для нахождения SVD разложения, но алгоритм Ланцоша может быть вызван и напрямую. 2005 год.
* [http://www.netlib.org/lanz/ LANZ] FORTRAN 77. Нет поддержки параллелизма. 1991 год.
* [http://www.netlib.org/svdpack/ SVDPACK] C, FORTRAN 77. Нет поддержки параллелизма. Предназначен для нахождения SVD разложения. 1992 год.

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2017-01-12T20:14:23Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

Макроструктура алгоритма Ланцоша с выборочной переортогонализацией представляется в виде совокупности метода Ланцоша с выборочной переортогонализацией построения крыловского подпространства и процедуры Рэлея-Ритца.

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

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

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков.

[[Файл:Algorithm_scheme_(v._3).svg|1350px|thumb|center| Рисунок 1. Детальный граф алгоритма Ланцоша с выборочной переортогонализацией. 
Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой. 
'?' - проверка условия на необходимость переортогонализации. 
На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]]

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2017-01-11T19:31:37Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

Макроструктура алгоритма Ланцоша с выборочной переортогонализацией представляется в виде совокупности метода Ланцоша с выборочной переортогонализацией построения крыловского подпространства и процедуры Рэлея-Ритца.

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

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

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков.

[[Файл:algorithm_scheme_svg.svg|1350px|thumb|center| Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой. 
'?' - проверка условия на необходимость переортогонализации. 
На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]]

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2016-11-21T21:25:27Z

Dlitvinov:

{{Assignment}}

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

В описанном алгоритме можно выделить следующие макрооперации:

* Умножение матрицы на вектор
* Сложение, перемножение и умножение на число для векторов
* Вычисление скалярного произведения и норм векторов
* Добавление столбца к матрице

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков.

[[Файл:algorithm_scheme_svg.svg|1350px|thumb|center| Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой. 
'?' - проверка условия на необходимость переортогонализации. 
На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]]

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2016-11-21T21:21:56Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

В описанном алгоритме можно выделить следующие макрооперации:

* Умножение матрицы на вектор
* Сложение, перемножение и умножение на число для векторов
* Вычисление скалярного произведения и норм векторов
* Добавление столбца к матрице

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

Ниже представлен граф одной итерации алгоритма с разделением на <math>n</math> потоков.

[[Файл:algorithm_scheme_svg.svg|1350px|thumb|center| Oper - простая операция сложения/вычитания/умножения/деления с числами с плавающей точкой. 
'?' - проверка условия на необходимость переортогонализации. 
На вход итерации подается матрица <math>A</math> и вектор <math>q_j</math>, на выходе получаем новый вектор <math>q_{j+1}</math>.]]

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2016-11-21T16:03:02Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

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

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

В описанном алгоритме можно выделить следующие макрооперации:

* Умножение матрицы на вектор
* Сложение, перемножение и умножение на число для векторов
* Вычисление скалярного произведения и норм векторов
* Добавление столбца к матрице

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2016-11-21T15:56:08Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

{{algorithm
| name = Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией
| serial_complexity = <math>0</math>
| pf_height = <math>0</math>
| pf_width = <math>0</math>
| input_data = <math>0</math>
| output_data = <math>0</math>
}}

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

В описанном алгоритме можно выделить следующие макрооперации:

* Умножение матрицы на вектор
* Сложение, перемножение и умножение на число для векторов
* Вычисление скалярного произведения и норм векторов
* Добавление столбца к матрице

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

* [http://www.caam.rice.edu/software/ARPACK/ ARPACK library]
* [https://www.nada.kth.se/~ruhe/2D5219/work/lanso.m KTH ROYAL INSTITUTE OF TECHNOLOGY]
* [https://www.nag.co.uk/content/nag-library THE NAG LIBRARY]

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

* Дж. Деммель, Вычислительная линейная алгебра. Теория и приложения. Пер. с англ. - М.: Мир, 2001.

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

2016-11-17T21:11:22Z

Dlitvinov:

Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описанные разделы.

{{algorithm
| name = Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией
| serial_complexity = <math>0</math>
| pf_height = <math>0</math>
| pf_width = <math>0</math>
| input_data = <math>0</math>
| output_data = <math>0</math>
}}

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

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

Пусть <math>Q = [Q_k,Q_u]</math> - ортогональная матрица порядка <math>n</math>, причем <math>Q_k</math> и <math>Q_u</math> имеют соответственно размеры <math>n \times k</math> и <math>n \times (n-k)</math>. Столбцы матрицы <math>Q_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>A = V \Lambda V^T</math> есть спектральное разложение матрицы <math>T_k</math>, то столбцы матрицы <math>Q_k V</math> являются приближениями соответствующих собственных векторов и называются векторами Ритца.

Числа и векторы Ритца являются ''оптимальными'' приближениями к собственным значениям и собственным векторам матрицы <math>A</math>: минимум величины <math>\Vert A P_k - P_k D \Vert_2</math> достигается при <math>P_k = Q_k V</math> и <math>D = \Lambda</math> и равен <math>\Vert T_ku \Vert_2</math>.

При применении метода Ланцоша в арифметике с плавающей точкой ортогональность в вычисленных векторах Ланцоша <math>q_k</math> пропадает из-за ошибок округления. Для того чтобы поддерживать ортогональность, проводится повторный процесс ортогонализации Грама-Шмидта. Доказано, что векторы <math>q_k</math> теряют ортогональность весьма систематическим образом, приобретая большие компоненты в отношении уже сошедшихся векторов Ритца. Поэтому можно выборочно ортогонализировать векторы <math>q_k</math> вместо того, чтобы на каждом шаге ортогонализировать вектор ко всем более ранним векторам, как это делается при полной ортогонализации. Таким образом достигается очень высокая степень ортогонализации векторов Ланцоша, и затрачивается очень небольшая дополнительная работа.

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

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

<math>
\begin{align}
q_1 = & b / \Vert b \Vert_2, \; \beta_0 = 0,\; q_0 = 0\\
for \; & j = 1 \; to \; k \\
& z = A\,q_j\\
& \alpha_j = q_j^T z\\
& z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\
& for \; i \le k \\
& \; \; \; if \; \beta_k \|v_i (k)\| \le \sqrt{\epsilon} \Vert T_k \Vert_2\\
& \; \; \; \; \; \; z = z - (y_{i,k} ^ T z) y_{i,k}\\
& end \; for\\
& \beta_j = \Vert z \Vert_2\\
& if \; \beta_j = 0 \; then\\
& \; \; \; Quit\; the\; algorithm\\
& q_{j+1} = z / \beta_j\\
end \; & for
\end{align}
</math>

Где <math>v_i</math> - это столбцы ортогональной матрицы <math>V</math> из спектрального разложения <math>T_k = V \Lambda V^T</math>, а <math>y_{k,i} = Q_k v_i</math> - столбцы матрицы <math>Q_k V</math> - векторы Ритца.

Согласно теореме Пейджа, <math>y_{k,i}^T q_{k+1} = O(\epsilon \Vert A \Vert_2) / \beta_k \Vert v_i(k) \Vert</math>, то есть компонента <math>y_{k,i}^T q_{k+1}</math> вычисленного вектора Ланцоша <math>q_{k+1}</math> в направлении вектора Ритца <math>y_{k,i} = Q_k v_i</math> обратно пропорциональна величине <math>\beta_k \Vert v_i(k) \Vert</math>, являющейся оценкой погрешности для соответствующего числа Ритца. Поэтому, когда число Ритца сходится, а его оценка погрешности приближается к нулю, вектор Ланцоша приобретает большую компоненту в направлении вектора Ритца <math>q_{k,i}</math>, и следует произвести процедуру переортогонализации.

Величина погрешности <math>\beta_k \Vert v_i(k) \Vert</math> считается малой, если она меньше, чем <math>\sqrt{\epsilon} \Vert A \Vert_2</math>, так как в этом случае, согласно теореме Пейджа, компонента <math>\Vert y_{i,k}^T q_{k+1} \Vert = \Vert y_{i,k}^T z \Vert / \Vert z \Vert_2</math> скорее всего превосходит уровень <math>\sqrt{\epsilon}</math> (на практике используется <math>\Vert T_k \Vert_2</math> вместо <math>\Vert A \Vert_2</math>, так как первое известно, а второе не всегда).

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

У описанного алгоритма два ядра, на которые приходится основное время работы алгоритма:

* Вычисление вектора <math>z = A q_j : O(n^2)</math> операций умножения
* Выборочная переортогонализация <math>z = z - (y_{i,k} ^ T z) y_{i,k} : O(n k^2)</math> операций умножения (на практике коэффициент маленький, что обеспечивает бОльшую скорость выполнения алгоритма, чем в случае с полной переортогонализацией)

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

В описанном алгоритме можно выделить следующие макрооперации:

* Умножение матрицы на вектор
* Сложение, перемножение и умножение на число для векторов
* Вычисление скалярного произведения и норм векторов
* Добавление столбца к матрице

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

На вход алгоритму подаются исходная матрица <math>A</math>, вектор начального приближения <math>b</math> и число <math>k</math>.
Далее последовательно выполняются итерационно для всех <math>j: 1 \le j \le k</math> следующие шаги:

# Умножается матрица на вектор для нахождения вектора <math>z</math>
# Вычисляется коэффициент <math>\Alpha_j</math> посредством скалярного перемножения векторов
# Вычисляется новое значение вектора <math>z</math>
# Производятся оценка погрешностей <math>\beta_k \Vert v_i(k) \Vert</math> и выборочная переортогонализация
# Вычисление очередного значения <math>\beta_j</math> - норма вектора <math>z</math>
# Вычисление и добавление к матрице <math>Q</math> нового столбца
# Вычисление собственных значений и собственных векторов матрицы <math>T_j</math>

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

* Вычисление вектора <math>z</math> - <math>n^2</math> операций умножения и <math>n(n-1)</math> операций сложения
* Вычисление скалярного произведения векторов - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление вектора <math>z</math> - <math>2n</math> умножений и <math>2n</math> сложений
* Оценка погрешностей - <math>kn</math> операций умножения
* Выборочная переортогонализация - <math>O(2nj)</math> умножений и <math>O(2j(n-1) + n)</math> сложений
* Вычисление нормы вектора - <math>n</math> операций умножения и <math>(n-1)</math> сложений
* Вычисление нового столбца матрицы <math>Q</math> - <math>n</math> операций умножения

Итого на каждой итерации:

* Умножений: <math>n^2 + n + 2n + kn + O(2nj) + n + n = n^2 + 5n + kn + O(2nj)</math>
* Сложений: <math>n(n-1) + n-1 + 2n + O(2j(n-1) + n) + n-1 = n^2 + 3n + O(2j(n-1) + n) -2</math>

Итого за все время выполнения алгоритма:

* Умножений: <math>kn^2 + 5kn + k^2n + O(k^2n + kn)</math>
* Сложений: <math>kn^2 + 3kn + O((k^2+k)(n-1)+kn) = kn^2 + 3kn + O(k^2(n-1)+2kn-k)</math>

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

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

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

На каждой из итераций алгоритма можно получить выигрыш за счет распараллеливания следующих шагов:

* Вычисление вектора <math>z</math> - умножение матрицы размера <math>n \times n</math> на вектор длины <math>n</math> - <math>n</math> ярусов сложений, <math>n</math> операций умножений в каждом
* Вычисление коэффициента <math>\alpha_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление нового значения вектора <math>z</math>: 2 яруса умножений длины <math>n</math> (умножение вектора на число), 2 яруса вычитаний длины <math>n</math>
* Выборочная переортогонализация: последовательно для всех <math>i \le k</math> выполняется <math>j</math> ярусов сложений, <math>n</math> операций умножений в каждом, <math>n</math> ярусов сложений, <math>j</math> операций умножений в каждом, один ярус вычитаний размера <math>n</math>
* Вычисление <math>\beta_j</math> - скалярное произведение векторов - <math>n</math> ярусов сложений, 1 операция умножения в каждом
* Вычисление <math>q_j</math> - деление вектора на число - один ярус делений размера <math>n</math>

Сложность алгоритма по ширине ЯПФ - <math>O(n)</math>.

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

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

Объем входных данных: <math>n^2 + n + 1</math>

Выходные данные: матрица <math>Q</math> размера <math>n \times k</math>

Объем выходных данных: <math>nk</math>

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

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

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

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

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

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

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

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

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

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

2016-11-16T21:09:32Z

Dlitvinov:

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

2016-11-16T20:24:39Z

Dlitvinov:

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

2016-11-16T13:35:18Z

Dlitvinov:

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

2016-11-14T21:07:21Z

Dlitvinov:

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

2016-11-14T00:45:43Z

Dlitvinov:

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

2016-11-13T23:04:39Z

Dlitvinov: Новая страница: «Автор статьи: Литвинов Д.А., студент, кафедра исследования операций, 611 группа - все описан…»