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

Авторы описания: Шерварлы В.Г.

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

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

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

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

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

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

Процесс ортогонализации описывается следующей последовательностью формул:

[math]\displaystyle{ \begin{align} b_1 & = a_1 \\ b_2 & = a_2 - \frac{\lt a_2, b_1\gt }{\|b_1\|^2} * b_1 \\ b_3 & = a_3 - \frac{\lt a_3, b_1\gt }{\|b_1\|^2} * b_1 - \frac{\lt a_3, b_2\gt }{\|b_2\|^2} * b_2 \\ ... \\ b_n & = a_n - \frac{\lt a_n, b_1\gt }{\|b_1\|^2} * b_1 - ... - \frac{\lt a_n, b_{n-1}\gt }{\|b_{n-1}\|^2} * b_{n-1} \\ \end{align} }[/math]

где: [math]\displaystyle{ a_1, ..., a_n }[/math] - исходное множество линейно-независимых векторов

[math]\displaystyle{ b_1, ..., b_n }[/math] - множество искомых векторов

[math]\displaystyle{ \lt a,b\gt }[/math] - скалярное произведение векторов

[math]\displaystyle{ \|b\| = \sqrt{\lt b,b\gt } }[/math]

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

Вычислительным ядром алгоритма является вычисление оператора проекции [math]\displaystyle{ proj_{b}a = \frac{\lt a, b\gt }{\|b\|^2} * b }[/math].

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

На уровне макроструктуры будем рассматривать [math]\displaystyle{ \frac{\lt a, b\gt }{\|b\|^2} * b }[/math] как отдельную макрооперацию [math]\displaystyle{ proj_{b}a }[/math].

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

Шаг 1:

[math]\displaystyle{ b_1 = a_1 }[/math]

Шаг i (i = 2, ..., n, где n - количество векторов):

[math]\displaystyle{ b_i = a_i - \sum^{i-1}_{k=1}{proj_{b_k}a_i} }[/math]

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

Определим входные и выходные данные (два множества векторов, n - количество векторов в множестве):

input: a = (a[1], ..., a[n])
output: b = (b[1], ..., b[n])

И далее последовательно находим каждый из векторов [math]\displaystyle{ b_i }[/math], используя полученные на предыдущих итерациях цикла векторы [math]\displaystyle{ b_1, ..., b_{i-1} }[/math]:

b[1] = a[1]
for i = 2,...,n
    b[i] = a[i]
    Comment: во вложенном цикле вычисляем [math]\displaystyle{ b_i = a_i - \sum^{i-1}_{k=1}{proj_{b_k}a_i} }[/math]
    for k = 1,...,i-1
        b[i] = b[i] - <a[i],b[k]>/<b[k],b[k]>*b[k]     

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

Последовательная сложность алгоритма выражается следующими формулами ([math]\displaystyle{ n }[/math] - количество векторов, [math]\displaystyle{ m }[/math] - количество компонент вектора):

• количество умножений: [math]\displaystyle{ \frac{3mn(n-1)}{2} }[/math]

• количество делений: [math]\displaystyle{ \frac{n(n-1)}{2} }[/math]

• количество сложений и вычитаний: [math]\displaystyle{ \frac{n(n-1)(2m - 1)}{2} }[/math]

Общая сложность алгоритма - [math]\displaystyle{ O(mn^2) }[/math]

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

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

• [math]\displaystyle{ a_i, i=1, ..., n }[/math] - операции присваивания [math]\displaystyle{ b_i = a_i }[/math]
• [math]\displaystyle{ s_{ij} = b_i - proj_{b_j}a_i }[/math]
• [math]\displaystyle{ a_i^{(*)} }[/math] и [math]\displaystyle{ f_{ij} }[/math] отличаются от [math]\displaystyle{ a_i }[/math] и [math]\displaystyle{ s_{ij} }[/math] соответственно тем, что в этот момент завершается вычисление вектора [math]\displaystyle{ b_i }[/math], являющегося и частью выходных данных, и необходимой информацией для последующего этапа. Так же в этих вершинах можно вычислять величину [math]\displaystyle{ \|b_i\|^2 }[/math], необходимую для вычисления соответствующих [math]\displaystyle{ proj_{b}a }[/math].
• Разными оттенками зеленого цвета показана ярусно-параллельная форма графа алгоритма (вершины одного яруса выделены одним цветом). Высота ЯПФ - [math]\displaystyle{ n }[/math], ширина ЯПФ - [math]\displaystyle{ n }[/math].

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

