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

Участник:Бобцов Борис/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки

2016-10-26T13:28:29Z

Sergeys79:

{{algorithm
| name = Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки
| serial_complexity = <math> O\left( n \right) </math>, где <math> n </math> - количество подотрезков, на которые разбился исходный отрезок
| pf_height = <math> log(n) </math>
| pf_width = <math> n </math>
| input_data = <math>1</math> функция <math>1</math> переменной, <math>3</math> числа (границы отрезка интегрирования, точность)
| output_data = <math>1</math> число
}}
Автор статьи: [[Участник:Stalles|Бобцов Борис]] (группа 619)

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

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

Пусть требуется вычислить определённый интеграл <math>I = \int_a^b{f(x)dx}</math> от непрерывной функции <math>f(x)</math>. Если может быть найдена первообразная <math>F(x)</math> подынтегральной функции, то по формуле Ньютона-Лейбница:

<math>\int_a^b{f(x)dx} = F(b) - F(a). </math>

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

Приближённые методы вычисления определённого интеграла в большинстве случаев основаны на том, что определённый интеграл численно равен площади криволинейной трапеции, ограниченной кривой <math> y =f(x)</math>, сегментом <math> \left[ a,b \right] </math> оси <math>Ox</math> и вертикальными прямыми <math>x = a</math> и <math>x = b</math>.

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

Разбивая сегмент интегрирования на несколько достаточно малых частей, на каждом из полученных маленьких сегментов мы можем приблизить интегрируемую функцию прямой. Тем самым мы заменим криволинейную трапецию прямолинейной, формула площади которой проста для вычисления.

=== Математическое описание алгоритма ===
----
==== Метод трапеций ====
Будем рассматривать проблему вычисления значения определенного интеграла с некоторой заданной точностью <math> \epsilon </math>:

<math>I(a,b) = \int_a^b{f(x)dx}</math>

Пусть на отрезке <math> [a,b] </math> задана равномерная сетка, содержащая <math> n+1 </math> узел:

<math> x_{i} = a + \frac{(b - a)} {n} i, i = 0,...,n </math>

Тогда, согласно методу трапеций, можно численно найти определенный интеграл от функции <math> f (x) </math> на отрезке <math>[a,b]</math>:

<math> \int_a^b{f(x)dx} = \frac{b - a}{n} \left( \frac{f(x_0)}{2} + \sum^{n-1}_{i=1} {f(x_i)} + \frac{f(x_n)}{2} \right)</math>

Будем полагать значение <math> I </math> найденным с точностью <math> \epsilon </math>, если выполнено условие:

<math> |I_{n_1} - I_{n_2}| \leq \epsilon |I_{n_2}|, n_1 > n_2, </math>
где <math> n_1 </math> и <math> n_2 </math> - количество узлов, образующих две разные сетки.

==== Метод рекурсивного деления ====

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

Метод рекурсивного деления был придуман, чтобы избежать этих недостатков. Суть метода заключается в рекурсивном вызове себя до момента, пока мы не достигнем нужной точности:
<math> I(a,b) = IntTrap(a, b, f(a), f(b)) = IntTrap(a, (b - a) / 2, f(a), f((b - a) / 2)) + IntTrap((b - a) / 2, b, f((b - a) / 2), f(b)) = ... </math>

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

====Метод локального стека ====

Данный алгоритм использует массив данных, доступ к которым осуществляется по принципу стека - первым удаляется последний из добавленных элементов. Теперь у разных процессов есть возможность обрабатывать отрезки независимо друг от друга. В рассматриваемом алгоритме вместо стека можно было использовать и другие способы хранения заданий, например очередь. Стек выбран исключительно из соображений простоты реализации. Алгоритм локального стека можно считать наилучшим из перечисленных последовательных методов для разработки его параллельной версии.

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

[[Файл:123321.png]]

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

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

int sp = 0; //указатель вершины стека

struct
{
double a;
double b;
double fa;
double fb;
double sab;
} stk[1000]; //массив структур, в которых хранятся отложенные значения. 1000 может быть заменена на необходимое число

//макроопределения доступа к данным стека
#define STACK_IS_NOT_FREE ( sp > 0 )

#define PUT_INTO_STACK( a, b, fa, fb, sab )
{
stk[sp].a = a;
stk[sp].b = b;
stk[sp].fa = fa;
stk[sp].fb = fb;
stk[sp].sab = sab;
++sp;
}

#define GET_FROM_STACK( a, b, fa, fb, sab )
{
--sp;
a = stk[sp].a;
b = stk[sp].b;
fa = stk[sp].fa;
fb = stk[sp].fb;
sab = stk[sp].sab;
}

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

Будем считать, что описана функция double f( double x ), возвращающая значение f(x) интегрируемой функции f в точке x.

double f( double x );

int IntTrap(a, b)
{
double I = 0; //значение интеграла
double fa = f( a ); //значение функции <math>y = f(x)</math> в точке <math>a</math>
double fb = f( b ); //значение функции <math>y = f(x)</math> в точке <math>b</math>
double c = 0; //середина отрезка интегрирования
double sab = 0; //
double sac = 0; //
double scb = 0; //
double sabc = 0; //переменные для подсчёта текущих сумм
while( 1 )
{
c = ( a + b ) / 2;
fc = f( c );
sac = ( fa + fc ) * ( c - a ) / 2;
scb = ( fc + fb ) * ( b - c ) / 2;
sabc = sac + scb;
if( abs( sab - sacb ) >= <math>\epsilon</math> * abs ( sacb ) ) //<math>\epsilon</math> - заданная точность
{
PUT_INTO_STACK( a, c, fa, fc, sac );
a = c;
fa = fc;
sab = scb;
}
else
{
I += sacb;
if( STACK_IS_NOT_FREE )
{
break;
}
GET_FROM_STACK( a, b, fa, fb, sab );
}
}
return I;
}

=== Последовательная сложность алгоритма ===
----
Пусть для работы алгоритма было необходимо разбить отрезок <math> \left[a, b\right] </math> на сегменты <math> \left[a=x_0 , x_1\right], \left[a=x_1 , x_2\right],... \left[a=x_{n-1} , x_n = b\right] </math>.
Так как каждый раз, когда нам приходится вычислять площадь трапеции на более коротком сегменте, мы делим отрезок пополам, получаем, что длина каждого отрезка равна <math> \left( x_i - x_{i - 1} \right) = \frac{b - a}{2^k} </math> для некоторого натурального <math> k </math>.

Назовём глубиной метода максимальное <math> k </math> из всех <math> k_i : \left( x_i - x_{i - 1} \right) = \frac{b - a}{2^{k_i}} </math>.

Очевидно, что для каждой функции <math> y </math> и для каждой точности <math> \epsilon </math> глубина алгоритма будет разной. То есть можно считать, что <math> k = k(f(x), \epsilon). </math> Также понятно, что, если <math> k </math>- глубина алгоритма при фиксированных <math> y </math> и <math> \epsilon </math>, то максимально возможное количество полученных отрезков, для которых считалась частичная сумма, и количество точек на отрезке <math> \left[a, b\right] </math> равны соответственно <math> 2^k </math> и <math> 2^k + 1 </math> при <math> k_i = k, i = 1,2,...,n </math>. В этом худшем с точки зрения количества операций случае алгоритм прокрутит цикл while( 1 ) <math> 1 + 2 + 4 + ... + 2^k = (2^{k + 1}-1) </math> раз.

Если не учитывать при расчёте сложности операции работы с памятью, а только сложения/вычитания и умножения*деления, получим, что за каждый выполненный цикл while( 1 ) алгоритм выполняет 7 сложений и 6 умножений. Итого получаем в худшем случае <math> 7*2^k </math> сложений и <math> 6*2^k </math> умножений. Кроме этого, алгоритм посчитает значение <math> f(x) </math> в <math> n+1 </math> точках. При этом понятно, что <math> k \leq log(n). </math>

Получаем, что при разбиении исходного отрезка на <math> n </math> подотрезков получаем глубину алгоритма <math> k \leq log(n) </math> и сложность алгоритма <math> O\left( n + 2^{log n} \right) = O\left( n + n \right) = O\left( n \right) </math>.

При этом из всего вышесказанного сразу видно, что при достаточном количестве вычислительных ресурсов распараллеливания высота ярусно-параллельной формы параллельной программы будет равна глубине <math> k \leq log(n), </math>, а ширина будет не больше, чем количество отрезков <math> n. </math>

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

[[Файл:12321.png]]

На рисунке овалом изображён стек, в котором будем хранить отрезки, сами отрезки представлены точками - вершинами графа. Стрелками - дугами обозначено дробление отрезков на два новых, вдвое меньших предыдущего. В целях наглядности отрезки одинаковой длины были сгруппированы в один слой (ярус) и помещены в один стек. Через <math> k </math> обозначена глубина алгоритма. В общем случае не из каждого отрезка выходят стрелки, это означает, что мы достигли нужной точности на отрезках не минимально возможной при глубине <math> k </math> длины.

По данному рисунку можно проследить ход работы последовательного алгоритма. Он начинает работу в самой левой точке, с ярусом под номером один. Далее он начинает выполнять свою программу, при этом станет двигаться по графу по следующему алгоритму:

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

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

Представленный алгоритм обладает хорошим ресурсом параллелизма. В самом деле, пусть в нашем распоряжении есть доступный всем параллельным процессам список отрезков интегрирования, организованный в виде стека. Назовём его глобальным стеком. Предположим, что перед запуском параллельных процессов в глобальный стек помещается единственная запись - координаты концов отрезка интегрирования - в дальнейшем - "отрезок". Тогда возможно реализовать следующую схему алгоритма, выполняемую каждым из параллельных процессов:

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

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

'''Входные данные''':
* функция <math> y = f(x)</math>;
* точность <math> \epsilon </math>;
* границы отрезка интегрирования <math> a, b.</math>

'''Выходные данные''':
* число <math>I </math>, приближённо равное <math> \int_a^b{f(x)dx}</math> с точностью <math> \epsilon </math>.

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

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

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

Соотношение последовательной и параллельной сложности алгоритма высоко и равно приблизительно <math> \frac{n}{log(n)} </math>.

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

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

На данный момент реализация алгоритма локального стека была указана только в книге Якобовского М.В., Кульковой Е.Ю. "Решение задач на многопроцессорных вычислительных системах с разделяемой памятью."

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

# Шнейдер В. Е. и др. Краткий курс высшей математики. Учеб. пособие для втузов. М., «Высш. школа», 1972. 640 с.
# Якобовский М.В., Кулькова Е.Ю. Решение задач на многопроцессорных вычислительных системах с разделяемой памятью. - М.: СТАНКИН, 2004. – 30 c.

Участница:Sannikovats/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (1)

2016-10-26T13:28:11Z

Sergeys79:

{{algorithm
| name = Вычисление определённого интеграла с использованием адаптивно сгущающейся сетки
| input_data = 1 функция,зависящая от 1 переменной, 3 числа ( точность ε и границы отрезка интегрирования[A,B])
| output_data = 1 число
}}

Основные авторы описания:[[Участница:Sannikovats| Санникова Татьяна Сергеевна]],[[Участник:Plukshin| Лукшин Петр Андреевич]].

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

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

Допустим, нам дана функция <math>f(x)</math> , определенная на отрезке, и возможность получать ее численное значение в любой из точек области определения за фиксированное время. Задача — вычислить определенный интеграл данной функции по заданной области и дать оценку погрешности вычисления. В случае, когда можно вычислить интеграл по формуле Ньютона-Лейбница, зная первообразную искомой функции, проблем не возникает. Однако, если первообразной нет, решение можно получить благодаря численным методам. Так определенному интегралу можно дать простую геометрическую интерпретацию - он равен площади криволинейной трапеции, ограниченной слева и справа вертикальными прямыми x=a, x=b, снизу - прямой y=0, а сверху – кривой y=f(x).

Существуют различные последовательные методы построения адаптивных подвижных сеток, которые позволяют значительно повысить точность вычисления определенного интеграла за счет учета характера изменения подынтегральной функции. Представим общее описание для методов рекурсивного деления и локального стека.

