Классический метод ортогонализации

Классическая ортогонализация Грама-Шмидта
Последовательный алгоритм
Последовательная сложность	$O(N^3)$
Объём входных данных	$N^2$
Объём выходных данных	$3N^2/2$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(N)$
Ширина ярусно-параллельной формы	$O(N^2)$

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

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

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

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

Пусть имеются линейно независимые векторы $\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N$.

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

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

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

${\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}}$

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

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

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

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

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

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

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

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

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

Следующий алгоритм реализует нормализацию Грамма-Шмидта. Векторы $v_1,...,v_k$ заменяются набором ортонормированных векторов, которые имеют ту же линейную оболочку.

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

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

1. $\mathbf {b}_{1}=\mathbf {a}_{1}$

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

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

1.5.1 Пример реализации на Python

Функция работает для произвольного количества векторов любой размерности. При этом если количество векторов $\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N$ больше их размерности или они линейно зависимы то функция возвращает максимально возможное число линейно независимых векторов $\mathbf{b}_1,\;\ldots,\;\mathbf{b}_n$, а остальные векторы $\mathbf{b}_{n+1},\;\ldots,\;\mathbf{b}_N$ нулевые.

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))

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

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

При k=N:

• $\frac{N(N-1)}{2}$ делений,

• $\frac{N(N-1)(3N-1)}{2}$ сложений (вычитаний),

• $\frac{3N^2 (N-1)}{2}$ умножений.

При k≠N

• $\frac{N(N-1)}{2}$ делений,

• $\frac{N(N-1)(3k-1)}{2}$ сложений (вычитаний),

• $\frac{3Nk (N-1)}{2}$ умножений.

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

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

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

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

Первая группа вершин расположена в двумерной области, соответствующая ей операция $\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} ,$ .

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

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

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

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

$a_i$: элементы входных данных, а именно $a_i$;

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

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

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

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

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

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

a:

j=1 - входные данные $a_j$

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

b:

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

Рис. 5. Граф алгоритма с отображением входных и выходных данных. Proj - вычисление оператора проекции, F - операция a-b, In - входные данные, Out - результаты.

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

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

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

• $N$ ярусов в каждом из которых будет выполнятся,
• по $N - 1$ операторов проекций (в каждом из ярусов линейное количество проекций, в зависимости от яруса — от $N-1$ до $1$),
• по $N - 1$ сложений (в каждом из ярусов линейное количество сложений, в зависимости от яруса — от $N-1$ до $1$),

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

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

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

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

Дополнительные ограничения:

• вектора ${\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}$ линейно независимы, поэтому $k \geqslant N$.

Объём входных данных: $Nk$

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

Объём выходных данных: $Nk$

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

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

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

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

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

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

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

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