Участник: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{ \frac{\lt a, b\gt }{\|b\|^2} * b }[/math] как отдельную макрооперацию [math]\displaystyle{ proj_{b}a }[/math].

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

Вычислительным ядром алгоритма является вычисление макроопераций [math]\displaystyle{ proj_{b}a }[/math], определенных выше.

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

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

b[1] = a[1]
for j = 2,...,n
    b[j] = a[j]
    for k = 1,...,j-1
        b[j] = b[j] - <a[j],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{ 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/Процесс_Грама_―_Шмидта