====Метод рекурсивного деления====

Простейший из них - это метод рекурсивного деления, выигрыш от применения которого достигается за счет возможности использования сеток с разным числом узлов на разных участках интервала интегрирования. Интервал интегрирования разбивается на две части и независимо интегрируется функция <math>f(x)</math> на каждой из них при соответствующих шагах интегрирования. Процедуру разбиения отрезков можно рекурсивно повторить до получения отрезков, на которых подынтегральная функция имеет простой вид и может быть аппроксимирована отрезком прямой с заданной точностью. Но при распараллеливании этого алгоритма (т.е. при запуске попарно подпрограмм, обрабатывающих половинки разбиваемого интервала) можно столкнуться со следующей проблемой: в реальной системе число процессоров ограничено и время, необходимое для порождения параллельных процессов существенно не равно нулю. А значит, запуск процессов, число которых значительно превышает число физических процессов системы при решении рассматриваемой задачи, приводит к увеличению времени ее решения по сравнению с последовательной программой. Основная сложность построения параллельного алгоритма на основе рассмотренного последовательного алгоритма, кроется в том, что, хотя несколько ветвей программы могли бы одновременно обрабатывать разные части интервала интегрирования, координаты концов этих отрезков хранятся в программном стеке процесса и фактически недоступны программисту. Если бы они были расположены в некотором явно описанном массиве, можно было бы передать их по частям для обработки разным нитям программы. Поэтому, необходимо рассмотреть метод локального стека, позволяющий распределить вычислительную работу между ограниченным числом процессов.

====Метод локального стека====

В основе этого алгоритма лежит описанный выше метод рекурсивного деления. Только данный метод использует массив данных, доступ к которым осуществляется по принципу стека - первым удаляется последний из добавленных элементов. В рассматриваемом алгоритме вместо стека можно было использовать и другие способы хранения заданий, например очередь. Стек выбран исключительно из соображений простоты реализации. Времена выполнения рассмотренных алгоритмов отличаются примерно на 5% в пользу локального стека, таким образом, алгоритм "локального стека" будем считать наилучшим из имеющихся в нашем распоряжении последовательных методов и на его примере построим дальнейшее описание.

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

''Входные данные:'' заданная функция <math>f(x)</math>, определенная на отрезке <math>[A,B]</math>, заданная точность решения ε. 
''Выходные данные:'' число, представляющее собой приближенное значение интеграла. 
В качестве модельной задачи рассматривается проблема вычисления с точностью ε значения определенного интеграла 
[[Файл:1.jpg]]

Пусть на отрезке <math> [A,B] </math> задана равномерная сетка, содержащая <math> n+1 </math> узел:

<math> x_{i} = A + \frac{(B - A)} {n} i, i = 0,...,n </math>

Тогда, согласно методу трапеций, можно численно найти определенный интеграл от функции <math> f(x) </math> на отрезке <math>[A,B]</math>:
<math> J_n(A,B)=\int_A^B{f(x)dx} = \frac{B - A}{n} \left( \frac{f(x_0)}{2} + \sum^{n-1}_{i=1} {f(x_i)} + \frac{f(x_n)}{2} \right)</math>

[[Файл: 2sann.jpg]] [[Файл:3sann.jpg]]

Будем полагать значение <math> J </math> найденным с точностью <math> \epsilon </math>, если выполнено условие:

<math> |J_{n_1} - J_{n_2}| \leq \epsilon |J_{n_2}|, n_1 > n_2, </math>
где <math> n_1 </math> и <math> n_2 </math> - количество узлов, образующих две разные сетки.

Используемые в алгоритме формулы: 
*Из метода трапеций <math>sAB=\frac{f(A)+f(B)}{2}*(B-A)</math>

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

Основное время работы алгоритма приходится на: 
#нахождение середины отрезка и значения функции в этой точке; 
#вычисление частичных сумм на подобластях методом трапеций; 
#определение выполнения условия точности:
#*если не выполняется, то продолжаем дробить область,"закладывая" в стек границы нового интервала,
#*иначе добавляем к результату вычисленную сумму, завершая цикл.

=== Макроструктура алгоритма ===
В реализации алгоритма используется такая структура данных, как стек, в котором хранятся координаты концов отрезков.
<source lang="cpp">
//Данные, описывающие стек
int sp = 0; // Указатель вершины стека

// Массив структур, в которых хранятся отложенные задания
struct mstk {
double A;//координата левой границы отрезка интегрирования
double B;//координата правой границы отрезка интегрирования
double fA;//значение функции на границе A
double fB;//значение функции на границе B
double s;//частичная сумма
} mstk[N];
</source>
В процедурах работы со стеком доступны следующие операции: 
1)Проверка наличия элементов в стеке:
<source lang="cpp">
#define STACK_IS_NOT_FREE (sp > 0)
</source>
2)Операция "положить в стек":
<source lang="cpp">
#define PUT_INTO_STACK(A, B, fA, fB, s) {
mstk[sp].A = A;
mstk[sp].B = B;
mstk[sp].fA = fA;
mstk[sp].fB = fB;
mstk[sp].s = s;
++sp;
}
</source>
3)Операция "взять из стека верхний элемент":

<source lang="cpp">
#define GET_FROM_STACK(A, B, fA, fB, s) {
--sp;
A = mstk[sp].A;
B = mstk[sp].B;
fA = mstk[sp].fA;
fB = mstk[sp].fB;
s = mstk[sp].s;
}

</source>

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

// Заданная функция f(x)
double f(double x);

Рассмотрим функцию JLS, вычисляющую определенный интеграл от f(x) на [A,B] с заданной точностью eps методом локального стека, которая является вычислительным ядром алгоритма:
<source lang="cpp">
int JLS(double A, double B, double eps) {
double J = 0;//значение интеграла
double fA = f(A);//значение функции f(x)на границе A
double fB = f(B);//значение функции f(x)на границе B
double C;//середина отрезка [A,B]
double fC;//значение функции f(x)в середине отрезка [A,B]
double sAB = (fA + fB)*(B - A)/2;//вычисление площади трапеции с основанием AB
double sAC;//площадь трапеции с основанием AС
double sCB;//площадь трапеции с основанием СВ
double sACB;//сумма sAC и sCB
while(1) {
C = (A + B)/2;//середина отрезка
fC = f(C);
sAC = (fA + fC)*(C - A)/2;
sCB = (fC + fB )*(B - C)/2;
sACB = sAC + sCB;
//проверка выполнения условия на порядок точности:
if(abs(sAB - sACB) >= eps*abs(sACB)) {
PUT_INTO_STACK(A, C, fA, fC, sAC);
A = C;
fA = fC;
sAB = sCB;
}
else {
J += sACB;
if(STACK_IS_NOT_FREE) {
break;
}
GET_FROM_STACK(A, B, fA, fB, sAB);
}
}
return J;
}
</source>

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

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

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

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

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

[[Файл:luksh.jpg]]
На рисунке изображен информационный граф, где k - это глубина алгоритма, а + обозначено суммирование площадей трапеций.

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

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

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

''Входные данные:'' заданная функция <math>f(x)</math>, определенная на отрезке <math>[A,B]</math>, заданная точность решения ε. 
''Выходные данные:'' число, представляющее собой приближенное значение интеграла. 

=== Свойства алгоритма ===
Рассмотренный алгоритм обладает свойством детерминированности, поскольку при определенных входных данных мы будем получать одинаковый результат. При этом он может адаптироваться в зависимости от сложности подынтегральной функции.
Кроме того, алгоритм обладает внутренним параллелизмом, так как является численно устойчивым как в последовательной, так и в параллельной версиях.

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

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

=== Существующие реализации алгоритма ===
В настоящее время проблемой вычисления определенного графа на адаптивной сетке по нашим сведениям занимались только М.В. Якобовский и Е.Ю. Кулькова, рассматривавшие в своей работе "РЕШЕНИЕ ЗАДАЧ НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ С РАЗДЕЛЯЕМОЙ ПАМЯТЬЮ" методы локального и глобального стеков.

==Литература==
Якобовский М.В.
Введение в параллельные методы решения задач
М.В. Якобовский и Е.Ю. Кулькова "РЕШЕНИЕ ЗАДАЧ НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ С РАЗДЕЛЯЕМОЙ ПАМЯТЬЮ"

Участник:Greno4ka/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки

2016-10-26T13:27:57Z

Sergeys79:

Авторы: [[User:Greno4ka|Устинов Г.А.]] [[User:Лиджиев_Данир_Владимирович|Лиджиев Д.В.]]

= ЧАСТЬ. Свойства и структура алгоритмов =

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

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

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

В качестве модельной задачи рассматривается проблема вычисления с точностью ε значение определенного интеграла:

<math>I(a,b) = \int_a^b{f(x)dx}</math>

Пусть на отрезке <math> [a,b] </math> задана равномерная сетка, содержащая <math> n+1 </math> узел:

<math> x_{i} = a + \frac{(b - a)} {n} i, i = 0,...,n </math>

Тогда, согласно методу трапеций, можно численно найти определенный интеграл от функции <math> f (x) </math> на отрезке <math>[a,b]</math>:

<math> \int_a^b{f(x)dx} = \frac{b - a}{n} \left( \frac{f(x_0)}{2} + \sum^{n-1}_{i=1} {f(x_i)} + \frac{f(x_n)}{2} \right)</math>

Будем полагать значение I найденным с точностью ε, если выполнено условие

<math> |I_{n_1} - I_{n_2}| \leq \epsilon |I_{n_2}|, n_1 > n_2, </math>
где <math> n_1 </math> и <math> n_2 </math> - количество узлов, образующих две разные сетки.

Данный метод неэффективен в силу ряда причин:
* в некоторых точках значение функции будет вычисляться несколько раз;
* на всём интервале интегрирования используется равномерная сетка, тогда как число точек, необходимых для достижения заданной точности, зависит от вида функции и может сильно различаться на разных участках сегмента интегрирования.

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

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

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

Метод локального стека
sp=0 // указатель вершины стека
struct {
A, B, fA, fB, s
} stk[1000]

// макроопределения доступа к данным стека

#define STACK_IS_NOT_FREE (sp>0)

#define PUT_INTO_STACK(A, B, fA, fB, s)
{
stk[sp].A=A
stk[sp].B=B
stk[sp].fA=fA
stk[sp].fB=fB
stk[sp].s=s
sp++
}

#define GET_FROM_STACK(A, B, fA, fB, s)
{
sp--
A=stk[sp].A
B=stk[sp].B
fA=stk[sp].fA
fB=stk[sp].fB
s=stk[sp].s
}

== Схема реализации последовательного алгоритма ==
IntTrap (A, B)
{
J=0
fA = f(A)
fB = f(B)
sAB = (fA + fB) * (B - A) / 2
while (1)
{
C = (A + B) / 2
fc = fun(C)
sAC = (fA + fC) * (C - A) / 2
sCB = (fC + fB ) * (B-C)/2
sACB = sAC + sCB
if (|sAB - sACB| ≥ ε|sACB|) {
PUT_INTO_STACK(A, C, fA, fC, sAC)
A = C
fA = fC
sAB = sCB
} else {
J += sACB
if (STACK_IS_NOT_FREE)
break
GET_FROM_STACK(A, B, fA, fB, sAB)
}
}
return J
}

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

.

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

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

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

Входные данные:
* функция <math>y = f(x)</math>;
* точность <math>\epsilon</math>;
* границы отрезка интегрирования <math>a</math>, <math>b</math>.

Выходные данные:
* число <math>I</math>, равное <math>\int_a^b{f(x)dx}</math> с точностью <math>\epsilon</math>.

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

1. Алгоритм адаптируется под конкретные особенности заданной функции

= ЧАСТЬ. Программная реализация алгоритма =

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

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

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

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

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

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

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