Возможно параллельное вычисление всех [math]\displaystyle{ proj_{b_j}a_i, i=1,...n }[/math] для уже найденных [math]\displaystyle{ b_j }[/math] вместо вычисления каждого вектора [math]\displaystyle{ b_j }[/math] по очереди. Таким образом после каждой итерации будет получен один посчитанный вектор, необходимый для расчетов на следующей итерации. Однако степень распараллеливания будет падать на каждой итерации, что вместе с необходимостью синхронизации и многочисленными обменами информации может снизить эффективность распараллеливания.

Сложность параллельного варианта равна [math]\displaystyle{ O(nm) }[/math], где [math]\displaystyle{ n }[/math] - количество векторов, [math]\displaystyle{ m }[/math] - количество компонент вектора. Можно сократить эту величину до [math]\displaystyle{ O(n\log{m}) }[/math], если еще параллельно проводить вычисление каждой макрооперации [math]\displaystyle{ proj_{b}a }[/math], однако, опять же, не стоит забывать о расходах на синхронизацию и обмен информацией.

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

Входные данные: N линейно независимых векторов (процесс Грама — Шмидта может применяться также к бесконечной последовательности линейно независимых векторов, а также к линейно зависимым векторам. В последнем случае вектор [math]\displaystyle{ a_j }[/math] получится нулевым, если он зависит от векторов [math]\displaystyle{ a_1, ..., a_{j-1} }[/math], и алгоритм должен сразу отбрасывать нулевые векторы). Объем входных данных: [math]\displaystyle{ mN }[/math], где [math]\displaystyle{ N }[/math] - количество векторов, [math]\displaystyle{ m }[/math] - количество компонент вектора

Выходные данные: N ортогональных векторов, линейная оболочка которых совпадает с линейной оболочкой входного множества векторов, объем выходных данных: [math]\displaystyle{ mN }[/math].

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

Вычислительная мощность алгоритма: [math]\displaystyle{ O(n) }[/math]. Отношение последовательной сложности к параллельной - [math]\displaystyle{ \frac{mn^2}{mn} }[/math] [math]\displaystyle{ ( }[/math]или [math]\displaystyle{ \frac{mn^2}{n\log{m}}) }[/math].

Классический алгоритм (CGS) является численно-неустойчивым из-за частой потери ортогональности, вызванной ошибками округления в процессе вычислений.

Существует модификации алгоритма Грама-Шмидта, которые являются более устойчивыми, например, модифицированный алгоритм Грама-Шмидта (MGS), блочные и итеративные варианты CGS и MGS. MGS численно эквивалентен методу Хаусхолдера QR-разложения, примененного к матрице, полученной из исходной путем добавления к ней сверху прямоугольной матрицы из нулевых элементов.

Как можно видеть на изображении информационного графа алгоритма, существуют вершины с большой степенью исхода. Эти вершины соответствуют (в том числе) вычислению вектора [math]\displaystyle{ b_i }[/math] и результат этого вычисления должен быть сообщен всем вершинам, вычисляющим [math]\displaystyle{ proj_{b_i}a_j }[/math] для [math]\displaystyle{ j=1,..., n }[/math].

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

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

• реализация для системы компьютерной алгебры maxima (GNU GPL)

load(eigen);
x: matrix ([-2,1,0],[-2,0,1],[-0.5,-1,1]);
y: gramschmidt(x);

• пример реализации для системы компьютерной алгебры mathematica (проприетарное программное обеспечение)

Projection[v1_, v2_] := (v1.v2*v2)/v2.v2
MultipleProjection[v1_, vecs_] := Plus @@ (Projection[v1, #1] &) /@ vecs
GramSchmidt[mat_] := Fold[Join[#1, {Normalize[#2 - MultipleProjection[#2, #1]]}] &, {}, mat]
GramSchmidt[{{-2, 1, 0}, {-2, 0, 1}, {-0.5, -1, 1}}]

• библиотека fplll

• Math.NET Numerics

• библиотека NAG

• библиотека PETSc

Литература

1. https://ru.wikipedia.org/wiki/Процесс_Грама_―_Шмидта
2. Канатников А.Н., Крищенко А.П. "Линейная алгебра."
3. Å. Björck. "Numerics of Gram-Schmidt orthogonalization"
4. B. Milde, M. Schneider. "Parallel Implementation of Classical Gram-Schmidt Orthogonalization on CUDA Graphics Cards"