Существует реализация в библиотеке [http://alglib.sources.ru/integration/adaptive.php ALGLIB].

= ЧАСТЬ. Литература =

Якобовский М.В., Кулькова Е.Ю. Решение задач на многопроцессорных вычислительных системах с разделяемой памятью. - М.: СТАНКИН, 2004. – 30 c.

Участник:Pandvik/Ортогонализация Грама - Шмидта

2016-10-26T13:27:47Z

Sergeys79:

Авторы описания алгоритма: Павлов Андрей, Филимонов Владимир.
{{algorithm
| name = Ортогонализация Грамма-Шмидта [https://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%BE%D1%86%D0%B5%D1%81%D1%81_%D0%93%D1%80%D0%B0%D0%BC%D0%B0_%E2%80%95_%D0%A8%D0%BC%D0%B8%D0%B4%D1%82%D0%B0]
| serial_complexity = <math>O(N^2M)</math>
| input_data = <math>NM</math>
| output_data = <math>NM</math>
}}

= ЧАСТЬ. Свойства и структура алгоритмов =

== Общее описание алгоритма ==
В конечномерном евклидовом пространстве существует ортонормированный базис.

Для доказательства этого факта требуется находить и строить такие базисы. Построить ортонормированный базис можно, отталкиваясь от некоторого исходного базиса, при помощи алгоритма, который называют процессом ортогонализации Грама — Шмидта<ref>А.А. Самарский, А.В. Гулин, Численные методы</ref>.

Процесс ортогонализации Грама-Шмидта используется для квадратных невырожденных матриц, которые преобразуются, либо уже преобразованы, к верхнему(нижнему) треугольному виду.

Процесс ортогонализации Грама-Шмидта нашёл применение в оптимизации оценивания параметров моделей управления объектом<ref>А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления</ref>, в протоколах безопасности<ref>А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности</ref>, в обработке сигналов<ref>МГУ ВМК, Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов</ref>, в вычислении локальных минимумов целочисленных решёток и многом другом.

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

== Математическое описание алгоритма ==
Исходные данные: квадратная матрица A с линейно независимыми векторами <math>\mathbf{a}_1,...,\mathbf{a}_N | \mathbf{a}_i \in \R^{M}</math>.

Определяется оператор проекции <math>\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b</math>, где <math>\left \langle a,b \right \rangle</math> - скалярное произведение векторов a и b. Данный оператор используется для проецирования вектора a коллинеарно вектору b.

<math>\mathbf{b}_1 = \mathbf{a}_1</math>

<math>\mathbf{b}_2 = \mathbf{a}_2 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_2 </math>

<math>\mathbf{b}_3 = \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_2} \mathbf{a}_3 </math>

<math>\mathbf{b}_4 = \mathbf{a}_4 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_4 - \mathbf{proj}_{\mathbf{b}_2} \mathbf{a}_4 - \mathbf{proj}_{\mathbf{b}_3} \mathbf{a}_4</math>

<math>\vdots</math>

<math>\mathbf{b}_N = \mathbf{a}_N - \sum_{j=1}^{N-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N</math>

На основе каждого вектора <math>\mathbf{b}_j (j = 1 \cdots N)</math> может быть получен нормированный вектор <math>\mathbf{e}_j = \frac{\mathbf{b}_j}{||\mathbf{b}_j||}</math>.

Результаты процесса ортогонализации Грама-Шмидта: <math>\mathbf{b}_1\cdots\mathbf{b}_N</math> - система ортогональных векторов либо система ортонормированных векторов <math>\mathbf{e}_1\cdots\mathbf{e}_N</math>.

Подробное описание можно найти в <ref>А.Н. Канатников, А.П.Крищенко, Линейная Алгебра</ref> <ref>The University of Southern Mississippi, Department of Mathematics, Jim Lambers, Gram-Schmidt Orthogonalization</ref>.

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

Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов <math>\sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N</math> для каждого <math>b_i</math>.

Подробное описание общей сложности алгоритма представлено в разделе [[#Последовательная сложность алгоритма]]

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

Из раздела описания, представленного в разделе [[#Математическое описание алгоритма]] можно сделать вывод, что основную часть метода составляют вычисления проекций векторов и сумм этих проекций.

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

Пусть дано множество векторов <math>{a_1, a_2, \cdots , a_N}</math>. Тогда:

Шаг 1: <math>\mathbf{b}_1 = \mathbf{a}_1</math>

Для вычисления остальных <math>\mathbf{b}_i </math> следует последовательно выполнять Шаг 2 для всех <math> i </math> таких что при <math> 2 \le i \le N </math>.

Шаг 2: <math>\mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{j}} \mathbf{a}_i</math>, где для каждого <math>\mathbf{proj}_b a</math> выполняется <math>\frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle}b</math>.

== Последовательная сложность алгоритма ==
Оценим сложность последовательной реализации алгоритма.

Для вычисления одной операции проекции требуется:

<math>{\color{Blue}\bullet} \quad M \cdot 3</math> - умножений;

<math>{\color{Blue}\bullet} \quad (M-1) \cdot 2</math> - сложений;

<math>{\color{Blue}\bullet} \quad 1 </math> - деление;

Общее число операций, необходимых для вычисления алгоритма следующее:

<math>{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot M \cdot 3</math> - умножений;

<math>{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot (M-1) \cdot 3 </math> - сложений и вычитаний;

<math>{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} </math> - делений.

Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма <math>O(N^2 \cdot M)</math>.

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

Опишем информационный граф алгоритма. В разделе [[#Схематичное описание информационного графа]] будет представлена схема с сокрытием подробностей реализации функции проекции. Такое ограничение позволяет продемонстрировать основные зависимости, не перегружая картинку лишними элементами. В разделе [[#Уточнение информационного графа]] будет раскрыта внутренняя структура операции проекции.

=== Схематичное описание информационного графа ===

Для вычисления <math>b_i</math> требуется найти <math>proj_{b_{j}}a_{i}</math> для всех <math>j \in [1, i]</math>. Следовательно для полного вычисления вектора <math>b_{i}</math> требуется знать все <math>b_{j}</math> с меньшим индексом. Такая зависимость по данным очень не удачна для параллелизма. Однако, если разбить процесс вычисления <math>b_{i}</math> на несколько этапов, соответствующих функциям проекции (<math>proj_{b_{j}}a_{i}</math>), то это позволит производить некоторые предварительные вычисления для <math>b_{j}</math> до момента, когда станут известны все предшествующие ей <math>b_{i}</math>.

На Рис. 1 изображена зависимость каждого из этапов от предыдущих вычислений. Первая строка овалов содержит все проекции, зависящие от <math>b_{1}</math>, вторая строка - все проекции, зависящие от <math>b_{2}</math>, <math>\;\cdots\;</math>, самая верхняя строка содержит проекцию, которая зависит от <math>b_{N-1}</math>.

Рис. 2 показывает зависимость по данным немного в другом формате. Каждая строка представляет собой набор данных, которые требуются для вычисления <math>b_{i}</math> из первого столбца. Второй и последующие столбцы группируют проекции, зависящие от одной из <math>b_{i}</math> и, с помощью стрелки, показывают эту зависимость.

[[Файл:Diagram_Gram-Schmidt_DataFlow.png|center|frame|Рис. 1. Информационный граф алгоритма Грама-Шмидта]]
[[Файл:Diagram_Gram-Schmidt_DataFlow_1.png|center|frame|Рис. 2. Информационный граф алгоритма Грама-Шмидта]]

=== Уточнение информационного графа ===
В этом разделе будет представлено уточнение информационной структуры графа изображенного на Рис. 1 с помощью раскрытия внутренней структуры операции проекции. Полученный граф несёт в себе новую информацию, которую можно использовать для улучшения алгоритма, но становится менее понятным.

На Рис. 3 изображен граф, аналогичный графу Рис. 1, но функция <math>proj</math> изображена более подробно. Оранжевым овалом обозначаются скалярные произведения двух векторов. Голубой треугольник - принимает на вход два числа (являющиеся результатами скалярных произведений) и вычисляет их частное. Зелёным шестиугольником обозначается произведение вектора <math>a_i</math> на число. Фиолетовый ромб складывает два числа, от которых имеется зависимость.

По сравнению с Рис. 1, более детальное описание информационной модели позволяет заметить повторяющиеся построчно скалярные произведения <math>(b_i, b_i)</math>. Это замечание позволяет уменьшить суммарное число операций, но увеличивает связность графа.

[[Файл:Diagram_Gram-Schmidt_DataFlow_proj.png|center|frame|Рис. 3. Уточненный информационный граф алгоритма Грама-Шмидта]]

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

На основе графа зависимости по данным, представленном на Рис. 1 или Рис. 2 раздела [[#Информационный граф]] можно выделить следующий ресурс для параллелизации. После вычисления <math>b_{i}</math> можно начинать параллельно вычислять все <math>proj_{b_i}a_{j}</math>. На Рис. 2 блоки, которые можно вычислять параллельно представлены столбцами.

Кроме того, при очень больших размерностях можно также распараллелить вычисление проекций. Так как в функции проекции участвует пара векторов <math>a</math> и <math>b</math>, а сама операция проекции определяется как <math>\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b</math>, где <math>\left \langle a,b \right \rangle</math> - скалярное произведение векторов <math>a</math> и <math>b</math>, то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.

=== Оценка ускорения параллельной реализации алгоритма ===

Рассмотрим случай, когда присутствует возможность запуска бесконечного числа потоков. Попробуем оценить сложность вычисления ортогональных векторов, при условии, что параллельно будут запускаться все проекции использующие одинаковую <math>b_i</math>. В этом случае, останется [[критический путь]] зависимости по данным: <math>b_1 \leftarrow proj_{b_2}a_{2} \leftarrow b_2 \leftarrow proj_{b_3}a_{3} \leftarrow \ldots \leftarrow proj_{b_{N-1}}a_{N-1} \leftarrow b_N</math>. Следовательно для оценки сложность распараллеленной версии алгоритма достаточно почитать сложность представленной цепочки вычислений.

<math>{\color{Blue}\bullet} \quad M \cdot 3</math> умножения в одной проекции;

<math>{\color{Blue}\bullet} \quad (M-1) \cdot 2</math> сложения;

<math>{\color{Blue}\bullet} \quad 1</math> деление;

В цепочке зависимости данных присутствует <math>N-1</math> проекция, значит всего:

<math>{\color{Blue}\bullet} \quad N \cdot M \cdot 3</math> - умножений;

<math>{\color{Blue}\bullet} \quad N \cdot (M-1) \cdot 2</math> - сложений;

<math>{\color{Blue}\bullet} \quad N \cdot 1</math> - делений.

Получаем итоговую сложность <math>O(N \cdot M)</math>.

Если дополнительно рассматривать возможность параллелизации векторных операций (скалярного произведения и умножения вектора на число), [[критический путь]] можно сократить до следующих операций:

<math>{\color{Blue}\bullet} \quad 3 \cdot N </math> умножений;

<math>{\color{Blue}\bullet} \quad N </math> делений;

<math>{\color{Blue}\bullet} \quad N \cdot logM </math> сложений;

Следовательно оценка сложности параллельного алгоритма станет равна <math>O(N \cdot logM)</math>.

== Входные и выходные данные алгоритма ==
Входные данные: прямоугольная матрица <math>A^{M \times N}</math> (элементы <math>\mathbf{a}_{ij}</math>).

Ограничение: матрица <math>A</math> должна быть невырожденной.

Объём входных данных: <math> N \cdot N</math>.

Выходные данные: матрица <math>B</math> (система ортогональных векторов <math>\mathbf{b}_{1} \cdots \mathbf{b}_{N}</math>).

Объём выходных данных: <math>N \cdot N</math>.

== Свойства алгоритма ==
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов : <math>\frac{N^2 \cdot M}{N \cdot log M} = \frac{N \cdot M}{log M}</math>.

Отношение числа операций к суммарному объему входных и выходных данных - линейное.

Существуют несколько реализаций рассматриваемого алгоритма. Статья <ref>Å. Björck, Numerics of Gram-Schmidt orthogonalization.</ref> описывает три версии данного алгоритма: строковый классический алгоритм Г-Ш (Грама - Шмидта), столбцовый классический алгоритм Г-Ш<ref>Massachusetts Institute of Technology, Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization</ref> и модифицированный алгоритм Г-Ш. Модифицированный алгоритм Грама-Шмидта<ref>С.А. Белов, Н.Ю. Золотых, Численные Методы Линейной Алгебры</ref> является вычислительно более устойчивым <ref>Massachusetts Institute of Technology, Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization</ref>, то есть позволяет уменьшить погрешность вычислений связанную с ошибкой выполнения операций с плавающей точкой.

Процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт <math> 0 </math> (нулевой вектор) на шаге <math> j </math>, если <math>a_{j}</math> является линейной комбинацией векторов <math>a_{1}, \ldots ,a_{j-1}</math>.

= ЧАСТЬ. Программная реализация алгоритма =

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

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

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

== Масштабируемость алгоритма и его реализации ==
TODO: (до 15 ноября)

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

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

== Существующие реализации алгоритма ==
Пример реализации базового алгоритма на языке Си с использованием технологий OpenMP и MPI представлен в статье <ref>Tirana University, Department of Mathematics, Genci Berati, Gram–Schmidt Process in Different Parallel Platforms (Control Flow versus Data Flow)</ref>.

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

<math>{\color{Blue}\bullet} \quad </math> Mathematica - Orthogonalize[{v1,v2,…}] [https://reference.wolfram.com/language/tutorial/VectorOperations.html]

<math>{\color{Blue}\bullet} \quad </math> Maxima - gramschmidt(x) [http://maxima.sourceforge.net/docs/manual/maxima_23.html]

<math>{\color{Blue}\bullet} \quad </math> MatLab - orth(A) [https://www.mathworks.com/help/symbolic/orth.html]

<math>{\color{Blue}\bullet} \quad </math> SciPy - scipy.linalg.orth(A) [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.orth.html] - Ортогонализация с использованием алгоритма SVD.

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

Участник:DVIN/Ортогонализация Грамма-Шмидта

2016-10-26T13:27:36Z

Sergeys79:

{{algorithm
| name = Ортогонализация Грама-Шмидта
| serial_complexity = <math>O(n^3)</math>
| pf_height = <math>O(n)</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>n^2</math>
| output_data = <math>n^2</math>
}}

Основные авторы описания: [[Участник:DVIN|Инжелевская Дарья Валерьевна]]

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===
Создатели алгоритма Эрхард Шмидт и Йорген Педерсен Грам.

'''Эрхард Шмидт''' (14 января 1876 — 6 декабря 1959) — немецкий математик, с 1917 года профессор Берлинского университета. В 1946—58 первый директор Института математики АН ГДР. Основные труды по теории функций, интегральным уравнениям, функциональному анализу. Определил и изучил геометрически гильбертово пространство, используя аналогии с геометрией Евклида. Занимался квадратичными формами. Известен оператор Гильберта-Шмидта и процесс Грама ― Шмидта.

'''Йорген Педерсен Грам''' (27 июня 1850 - 29 апреля 1916) датский математик, родился в герцогстве Шлезвиг, Дания и умер в Копенгагене, Дания. Основные его работы "On series expansions determined by the methods of least squares", "Investigations of the number of primes less than a given number". Математический метод, который носит его имя, процесс Грама-Шмидта, впервые был опубликован в 1883. Теорема Грама и Определитель Грама также назван в его честь.

Процесс Грама-Шмидта — это один из алгоритмов, в которых на основе счётного множества линейно независимых векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math> строится множество ортогональных векторов <math>{\displaystyle \mathbf {b}_{1},\;\ldots ,\;\mathbf {b} _{N}} </math> или ортонормированных векторов <math>{\displaystyle \mathbf {e} _{1},\;\ldots ,\;\mathbf {e}_{N}} </math>, причём так, что каждый вектор <math>{\displaystyle \mathbf {b} _{j}} </math> или <math>{\displaystyle \mathbf {e} _{j}}</math> может быть выражен линейной комбинацией векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\; \mathbf {a} _{j}}</math>.

=== Математическое описание алгоритма ===
==== Классический процесс Грама — Шмидта ====

===== Алгоритм =====
Пусть имеются линейно независимые векторы <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N</math>.

Определим оператор проекции следующим образом: <math>\mathbf{proj}_{\mathbf{b}}\,\mathbf{a} = {\langle \mathbf{a}, \mathbf{b} \rangle \over \langle \mathbf{b}, \mathbf{b}\rangle} \mathbf{b} ,</math>

где <math>\langle \mathbf{a}, \mathbf{b} \rangle</math> — скалярное произведение векторов <math>\mathbf{a}</math> и <math>\mathbf{b}</math>.

Скалярное произведение для двух векторов <math>\mathbf{ a= [a_1, a_2, ...,a_k]}</math> и <math>\mathbf{ b= [b_1, b_2, ..., b_k]}</math> в '''''k'''''-мерном действительном пространстве определяется как:
:<math>\mathbf{a}\cdot\mathbf{b}=\sum_{i=1}^k a_ib_i=a_1b_1+a_2b_2+\cdots+ a_kb_k</math>.

Этот оператор проецирует вектор <math>\mathbf{a}</math> коллинеарно вектору <math>\mathbf{b}</math>.

Ортогональность векторов <math>\mathbf{a}</math> и <math>\mathbf{b}</math> достигается на шаге (2).

Классический процесс Грама — Шмидта выполняется следующим образом:

: <math>
{\begin{array}{lclr}
{\mathbf {b}}_{1}&=&{\mathbf {a}}_{1}&(1)\\
{\mathbf {b}}_{2}&=&{\mathbf {a}}_{2}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{2}&(2)\\
{\mathbf {b}}_{3}&=&{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{2}}}\,{\mathbf {a}}_{3}&(3)\\
{\mathbf {b}}_{4}&=&{\mathbf {a}}_{4}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{4}-{\mathbf {proj}}_{{{\mathbf {b}}_{2}}}\,{\mathbf {a}}_{4}-{\mathbf {proj}}_{{{\mathbf {b}}_{3}}}\,{\mathbf {a}}_{4}&(4)\\
&\vdots &&\\{\mathbf {b}}_{N}&=&{\mathbf {a}}_{N}-\displaystyle \sum _{{j=1}}^{{N-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{N}&(N)
\end{array}}
</math>

На основе каждого вектора <math>\mathbf{b}_j \;(j = 1 \ldots N)</math> может быть получен нормированный вектор: <math>\mathbf{e}_j = {\mathbf{b}_j\over \| \mathbf{b}_j \|}</math> (у нормированного вектора направление будет таким же, как у исходного, а длина — единичной).

Результаты процесса Грама — Шмидта:

<math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math> — система ортогональных векторов либо

<math>\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N</math> — система ортонормированных векторов.

Вычисление <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math> носит название ортогонализации Грама — Шмидта, а <math>\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N</math> — ортонормализации Грама — Шмидта.

===== Доказательство =====
Докажем ортогональность векторов <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math>.

Для этого вычислим скалярное произведение <math>\langle \mathbf{b}_1, \mathbf{b}_2 \rangle</math>, подставив в него формулу (2). Мы получим ноль. Равенство нулю скалярного произведения векторов означает, что эти векторы ортогональны. Затем вычислим скалярное произведение <math>\langle \mathbf{b}_1, \mathbf{b}_3 \rangle</math>, используя результат для <math>\langle \mathbf{b}_1, \mathbf{b}_2 \rangle</math> и формулу (3). Мы снова получим ноль, то есть векторы <math>\mathbf{b}_1</math> и <math>\mathbf{b}_3</math> ортогональны. Общее доказательство выполняется методом математической индукции.

===== Геометрическая интерпретация — вариант 1 =====
[[Файл:Gram-Schmidt-step2.svg|right|frame|Рис. 1. Второй шаг процесса Грама — Шмидта]]
Рассмотрим формулу (2) — второй шаг алгоритма.
Её геометрическое представление изображено на рис. 1:

1 — получение проекции вектора <math>\mathbf{a}_2</math> на <math>\mathbf{b}_1</math>;

2 — вычисление <math>\mathbf{a}_2-\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_2</math>, то есть [[перпендикуляр]]а, которым выполняется проецирование конца <math>\mathbf{a}_2</math> на <math>\mathbf{b}_1</math>.

Этот перпендикуляр — вычисляемый в формуле (2) вектор <math>\mathbf{b}_2</math>;

3 — перемещение полученного на шаге 2 вектора <math>\mathbf{b}_2</math> в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (2).

На рисунке видно, что вектор <math>\mathbf{b}_2</math> ортогонален вектору <math>\mathbf{b}_1</math>, так как <math>\mathbf{b}_2</math> является перпендикуляром, по которому <math>\mathbf{a}_2</math> проецируется на <math>\mathbf{b}_1</math>.

Рассмотрим формулу (3) — третий шаг алгоритма — в следующем варианте:
: <math>
\begin{array}{lcr}
\mathbf{b}_3 = \mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3+\mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3) & & (6) \\
\end{array}
</math>

Её геометрическое представление изображено на рис. 2:
[[Файл:Gram-schmidt-step3base1.png|center|frame|Рис. 2. Третий шаг процесса Грама — Шмидта]]

1 — получение проекции вектора <math>\mathbf{a}_3</math> на <math>\mathbf{b}_1</math>;

2 — получение проекции вектора <math>\mathbf{a}_3</math> на <math>\mathbf{b}_2</math>;

3 — вычисление суммы <math>\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3</math>, то есть проекции вектора <math>\mathbf{a}_3</math> на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>. Эта плоскость закрашена на рисунке серым цветом;

4 — вычисление <math>\mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3)</math>, то есть перпендикуляра, которым выполняется проецирование конца <math>\mathbf{a}_3</math> на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>. Этот перпендикуляр — вычисляемый в формуле (6) вектор <math>\mathbf{b}_3</math>;

5 — перемещение полученного <math>\mathbf{b}_3</math> в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (6).

На рисунке видно, что вектор <math>\mathbf{b}_3</math> ортогонален векторам <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>, так как <math>\mathbf{b}_3</math> является перпендикуляром, по которому <math>\mathbf{a}_3</math> проецируется на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>.

Таким образом, в процессе Грама — Шмидта для вычисления <math>\mathbf{b}_j</math> выполняется проецирование <math>\mathbf{a}_j</math> ортогонально на [[гиперплоскость]]?, формируемую векторами <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}</math>. Вектор <math>\mathbf{b}_j</math> затем вычисляется как разность между <math>\mathbf{a}_j</math> и его проекцией. То есть <math>\mathbf{b}_j</math> — это перпендикуляр? от конца <math>\mathbf{a}_j</math> к гиперплоскости?, формируемой векторами <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}</math>. Поэтому <math>\mathbf{b}_j</math> ортогонален векторам, образующим эту гиперплоскость?.

===== Численная неустойчивость =====
При вычислении на ЭВМ по формулам (1) — (5) векторы <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{N-1}</math> часто не точно ортогональны из-за [[Ошибка квантования|ошибок округления]]. Из-за потери ортогональности в процессе вычислений классический процесс Грама — Шмидта называют [[Численная устойчивость|численно неустойчивым]].

==== Модифицированный процесс Грама — Шмидта ====
Процесс Грама — Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Вместо вычисления <math>\mathbf{b}_j</math> как: <math> \mathbf{b}_j = \mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_j - \ldots - \mathbf{proj}_{\mathbf{b}_{j-1}}\,\mathbf{a}_j \; (7) </math> этот вектор вычисляется следующим образом:
: <math> \begin{array}{lclr}
\mathbf{a}_j^{(1)} & = &\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j & (8) \\
\mathbf{a}_j^{(2)} & = &\mathbf{a}_j^{(1)} - \mathbf{proj}_{\mathbf{b}_2} \, \mathbf{a}_j^{(1)} & (9) \\& \vdots & & \\
\mathbf{a}_j^{(j-2)} & = & \mathbf{a}_j^{(j-3)} - \mathbf{proj}_{\mathbf{b}_{j-2}} \, \mathbf{a}_j^{(j-3)} & (10) \\
\mathbf{b}_j & = & \mathbf{a}_j^{(j-2)} - \mathbf{proj}_{\mathbf{b}_{j-1}} \, \mathbf{a}_j^{(j-2)} & (11)
\end{array} </math>

=== Вычислительное ядро алгоритма ===
'''Вычислительное ядро''' (computational kernel) алгоритма - это часть алгоритма, на которую приходится основное время его работы.

Вычислительное ядро последовательной версии метода ортогонализации Грамма-Шмидта можно составить из множественных (всего их <math>\frac{N(N-1)}{2}</math>) вычислений проекций : <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math>

=== Макроструктура алгоритма ===
Данный алгоритм использует в качестве составных частей другие алгоритмы. В дальнейшем имеет смысл описывать алгоритм не в максимально детализированном виде (т.е. на уровне арифметических операций), а давать только его макроструктуру, то есть описывает структуру и состав макроопераций. Типичной макрооперацией, часто встречающиеся в алгоритме является оператор проекции векторов. Как записано и в описании ядра алгоритма, основную часть метода ортогонализации Грамма-Шмидта составляют множественные (всего их <math>\frac{N(N-1)}{2}</math>) вычисления оператора проекции.

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

Следующий алгоритм реализует нормализацию Грамма-Шмидта. Векторы <math>v_1,...,v_k</math> заменяются набором ортонормированных векторов, которые имеют ту же линейную оболочку.
[[File:Gram-Schmidt ortho.png|none|thumb]]

Вычислительная сложность этого <math>2Nk^2</math> операции с плавающей точкой, где N - размерность векторов.

Последовательность исполнения метода следующая:

1. <math>\mathbf {b}_{1}=\mathbf {a}_{1}</math>

2. Далее для всех векторов <math>\mathbf {b}_{i}</math> для <math>i=2 ... N</math> производим вычисление по следующей формуле:
<math>{\mathbf {b}}_{i}={\mathbf {a}}_{i}-\displaystyle \sum _{{j=1}}^{{i-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}</math>.

В ней на каждом шаге <math>i</math> по очереди вычисляются все <math>{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}</math> для <math>j=1 ... i-1</math>

Пример реализации на Python. Функция работает для произвольного количества векторов любой размерности. При этом если количество векторов <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N</math> больше их размерности или они линейно зависимы то функция возвращает максимально возможное число линейно независимых векторов <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_n</math>, а остальные векторы <math>\mathbf{b}_{n+1},\;\ldots,\;\mathbf{b}_N</math> нулевые.
<source lang="python">
import random

def GramSchmidt(*a):
k=len(a[0])
N=len(a);
b = [[0] * k for i in range(N)]
b[0]=a[0]
for i in range(1,N):
sum=a[i]
for j in range(0,i):
scolar_ab=0
scolar_bb=0
proj=[i for i in range(k)]
for n in range(k):
scolar_ab+=b[j][n]*a[i][n]
scolar_bb+=b[j][n]*b[j][n]
for n in range(k):
proj[n]=(scolar_ab/scolar_bb)*b[j][n]
for n in range(k):
sum[n]-=proj[n]
b[i]=sum
return b;

l1=[random.randrange(0,10) for i in range(3)]
l2=[random.randrange(0,10) for i in range(3)]
l3=[random.randrange(0,10) for i in range(3)]
print(l1,l2,l3)
print(GramSchmidt(l1,l2,l3))
</source>

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

Для построения ортогонального набора векторов в последовательном варианте требуется:

При k=N:

* <math>\frac{N(N-1)}{2}</math> делений,

* <math>\frac{N(N-1)(3N-1)}{2}</math> сложений (вычитаний),

* <math>\frac{3N^2 (N-1)}{2}</math> умножений.

При k≠N

* <math>\frac{N(N-1)}{2}</math> делений,

* <math>\frac{N(N-1)(3k-1)}{2}</math> сложений (вычитаний),

* <math>\frac{3Nk (N-1)}{2}</math> умножений.

При классификации по последовательной сложности, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам с кубической сложностью.

=== Информационный граф ===
Опишем граф алгоритма как аналитически, так и в виде рисунка.

Граф алгоритма состоит из двух групп вершин, расположенных в целочисленных узлах двух областей.

'''Первая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i} = {\langle \mathbf{a_i}, \mathbf{b_j} \rangle \over \langle \mathbf{b_j}, \mathbf{b_j}\rangle} \mathbf{b_j} ,</math> .

Естественно введённые координаты области таковы:

i — меняется в диапазоне от 2 до N, принимая все целочисленные значения;

j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения.

Аргументы операции следующие:

<math>a_i</math>:
элементы входных данных, а именно <math>a_i</math>;

<math>b_j</math>:
результат срабатывания операции, соответствующей вершине из второй группы, с координатой j.

'''Вторая''' группа вершин расположена в двухмерной области, соответствующая ей операция <math>a - b</math>.

Естественно введённые координаты области таковы:

i — меняется в диапазоне от 2 до N, принимая все целочисленные значения;

j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения.

Аргументы операции следующие:

a:

j=1 - входные данные <math>a_j</math>

j>1 - результат срабатывания операции, соответствующей вершине из второй группы, с координатой j-1

b:

результат срабатывания операции, соответствующей вершине из первой группы, с координатой j
[[Файл:picture00001.jpg|мини|upright=2|centre|frame|Рис. 5. Граф алгоритма с отображением входных и выходных данных. Proj - вычисление оператора проекции, F - операция a-b, In - входные данные, Out - результаты.]]

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

В параллельном варианте, в отличие от последовательного, вычисления <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> можно выполнять параллельно. То есть после каждого вычисления <math>b_j</math> можно запускать параллельное вычисление <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> для всех <math>i=j+1..N</math>

Для ортогонализации в параллельном варианте требуется последовательно выполнить следующие ярусы:
* <math>N</math> ярусов в каждом из которых будет выполнятся,
* по <math>N - 1</math> операторов проекций (в каждом из ярусов линейное количество проекций, в зависимости от яруса — от <math>N-1</math> до <math>1</math>),
* по <math>N - 1</math> сложений (в каждом из ярусов линейное количество сложений, в зависимости от яруса — от <math>N-1</math> до <math>1</math>),

При классификации по высоте (количество ярусов в ЯПФ ) ЯПФ, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам со сложностью <math>O(N)</math>.

При классификации по ширине (максимальное количество вершин в ярусе) ЯПФ его сложность будет <math>O(N^2)</math>.

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

'''Входные данные''': множества линейно независимых векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math>, каждый из которых описывается <math>\mathbf{ a_i= [a^i_1, a^i_2, ..., a^i_k]}</math> .

Дополнительные ограничения:
* вектора <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math> линейно независимы, поэтому <math>k \geqslant N</math>.

'''Объём входных данных''': <math>Nk</math>

'''Выходные данные''': множество ортогональных векторов <math>{\displaystyle \mathbf {b} _{1},\;\ldots ,\;\mathbf {b} _{N}} </math>, каждый из которых описывается <math>\mathbf{ b_i= [b^i_1, b^i_2, ..., b^i_k]}</math> .

'''Объём выходных данных''': <math>Nk</math>

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение кубической к квадратичной).

При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, также линейна.

Алгоритм почти полностью детерминирован — единственность результата выполнения гарантирована, однако возможно накопление ошибок округления при использовании классического процесса Грама-Шмидта. При этом модифицированный алгоритм дает более хороший результат.

Алгоритм является численно неустойчивым — ошибки округления могут привести к неортогональности полученных векторов.

Процесс Грама — Шмидта может применяться также к бесконечной последовательности линейно независимых векторов.

Кроме того, процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт <math>\mathbf{0}</math> (нулевой вектор) на шаге <math>j</math>, если <math>\mathbf{a}_j</math> является линейной комбинацией векторов <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_{j-1}</math>. Если это может случиться, то для сохранения ортогональности выходных векторов и для предотвращения деления на ноль при ортонормировании алгоритм должен делать проверку на нулевые векторы и отбрасывать их. Количество векторов, выдаваемых алгоритмом, будет равно размерности подпространства, порождённого векторами (т.е. количеству линейно независимых векторов, которые можно выделить среди исходных векторов).

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

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

Примеры библиотек, содержащих реализации алгоритма ортогонализации Грама-Шмидта:
# [http://numerics.mathdotnet.com/Matrix.html C#]
# [http://www.mathworks.com/help/symbolic/mupad_ref/linalg-orthog.html?searchHighlight=orthogonalization Matlab]
# [https://reference.wolfram.com/language/ref/Orthogonalize.html Wolfram]
# [http://www.nag.co.uk/numeric/Fl/manual/pdf/F05/f05aaf.pdf NAG Fortran]

====Реализация для пакета Mathematica====
Данный скрипт, предназначенный для пакета Mathematica, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в фигурных скобках последней строки. Количество векторов и их координат могут быть произвольными. В данном случае для примера взяты векторы {-2,1,0}, {-2,0,1}, {-0.5,-1,1}.

<source lang='mathematica'>
Projection[v1_, v2_] := (v1.v2*v2)/v2.v2
MultipleProjection[v1_, vecs_] := Plus @@ (Projection[v1, #1] &) /@ vecs
GramSchmidt[mat_] := Fold[Join[#1, {#2 - MultipleProjection[#2, #1]}] &, {}, mat]
GramSchmidt[{{-2, 1, 0}, {-2, 0, 1}, {-0.5, -1, 1}}]
</source>
====Реализация для Maxima====
Данный скрипт, предназначенный для пакета Maxima, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в квадратных скобках. Количество векторов и их координат могут быть произвольными.

В данном случае для примера взяты векторы [-2,1,0], [-2,0,1], [-0.5,-1,1].
<source lang="tsql">
load(eigen);
x: matrix ([-2,1,0],[-2,0,1],[-0.5,-1,1]);
y: gramschmidt(x);
</source>
====Реализация для MATLAB====
В MATLAB существует встроенная функция [Q,R]=qr(A) которая осуществляет ортогонализацияю Грамма-Шмидта столбцов матрицы А. При этом столбцы матрицы предполагаются линейнонезависимыми.

[Q, R] = qr(A) возвращает матрицы Q с ортонормированными столбцами и обратимую верхнетреугольную матрицу R, так что A=Q*R

== Литература ==
* Беклемишев Д. В. Курс аналитической геометрии и линейной алгебры. — М.: Наука.
* Гельфанд И. М. Лекции по линейной алгебре. — 5-е. — М.: Добросвет, МЦНМО, 1998. — 319 с. — ISBN 5-7913-0015-8.
*Ильин В. А., Позняк Э. Г. Линейная алгебра. 6-е изд. — М.: Физматлит, 2010. — 280 с. — ISBN 978-5-9221-0481-4.
*Халмош П. Конечномерные векторные пространства = Finite-Dimensional Vector Spaces. — М.: Физматгиз, 1963. — 263 с.
*Фаддеев Д. К. Лекции по алгебре. — 5-е. — СПб.: Лань, 2007. — 416 с.

== Ссылки ==
* https://ru.wikipedia.org/wiki/Процесс_Грама_―_Шмидта ''Процесс Грама ― Шмидта''

Участник:KibAndrey/Ортогонализация Грама-Шмидта

2016-10-26T13:27:24Z

Sergeys79:

{{algorithm
| name = Ортогонализация Грама-Шмидта
| serial_complexity = <math>O\left(k^2n\right)</math>
| pf_height = <math>O(k\sqrt{n})</math>
| pf_width = <math>O(kn)</math>
| input_data = <math>kn</math>
| output_data = <math>kn</math>
}}

Основные авторы описания: [[Участник:KibAndrey| А.В.Кибанов]], Т.З.Аджиева.

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

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

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

'''Процесс ортогонализации Грама–Шмидта'''<ref>Математическая энциклопедия / Гл. ред.
И. М. Виноградов — М: Советская энциклопедия — Т.4 Ок–Сло, 1984. — 215 с.</ref> — алгоритм построения для данной линейно независимой системы векторов евклидова или эрмитова пространства <math>V</math> ортогональной системы ненулевых векторов, при котором каждый вектор построенной системы линейно выражается через векторы данной системы.

Исторически это первый алгоритм, описывающий процесс ортогониализации. Традиционно его связывают с именами [https://en.wikipedia.org/wiki/J%C3%B8rgen_Pedersen_Gram Йоргена Педерсена Грама] и [https://en.wikipedia.org/wiki/Erhard_Schmidt Эрхарда Шмидта]. Йорган Педерсен Грам<ref name="multiple">Meyer C. D. Matrix analysis and applied linear algebra. — Siam, 2000. — Т.2.</ref> был датским актуарием, который в неявном виде представил суть процесса ортогонализации в 1883 году<ref>[http://gdz.sub.uni-goettingen.de/dms/load/img/?PID=GDZPPN002158604&physid=PHYS_0047 Gram J. P. Ueber die Entwickelung reeller Functionen in Reihen mittelst der Methode der kleinsten Quadrate.] — ''Journal für die reine und angewandte Mathematik.'' — 1883. — Т. 94. — pp. 41–73.</ref>. Нельзя не сказать о том, что метод ортогонализации в несколько ином, модифицированом виде, ранее использовал [https://en.wikipedia.org/wiki/Pierre-Simon_Laplace Пьер-Симон Лаплас] при нахождении массы Сатурна и Юпитера из систем нормальных уравнений и расчете распределения ошибок при этих вычислениях<ref> [https://books.google.ru/books?id=fjAVAAAAQAAJ&printsec=frontcover&hl=ru&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false Marquis de Laplace P. S. Théorie analytique des probabilités.] — V. Courcier, 1820.</ref>. Ясно, что Й. Грам не знал о том, что Лаплас ранее использовал этот метод. Эрхард Шмидт<ref name="multiple">Meyer C. D. Matrix analysis and applied linear algebra. — Siam, 2000. — Т.2.</ref> был учеником великих математиков [https://en.wikipedia.org/wiki/Hermann_Schwarz Германа Шварца] и [https://en.wikipedia.org/wiki/David_Hilbert Давида Гильберта]. Алгоритм ортогонализации был опубликован Э. Шмидтом в 1907 году в его исследовании интегральных уравнений<ref name=Shmidt>Erchard Shmidt. Mathematische Annalen Zur Theorie der linearen und nichtlinearen Integralgleichungen, 1. Teil: Entwicklung willkürlicher
Funktionen nach Systemen vorgeschriebener. — Mathematische Annalen. — 1908. — V. 65. — №. 3. — pp. 370-399.</ref>, которое в свою очередь дало основание для развития теории гильбертовых пространств. Шмидт использовал процесс ортогонализации применительно к геометрии гильбертова пространства решений интегральных уравнений отмечал, что процесс ортогонализации принципиально такой же, какой прежде использовал Й. Грам. В отечественной литературе иногда этот процесс называют ''ортогонализацией Сонина–Шмидта''<ref>Краснов М.Л. Интегральные уравнения. Введение в теорию. — М.: Наука, 1975. — 302 с.</ref><ref>Марченко В.А. Введение в теорию обратных задач спектрального анализа. — Харьков: Акта, 2005. — 143 с.</ref>. Это название связано с тем, что разработанный Сониным<ref>О некоторых неравенствах, относящихся к определенным интегралам. Записки Академии наук по физико-математическому отделению. Серия 8. Т. 6 № 6. — C. 1–53.</ref> метод ортогонализации системы функций для частного случая по существу не отличается от метода ортогонализации системы функций, предложенного ранее Шмидтом<ref name=Shmidt>Erchard Shmidt. Mathematische Annalen Zur Theorie der linearen und nichtlinearen Integralgleichungen, 1. Teil: Entwicklung willkürlicher
Funktionen nach Systemen vorgeschriebener. — Mathematische Annalen, 63 (1907), pp. 433–476.</ref>.

Классический процесс ортогонализации Грама–Шмидта применяется в таких совменных областях, как анализ изображений<ref>Nussbaum S., Menz G. Object-based image analysis and treaty verification: new approaches in remote sensing-applied to nuclear facilities in Iran. — Springer Science & Business Media, 2008.</ref>, цифровая обработка сигналов<ref>Gopi E. S. Algorithm collections for digital signal processing applications using matlab. — Springer Science & Business Media, 2007.</ref>, нейронные сети<ref>[http://thirdworld.nl/application-of-sequential-learning-neural-networks-to-civil-engineering-modeling-problems Rajasekaran S., Suresh D., Pai G. A. V. Application of sequential learning neural networks to civil engineering modeling problems.] — ''Engineering with Computers.'' — 2002. — Т. 18. – №. 2. — pp. 138–147.</ref>.

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

'''Bходные данные:''' <math>k</math> линейно независимых векторов <math>\mathbf{a_1},\mathbf{a_2},\ldots,\mathbf{a_k}</math> длины <math>n</math> <math>\left(\alpha_{ij}\right.</math>, <math>j=1,2, \ldots,n</math>, — координаты вектора <math>\left.\mathbf{a_i}\right)</math> .

'''Вычисляемые данные:'''

<math>k</math> ортогональных векторов <math>\mathbf{b_1},\mathbf{b_2},...,\mathbf{b_k}</math> длины <math>n</math>, причем <math>\mathbf{b_1}=\mathbf{a_1}</math>, либо

<math>k</math> ортонормированных векторов <math>\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}</math> длины <math>n</math>, причем <math>\mathbf{q_i}=\frac{\mathbf{b_i}}{|\mathbf{b_i}|}</math>, <math>i=1,2, \ldots,k</math>.

'''Формулы процесса ортогонализации:'''
:<math>
\begin{align}
\mathbf{b_{1}} & =\mathbf{a_{1}}, \\
\mathbf{b_{2}}& =\mathbf{a_{2}}-\mathbf{b_{1}} proj_{\mathbf{b_1}}\mathbf{a_{2}},\\
& ...\\
\mathbf{b_{i}} & =\mathbf{ a_{i}}-\sum\limits_{j=1}^{i-1}\mathbf{b_{i}} proj_{\mathbf{b_j}}\mathbf{a_{i}}.\\
& ...\\
\mathbf{b_{k}} & =\mathbf{ a_{k}}-\sum\limits_{j=1}^{k-1}\mathbf{b_{k}} proj_{\mathbf{b_j}}\mathbf{a_{k}}.\\
\end{align}
</math>

Здесь <math>proj_{\mathbf{b_j}}\mathbf{a_{i}}</math>, для <math>j=1,...,k-1</math> — это число, равное по величине проекции вектора <math>\mathbf{a_{j}}</math> на ось, проходящую через вектор <math>\mathbf{b_j}</math>.

Формула для ее вычисления, полученная из определения скалярного произведения:
<math>proj_{\mathbf{b_j}}\mathbf{a_{i}}=\dfrac{(ai,bj)}{\mathbf{\left \|b_j\right \|}}</math>, для <math>j=1,...,i-1</math>

Произведение длин <math>\mathbf{|b_1|},\mathbf{|b_2|},\ldots ,\mathbf{|b_k|}</math> равно объему параллелепипеда, построенного на векторах системы <math>\left\{\mathbf{a_i}\right\}</math>, <math>i=1,2,\ldots ,k</math>, как на ребрах.

Явное выражение векторов <math>\mathbf{b_i}</math> для <math>i=1,...,k</math> через <math>\mathbf{a_1},\mathbf{a_2},...,\mathbf{a_k}</math> дает формула

<math> \mathbf{b_i}=
\begin{vmatrix}
(a_1,a_1) & \cdots & (a_1,a_i-1) &a_1 \\
(a_{2},a_1) & \cdots & (a_2,a_i-1) &a_2 \\
\cdots & \cdots & \cdots& \cdots \\
(a_{i},a_1) & \cdots & (a_i,a_i-1) &a_i \\
\end{vmatrix}
</math>.
(В правой части этого равенства определитель следует формально разложить по последнему столбцу).

Иногда полученные векторы нормируются сразу после их нахождения и находится система ортонормированных векторов <math>\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}</math>. В этом случае знаменатель в формуле для вычисления проекции <math>|\mathbf{q_i}|=1</math> для <math>i=1,...,k-1</math>, что существенно упрощает вычисления алгоритма. В настоящей статье применяя процесс ортогонализации к системе линейно независимых векторов, мы так и будем поступать.

Если в описанном алгоритме в формулах для вычисления векторов <math>\mathbf{b_i}\,</math>, для <math>i=1,...,k</math> процесса ортогонализации заменить <math>proj_{\mathbf{b_j}}\mathbf{a_{i}}</math>, для <math>j=1,...,i-1</math> на <math>proj_{\mathbf{b_j}}\mathbf{b_{i}}</math>, для <math>j=1,...,i-1</math>,алгоритм, полученный таким образом называется ''модифицированный алгоритм Грама–Шмидта''.

Чаще всего, процесс Грама–Шмидта приводит к значительному нарушению ортогональности. Что касается модифицированного процесса Грама–Шмидта, он проявляет себя достаточно неусточиво, хоть и чуть менее, чем классический процесс.

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

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

*<math>\dfrac{k(k+1)}{2}</math> [[Скалярное произведение векторов, вещественная версия, последовательный вариант| скалярных произведений]] векторов,
включая нахождения скалярных произведений одинаковых векторов дальнейшего нахождения их длин;

*<math>\dfrac{(k-1)k}{2}</math> умножений векторов на число после вычисления нового базисного вектора на итерации.

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

Как уже отмечено в [[#Вычислительное ядро алгоритма| описании вычислительного ядра алгоритма]],
основную часть метода составляют вычисления [[Скалярное_произведение_векторов,_вещественная_версия,_последовательно-параллельный_вариант| скалярных произведений]],
как для вычислении длины очередного вектора так и для расчета проекций прочих векторов на новый вектор,
а также умножения векторов на числа, используемые при корректировке
значений остальных векторов после вычисления очередного вектора, ортогонального предыдущим.

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

'''Входные данные:''' <math>k</math> линейно независимых векторов <math>\mathbf{a_1},\mathbf{a_2},\ldots,\mathbf{a_k}</math> длины <math>n</math> <math>\left(\alpha_{ij}\right.</math>, где <math>j=1,2, \ldots,n</math> — координаты вектора <math>\mathbf{a_i}</math>, для <math>\left.i=1,2, \ldots,k\right)</math>.

'''Вычисляемые данные:'''

<math>k</math> ортогональных векторов <math>\mathbf{b_1},\mathbf{b_2},...,\mathbf{b_k}</math> длины <math>n</math> <math>\left(\beta_{ij}\right.</math>, где <math>j=1,2, \ldots,n</math> — координаты вектора <math>\mathbf{b_i}</math>, для <math>\left.i=1,2, \ldots,k\right)</math>, причем <math>\mathbf{b_1}=\mathbf{a_1}</math>, либо

<math>k</math> ортонормированных векторов <math>\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}</math> длины <math>n</math> (<math>\gamma_{ij}</math>, <math>j=1,2, \ldots,n</math> — координаты вектора <math>\mathbf{q_i}</math>, для <math>\left.i=1,2, \ldots,k\right)</math>, причем <math>\mathbf{q_1}=\frac{\mathbf{b_i}}{|\mathbf{b_i}|}</math>, для <math>i=1,2, \ldots,k</math>.

'''Последовательность исполнения метода''':

1. <math>len(a_1)=\sqrt{\alpha_{11}^2+\alpha_{12}^2+\ldots+\alpha_{1n}^2}</math>.

2. При <math>i</math> от <math>1</math> до <math>n</math>:

::<math>\gamma_{1j}= \frac{\alpha_{1j}}{len(a_1)}</math>.

3. При <math>i</math> от <math>2</math> до <math>k</math>:

:::при <math>j</math> от <math>1</math> до <math>n</math>:

:::::<math>\beta_{ij}=\alpha_{ij}-\sum\limits_{k=1}^{j-1}\left(\sum\limits_{m=1}^n \alpha_{im}\gamma_{km}\right)\gamma_{kj} </math>
::: <math>len(b_i)=\sqrt{\beta_{i1}^2+\beta_{i2}^2+\ldots+\beta_{in}^2}</math>

:::при <math>j</math> от <math>1</math> до <math>n</math>:
::::::<math>\gamma_{ij}= \frac{\beta_{ij}}{len(b_{i})}</math>.

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

Для ортогонализации <math>n</math> линейно независимых векторов длины <math>n</math> в последовательном варианте с помощью классического метода ортогонализации Грама-Шмидта потребуется:

* <math>k</math> вычислений квадратного корня
* <math>\frac{k(k-1)(n-1)}{2}</math> сложений,
* <math>\frac{k(k-1)(n-1)}{2}</math> вычитаний,
* <math>k^2 n</math> умножений,
* <math>k</math> делений.

Умножения, сложения и вычитания составляют ''основную часть алгоритма''.

При классификации по последовательной сложности, таким образом, метода ортогонализации Грама-Шмидта ''с кубической сложностью''.

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

Представим граф алгоритма с помощью текстового описания, а также с помощью изображения.

Все вершины графа можно разбить на пять групп, каждой из которых соответствует вид исполняемой операции. Каждая вершина расположена в точках с целочисленными координатами, в пространствах с числом измерений от <math>1</math> до <math>3</math>. Дадим описание каждой из групп.

'''Первая''' группа — вершины, которым соответствует операция скалярного умножения вектора на себя. На графе они располагаются в одномерной области. <math>i</math> изменяется от <math>1</math> до <math>k</math>.

Значениями аргументов для этих функций являются:
* при <math>i = 1</math> — входные данные <math>\{a_{1p}\}, p = 1..n</math>;
* при <math>i > 1</math> — результат выполнения функции в вершине пятой группы с координатами <math>\{i-1,1,p\}, p = 1..n</math>.

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

'''Вторая''' группа — вершины, которым соответствует операция SQRT извлечения квадратного корня. Они расположены в одномерной области, <math>i </math> изменяется от <math>1</math> до <math>k</math>.

Значением аргумента для этих функций является результат выполнения функции в вершине первой группы, с координатой <math>i</math>.

Результат является промежуточным данным алгоритма.

'''Третья''' группа – вершины, которым соответствует операция деления. Располагаются в двумерной области, <math>i</math> изменяется от <math>1</math> до <math>k</math>, <math>p</math> изменяется от <math>1</math> до <math>n</math>.

Значениями аргументов для этих функций являются:
* делимое:
** при <math>i = 1</math> — входные данные <math>a_{1p}</math>;
** при <math>i > 1</math> — результат выполнения функции в вершине пятой группы с координатами <math>i-1,1,p</math>;
* делитель - результат выполнения функции в вершине второй группы с координатой <math>i</math>.

Результат является выходным данным алгоритма <math>q_{ip}</math>.

'''Четвертая''' группа — вершины, которым соответствует операция скалярного умножения векторов <math>\mathbf{a}\cdot\mathbf{b}</math>. Располагаются в двумерной области, <math>i</math> изменяется от <math>1</math> до <math>k-1</math>, <math>j</math> изменяется от <math>1</math> до <math>k-i</math>.

Значениями аргументов для этих функций являются:
* <math>\mathbf{a}</math> — входные данные <math>\{a_{1+j,p}\},p = 1..n </math>;
* <math>\mathbf{b}</math> — результат выполнения функции в вершине третьей группы с координатой <math>\{i,p\}, p = 1..n </math>.

Результат является промежуточным данным алгоритма.

'''Пятая''' группа – вершины которые соответствуют операции <math>a - b c</math>. Располагаются в трехмерной области, <math>i</math> изменяется от <math>1</math> до <math>k-1</math>, <math>j</math> изменяется от <math>1</math> до <math>k-i</math>, <math>p</math> изменяется от <math>1</math> до <math>n</math>.

Значениями аргументов для этих функций являются:
* <math>a</math>:
** при <math>i = 1</math> — входные данные <math>a_{1+j,p}</math>;
** при <math>i > 1</math> — результат выполнения функции в вершине пятой группы с координатами <math>i-1,j+1,p</math>;
* <math>b</math> — результат выполнения функции в вершине четвертой группы, с координатами <math>i,j</math>;
* <math>c</math> — результат выполнения функции в вершине третьей группы с координатой <math>i,p</math>.

Результат является промежуточным данным алгоритма.

Описанный граф для случая <math>k = 3</math>, <math>n = 3</math> изображен на рисунке. Каждая группа вершин отличается цветом и буквенным обозначением. "S" — первая группа, "SQ" — вторая группа, "Div" — третья группа, "Dot" — четвертая группа, "f" — пятая группа. Квадратные метки "In" и "out" обозначают передачу входных и выходных данных соответственно.
[[File:GSParallelDiagram.png|thumb|center|1000px|Рисунок 1. Граф алгоритма с отображением входных и выходных данных]]

== Ресурс параллелизма алгоритма ==
Для построения ортонормированного базиса методом Грама–Шмидта необходимо выполнить следующие ярусы операций:
*<math>k</math> ярусов вычисления скалярных произведений векторов на самих себя (по одной операций на каждом, сложность операции по высоте и по ширине ЯПФ - корень квадратный);
*<math>k</math> ярусов вычисления квадратного корня (по одной операции на каждом);
*<math>k</math> ярусов деления (по <math>n</math> операций на каждом);
*<math>k-1</math> ярус вычисления скалярных произведений пары векторов (от <math>1</math> до <math>(k-1)</math> операций на каждом,сложность операции по высоте и по ширине ЯПФ - корень квадратный);
*<math>k-1</math> ярус умножения/вычитания чисел (от <math>n</math> до <math>(k-1)n</math> операций на каждом).
В этом случае сложность алгоритма при классификации по высоте ЯПФ — <math>O(k\sqrt{n})</math>, при классификации по ширине ЯПФ — <math>O(kn)</math>.

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

'''Входные данные''': матрица <math>A</math> с элементами <math>\alpha_{ij}</math>,
где <math>i</math> принимает значения <math>1,\ldots,k</math>, <math>j</math> принимает значения <math> 1,\ldots,n</math>.
Здесь <math>k\leqslant n</math>. По строкам этой матрицы записаны <math>k</math>
линейно-независимых векторов некоторого пространства размерности <math>n</math>.

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

'''Выходные данные''': матрица <math>Q</math> с элементами <math>\gamma_{ij}</math>. Строки этой матрицы также задают <math>k</math> векторов единичной длины размерности <math>n</math>, причем всякое скалярное произведение двух различных векторов этой системы равно <math>0</math>.

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

Замечание: Алгоритм может работать с вырожденной матрицей, он выдаёт нулевую строку на шаге <math>j</math>,
если строка <math>j</math> матрицы <math>A</math> является линейной комбинацией предыдущих строк.
В этом случае для корректного продолжения работы алгоритма необходима
дополнительная проверка на равенство строки нулю и пропуск соответствующих строк.
Количество построенных векторов будет равняться размерности пространства, порождаемого строками матрицы.

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

Отношение последовательной и параллельной сложности алгоритма в случае неограниченных ресурсов равно <math>\dfrac{k^2n}{k\sqrt{n}}=k\sqrt{n}</math>.
Вычислительная мощность алгоритма линейная.

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

Алгоритм является численно неустойчивым — ошибки округления могут привести к неортогональности полученных векторов.

Большинство дуг информационного графа являются локальными. Исключение составляют дуги исходящие из операций деления и извлечения корня — для первой группы количество дуг пропорционально квадрату количества векторов, а для второй — пропорционально размерности пространства, в котором заданы эти вектора. Возникающие при это длинные дуги могут быть заменены несколькими короткими пересылками между соседями.

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

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

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

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

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

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

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

== Существующие реализации алгоритма ==
Для классического метода ортогонализации Грама–Шмидта существуют реализации в математических программных пакетах, таких как
:*[https://www.mathworks.com/help/matlab/ref/cgs.html Matlab],
:*[http://maxima.sourceforge.net/docs/manual/maxima_23.html#gramschmidt Maxima],
:*[https://reference.wolfram.com/language/ref/Orthogonalize.html# Mathematica].

Также метод реализован в библиотеке [http://www.spectralpython.net/class_func_ref.html#spectral.algorithms.algorithms.orthogonalize SpectralPython] для языка Python. Для устранения проблемы ошибок численных округлений иногда используется повторное применение метода к уже построенному ортогональному базису.

Поскольку метод ортогонализации Грама–Шмидта гарантирует получение <math>j</math> ортогональных векторов к концу шага <math>j</math>, а не в самом конце работы алгоритма, используется в итерационных методах, таких как метод Арнольди. Для поддержки работы этого метода реализация ортогонализации Грама–Шмидта включена в библиотеку [http://trac.alcf.anl.gov/projects/libs/browser/ARPACK/PARPACK/SRC/MPI/pzgetv0.f/ ARPACK].

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

[[en: Gram–Schmidt Orthogonalization]]

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

2016-10-26T13:27:04Z

Sergeys79:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2016-10-26T13:26:42Z

Sergeys79:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2016-10-26T13:26:26Z

Sergeys79:

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

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

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

=== Математическое описание алгоритма ===
Исходные данные: положительно определённая симметрическая матрица <math>A</math>, вектор <math>b</math>, количество итераций <math>k</math>.

Вычисляемые данные: трехдиагональная матрица <math>T_k</math>(элементы <math>t_{ij}</math>) размерности <math>k</math>.

Итак, пусть <math>A=A^t \in R^{n{\times}n}</math>, тогда

:<math>T_k</math><math>=Q_k^tAQ_k, </math>

где <math>T_k</math>-верхняя симметричная трехдиагональная матрица, <math>Q_k=[q_1,q_2,...,q_k]</math>, столбцы которой образуют
ортонормированный базис подпространства

Крылова <math>K_m(A,b)</math> для заданной матрицы <math>A</math> и вектора <math>b</math>.

Обозначим:

<math>\alpha_i=t_{i,i},i=1,k </math>

<math>\beta_i=t_{i-1,i},i=2,k</math>

Тогда матрица <math>T_k=\begin{bmatrix}
\alpha_1 & \beta_2 \\
\beta_2 & \alpha_1 & \beta_3 &\\
&. & . & .\\
&&\beta_{k-1} & \alpha_{k-1} & \beta_k\\
&&&\beta_k & \alpha_k

\end{bmatrix}</math>,

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

:<math>
\begin{align}
q_1&=b/||b||,\beta_0=0,q_o=0. \\
z&=Aq_j, \\
\alpha_j&=q_j^Tz, \\
z&=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, \\
\beta&=||z||,\\
q_{j+1}&=z/\beta_j, \quad j \in [1, k].
\end{align}
</math>
В точной арифметике на <math>j</math>-ой итерации <math>z</math> будет ортогонален векторам <math>q_1,q_2..q_{j-1}</math>, но к сожалению это свойство теряется из-за округлений в арифметике с плавающей точкой. Чтобы бороться с этим, проводится полная переортоганализация, что соответствует повторному проведению процесса ортогонализации Грамма-Шмидта <math>z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i.</math> Однако добавление этих двух операций делает реализацию алгоритма Ланцоша гораздо дорогостоящей.

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

*1. <math>Aq=( \sum\nolimits_{i=^n}a_{1i}q_i, \sum\nolimits_{i=2}^na_{2i}q_i, ..., \sum\nolimits_{i=1}^na_{ni}q_i)</math>, как известно в операции задействуется <math>n^2</math> умножений и <math>n(n-1)</math>сложений.

*2. <math>z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i,
z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i,</math> сложность переортогонализации:<math>k(k+1)n+2k</math> умножений и <math>k(k+1)n+2+k</math> сложений.

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

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

:<math>(x,y)=\sum\nolimits_i^nx_iy_i</math>

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

=== Схема реализации последовательного алгоритма ===
Последовательность исполнения метода следующая:

<math>
\begin{align}
\\
1 .& q_1=b/||b||,\\
2 .& \beta_0=0,\\
3 .&q_o=0. \\
\end{align}
</math>

Далее для всех <math>j</math> от <math>1</math> до <math>k</math> по нарастанию выполняются

<math>
\begin{align}
4&.z=Aq_j, \\
5&.\alpha_j=q_j^Tz, \\
6&.z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i, \\
7&.z=z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i, \\
8&.z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, \\
9&.\beta=||z||,\\\end{align}

</math>

Если <math>\beta=0</math>, то алгоритм завершается,

<math>10.q_{j+1}=z/\beta_j.</math>

=== Последовательная сложность алгоритма ===
Для выполнения алгоритма Ланцоша в последовательном варианте требуется:

*<math>k</math> вычислений квадратного корня,

*<math>k((n+1)^2+nk-1)</math> сложений (вычитаний),

*<math>k((n+1)^2+2kn)</math> умножений

Умножения и сложения (вычитания) составляют основную часть алгоритма.

=== Информационный граф===
Опишем граф алгоритма в виде рисунка.
* Вершина черного цвета соответствует операции <math>Aq_j,</math>.
* Вершина красного цвета соответствует операции скалярного произведения <math>q_j^Tz</math>.
* Вершина сиреневого цвета соответствует операции <math>z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i</math>.
* Вершина желтого цвета соответствует операции <math>z-\alpha_jq_j-\beta_{j-1}q_{j-1}</math>.
* Вершина синего цвета соответствует операции <math>||z||</math>.
* Вершина зеленого цвета соответствует операции <math>z/\beta_j</math>.
[[Файл:photo.jpg]]Здесь показаны первые три итерации алгоритма, при больших итерациях граф строится аналогичным образом.

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

=== Ресурс параллелизма алгоритма ===
Как видно из графа алгоритма каждая операция передает результат своих вычислений последующей операции. Соответственно никаких ресурсов параллелизма на уровне взаимодействия этих базовых операции нету. Они могут выполнятся строго последовательно. Однако, шесть операций содержат вычисления множества скалярных произведений.Стоит отметить, что в каждой итерации есть операция умножения матрицы на вектор, что представляет из себя <math>n</math> информационно независимых скалярных произведений. Также хорошие возможности для параллельной реализации представляет повторная переортогонализация, так как на <math>j</math> итерации обе операции ортогонализации вычисляют <math>j-1</math> независимых скалярных произведения. Поскольку каждое вычисление скалярного произведения векторов длины <math>n</math> требует выполнения <math>n</math> операций умножения и <math>n-1</math> операций сложения, его сложность <math>n</math>. Если воспользоваться разделением данных по строкам, то в случае параллельных вычислений каждый процессор производит умножение только части (полосы) первого вектора на второй, размер этих полос равен <math>n/p</math>

(<math>p</math>-количество процессоров). Следовательно, вычислительная сложность параллельного алгоритма скалярного произведения определяется выражением:<math>n/p</math>, а мультипликативная сложность всего алгоритма <math>O(n^2k+nk^2)/p)</math>. При <math>p=n</math> получаем <math>O(nk+k^2)</math>. Однако этот метод не учитывает, то что есть скалярные произведения, которые могут вычисляться параллельно. Если учесть весь ресурс параллелизма можно получить сложность <math>O(k)</math>.

=== Входные и выходные данные алгоритма ===
'''Входные данные''': плотная матрица <math>A</math> (элементы <math>a_{ij}</math>),произвольный вектор ненулевой <math>b</math>, количество итераций <math>k</math>.

Дополнительные ограничения:
* <math>A</math> – симметрическая матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>..
'''Объём входных данных''': <math>\frac{n (n + 3)}{2}</math> (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы).

'''Выходные данные''': трехдиагональная симметрическая матрица <math>T</math> (элементы <math>_{ij}</math>), либо два вектора <math>a</math> и <math>b</math> размерности <math>k</math>.

'''Объём выходных данных''': <math>2k</math> (в силу трехдиагональности симметричности достаточно хранить только ненулевые элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом.

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

Количество операций умножения и количество операций сложения имеют один и тот же порядок, но они явно преобладают над операцией взятия корня. Также сбалансированы между собой арифметические операции и операции обращений в память. Алгоритм обрабатывает матрицу и вектора и считает множество скалярных произведений, следовательно обеспечивается хорошая сбалансированность операций между параллельными ветвями алгоритма.
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов равняется <math>O(n^2)</math>.

При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных равняется <math>O(k)</math>. Алгоритм недетерминированный, так как если в какой-то итерации <math>b_j=0</math>, то алгоритм останавливается. Также на недетерминированность может повлиять изменение порядка выполнения ассоциативных операций, в частности при вычислении скалярного произведения.

"Длинных" дуг в информационном графе нет, запись данных происходит предсказуемо(на каждом шаге заполняются следующие элементы выходных векторов), следовательно нет никаких потенциальных сложностей с размещением данных в иерархии памяти компьютера на этапе выполнения программы.

Основным недостатком алгоритма заключается в том, что высокая цена которую мы платим за переортогонализацию, ведь за <math>k</math> итераций алгоритм потребуют <math>O(k^2n)</math> флопов вместо <math>O(kn)</math> и <math>O(kn)</math> слов памяти вместо <math>O(n)</math> при реализации без ортогонализации.

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

=== Существующие реализации алгоритма ===
NAG Library имеет множество процедур для решения система алгебраических уравнений и нахождения собственных значений, которые используют алгоритм Ланцоша. MATLAB содержит реализацию алгоритма в пакете Gaussian Belief Propagation Matlab Package. GraphLab имеет библиотеку, в которой есть огромное количество параллельных реализаций алгоритма Ланцоша на языке С++. Также в библиотеках PRIMME есть реализации алгоритмов подобных алгоритму Ланцоша.
==Литература ==
ДЖ.Деммель Вычислительная линейная алгебра. М.:Мир,2001

В. В. Вербицкий, В. В. Реут ВВЕДЕНИЕ В ЧИСЛЕННЫЕ МЕТОДЫ АЛГЕБРЫ 2015

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

2016-10-26T08:37:32Z

Sergeys79:

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (2)]] ==

* Нет краткой сводки свойств алгоритма в начале страницы
* Нумерация подразделов не соответствует шаблону
* Разделы 1.6 и 1.8 (6 и 8 в оригинальной статье) не заполнены
* Картинка в разделе 1.7 (7) не является информационным графом
* Раздел 2.7 отсутствует
* Список литературы отсутствует

Обсуждение участницы:Sannikovats

2016-10-26T08:37:10Z

Sergeys79:

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (1)]] ==

* Не указан вклад обоих авторов
* Информационный граф выглядит странно

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

2016-10-26T08:36:50Z

Sergeys79:

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки]] ==

* Нет краткой сводки свойств алгоритма в начале страницы
* В разделе 1.1 по сути нет описания алгоритма
* Разделы 1.6-1.8 пусты
* В разделах 1.3 и 1.4 хотелось бы видеть понятное словесное описание
* Содержимое раздела 1.4 должно быть внесено в раздел 1.5
* Не указан вклад каждого автора

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

2016-10-26T08:36:05Z

Sergeys79: Содержимое страницы заменено на «== Статья Ортогонализация Грама-Шмидта ==»

== Статья [[Ортогонализация Грама-Шмидта]] ==

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

2016-10-26T08:34:55Z

Sergeys79:

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

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

2016-10-26T08:34:22Z

Sergeys79:

== Статья [[Алгоритм Ланцоша в арифметике с плавающей точкой.]] ==
* Точку в конце названия статьи лучше убрать
* П. 2.1. должен называться 2.7
* Не приведена краткая сводка свойств алгоритма в начале статьи

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

2016-10-25T13:22:18Z

Sergeys79: /* Статья Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (2) */

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (2)]] ==

* Статья размещена не в личном пространстве автора
* Нет краткой сводки свойств алгоритма в начале страницы
* Нумерация подразделов не соответствует шаблону
* Разделы 1.6 и 1.8 (6 и 8 в оригинальной статье) не заполнены
* Картинка в разделе 1.7 (7) не является информационным графом
* Раздел 2.7 отсутствует
* Список литературы отсутствует

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

2016-10-25T13:20:31Z

Sergeys79: Новая страница: «== Статья Вычисление определенного интеграла с использованием адаптивно сгущающейся с…»

Обсуждение участницы:Sannikovats

2016-10-25T13:16:37Z

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки (1)]] ==

* Статья размещена не в личном пространстве автора
* Не указан вклад обоих авторов
* Информационный граф выглядит странно

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

2016-10-25T13:15:58Z

Sergeys79: /* Статья Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки */

== Статья [[Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки]] ==

* Статья размещена не в личном пространстве автора
* Нет краткой сводки свойств алгоритма в начале страницы
* В разделе 1.1 по сути нет описания алгоритма
* Разделы 1.6-1.8 пусты
* В разделах 1.3 и 1.4 хотелось бы видеть понятное словесное описание
* Содержимое раздела 1.4 должно быть внесено в раздел 1.5
* Не указан вклад каждого автора

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

2016-10-25T13:10:11Z

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

2016-10-25T12:48:50Z

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

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

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

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

2016-10-25T12:47:32Z

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

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

* В разделе 1.3 не приведено описание вычислительного ядра
* Не указано ФИО автора
* Не приведена краткая сводка свойств алгоритма в начале статьи

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

2016-10-25T12:47:09Z

Sergeys79:

== Статья [[Алгоритм Ланцоша в арифметике с плавающей точкой.]] ==
* Точку в конце названия статьи лучше убрать
* Статья не в личном пространстве автора
* П. 2.1. должен называться 2.7
* Не приведена краткая сводка свойств алгоритма в начале статьи

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

2016-10-25T12:43:10Z

Sergeys79: Новая страница: «== Статья Участник:DVIN/Ортогонализация Грамма-Шмидта == * В разделе 1.3 определение вычисл…»

== Статья [[Участник:DVIN/Ортогонализация Грамма-Шмидта]] ==

* В разделе 1.3 определение вычислительного ядра можно опустить

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

2016-10-25T12:38:43Z

Sergeys79: Новая страница: «== Статья Ортогонализация Грама-Шмидта == * Статья не в личном пространстве автора * Стр…»

== Статья [[Ортогонализация Грама-Шмидта]] ==

* Статья не в личном пространстве автора
* Странное содержание в подразделе 2.4

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

2016-10-25T12:31:09Z

Sergeys79: Новая страница: «== Статья Участник:Yastrebovks/Алгоритм Ланцоша с полной переортогонализацией == * Неправиль…»

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

* Неправильная нумерация подразделов, искажены названия нескольких подразделов
* Указано два автора, но не указан вклад каждого
* В информационном графе не ясно, что означают троеточия

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

2016-10-25T12:25:22Z

Sergeys79:

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

2016-10-25T12:23:49Z

Sergeys79:

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

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

2016-10-25T12:23:33Z

Sergeys79: Новая страница: «== Статья Алгоритм Ланцоша для арифметики с плавающей точкой с полной переортогонализа…»

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

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

2016-10-25T12:21:28Z

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

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

* Статья размещена не в личном пространстве автора

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

2016-10-25T12:21:10Z

Sergeys79: Новая страница: «== Статья Алгоритм Ланцоша с полной переортогонализацией — история изменений == * Стат…»

== Статья [[Алгоритм Ланцоша с полной переортогонализацией — история изменений]] ==

* Статья размещена не в личном пространстве автора

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

2016-10-25T12:18:05Z

Sergeys79: Новая страница: «== Статья Алгоритм Ланцоша в арифметике с плавающей точкой. == * Точку в конце названия с…»

== Статья [[Алгоритм Ланцоша в арифметике с плавающей точкой.]] ==
* Точку в конце названия статьи лучше убрать
* Статья не в личном пространстве автора
* П. 2.1. должен называться 2.7
* Раздел 2 представлен единственным пунктом