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

Файл:Symbol wait.svg

Эта работа прошла предварительную проверку Дата последней правки страницы: 23.11.2016 Данная работа соответствует формальным критериям. Проверено ASA.

Ортогонализация Грама-Шмидта
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(N^2k) }[/math]
Объём входных данных	[math]\displaystyle{ Nk }[/math]
Объём выходных данных	[math]\displaystyle{ Nk }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(N) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(Nk) }[/math]

Основные авторы описания: Инжелевская Дарья Валерьевна

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

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

Создатели алгоритма Эрхард Шмидт и Йорген Педерсен Грам.

Эрхард Шмидт (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{ {\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}} }[/math] строится множество ортогональных векторов [math]\displaystyle{ {\displaystyle \mathbf {b}_{1},\;\ldots ,\;\mathbf {b} _{N}} }[/math] или ортонормированных векторов [math]\displaystyle{ {\displaystyle \mathbf {e} _{1},\;\ldots ,\;\mathbf {e}_{N}} }[/math], причём так, что каждый вектор [math]\displaystyle{ {\displaystyle \mathbf {b} _{j}} }[/math] или [math]\displaystyle{ {\displaystyle \mathbf {e} _{j}} }[/math] может быть выражен линейной комбинацией векторов [math]\displaystyle{ {\displaystyle \mathbf {a} _{1},\;\ldots ,\; \mathbf {a} _{j}} }[/math].

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

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

Алгоритм

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

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

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

Скалярное произведение для двух векторов [math]\displaystyle{ \mathbf{ a= [a_1, a_2, ...,a_k]} }[/math] и [math]\displaystyle{ \mathbf{ b= [b_1, b_2, ..., b_k]} }[/math] в k-мерном действительном пространстве определяется как:

[math]\displaystyle{ \langle \mathbf{a}, \mathbf{b} \rangle=\sum_{i=1}^k a_ib_i=a_1b_1+a_2b_2+\cdots+ a_kb_k }[/math].

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

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

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

[math]\displaystyle{ {\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]\displaystyle{ \mathbf{b}_j \;(j = 1 \ldots N) }[/math] может быть получен нормированный вектор: [math]\displaystyle{ \mathbf{e}_j = {\mathbf{b}_j\over \| \mathbf{b}_j \|} }[/math] (у нормированного вектора направление будет таким же, как у исходного, а длина — единичной). Причем берется норма согласованная со скалярным произведением [math]\displaystyle{ \| x \| = \sqrt{\langle x, x \rangle} }[/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}_1,\;\ldots,\;\mathbf{b}_N }[/math] носит название ортогонализации Грама — Шмидта, а [math]\displaystyle{ \mathbf{e}_1,\;\ldots,\;\mathbf{e}_N }[/math] — ортонормализации Грама — Шмидта.

Доказательство

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

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

Геометрическая интерпретация — вариант 1

Рис. 1. Второй шаг процесса Грама — Шмидта

Рассмотрим формулу (2) — второй шаг алгоритма. Её геометрическое представление изображено на рис. 1:

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

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

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

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

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

Рассмотрим формулу (3) — третий шаг алгоритма — в следующем варианте:

[math]\displaystyle{ \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:

Рис. 2. Третий шаг процесса Грама — Шмидта

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

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

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

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

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

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

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

Численная неустойчивость

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

Модифицированный процесс Грама — Шмидта

Процесс Грама — Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Вместо вычисления [math]\displaystyle{ \mathbf{b}_j }[/math] как: [math]\displaystyle{ \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]\displaystyle{ \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]

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

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

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

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

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

Следующий алгоритм реализует нормализацию Грамма-Шмидта. Векторы [math]\displaystyle{ v_1,...,v_k }[/math] заменяются набором ортонормированных векторов, которые имеют ту же линейную оболочку.

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

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

1. [math]\displaystyle{ \mathbf {b}_{1}=\mathbf {a}_{1} }[/math]

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

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

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

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

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

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

При k=N:

• [math]\displaystyle{ \frac{N(N-1)}{2} }[/math] делений,

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

• [math]\displaystyle{ \frac{3N^2 (N-1)}{2} }[/math] умножений.

При k≠N

• [math]\displaystyle{ \frac{N(N-1)}{2} }[/math] делений,

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

• [math]\displaystyle{ \frac{3Nk (N-1)}{2} }[/math] умножений.

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

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

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

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

Первая группа вершин расположена в двумерной области, соответствующая ей операция [math]\displaystyle{ \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]\displaystyle{ a_i }[/math]: элементы входных данных, а именно [math]\displaystyle{ a_i }[/math];

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

Вторая группа вершин расположена в двухмерной области, соответствующая ей операция [math]\displaystyle{ a - b }[/math].

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

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

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

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

a:

j=1 - входные данные [math]\displaystyle{ a_j }[/math]

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

b:

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

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

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

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

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

• [math]\displaystyle{ N }[/math] ярусов в каждом из которых будет выполнятся,
• по [math]\displaystyle{ N - 1 }[/math] операторов проекций (в каждом из ярусов линейное количество проекций, в зависимости от яруса — от [math]\displaystyle{ N-1 }[/math] до [math]\displaystyle{ 1 }[/math]),
• по [math]\displaystyle{ N - 1 }[/math] сложений (в каждом из ярусов линейное количество сложений, в зависимости от яруса — от [math]\displaystyle{ N-1 }[/math] до [math]\displaystyle{ 1 }[/math]),

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

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

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

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

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

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

Объём входных данных: [math]\displaystyle{ Nk }[/math]

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

Объём выходных данных: [math]\displaystyle{ Nk }[/math]

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

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

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

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

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

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

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

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

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

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

Проведём исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Исследовалась собственная параллельная реализация алгоритма, написанная на языке C++ с использованием стандарта MPI.

Версия MPI - intel opemmpi/1.8.4-icc

Модель используемого CPU - Intel Xeon X5570 2.93GHz

Набор и границы значений изменяемых реализации алгоритма:

• число процессоров [10:100] с шагом 10;
• размерность векторов [50 : 500] с шагом 50.

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

Рис. 6 Время работы алгоритма в зависимости от числа процессоров и размерности задачи

Так как алгоритм не поддается полному распараллеливанию, то при увеличении числа процессоров не происходит значительного ускорения. Это связано с тем, что параллельно можно вычислить лишь [math]\displaystyle{ b_j }[/math]. То есть после каждого вычисления [math]\displaystyle{ b_j }[/math] можно запускать параллельное вычисление [math]\displaystyle{ \mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i} }[/math] для всех [math]\displaystyle{ i=j+1..N }[/math]. Кроме того можно параллельно производить вычитание соответствующего [math]\displaystyle{ \mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i} }[/math] для всех [math]\displaystyle{ i=j+1..N }[/math]. Но после этого все процессоры должны синхронизоваться и дождаться когда получат следующий вектор [math]\displaystyle{ b_j }[/math].

Проведём еще одно исследование исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта для меньшего числа процессоров. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Исследовалась параллельная реализация алгоритма, написанная с использованием стандарта MPI.

Набор и границы значений изменяемых реализации алгоритма:

• число процессоров [1:10] с шагом 1;
• размерность векторов [20 : 400] с шагом 20.

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

Рис. 7 Время работы программы в зависимости от числа процессоров и размерности задачи

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

Рис. 8 Производительность программы в зависимости от количество процессоров и размерности задачи

На графике 3 приведена эффективность распараллеливания работы программы для разных значений размерности матрицы и числа процессоров.

Рис. 9 Эффективность распараллеливания алгоритма

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

• По числу процессоров: -0.009265. С увеличением числа процессоров наблюдается слабое снижение эффективности работы параллельной реализации алгоритма.

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

• По размеру задачи: 0.00002555. При увеличении размерности задачи эффективность алгоритма очень медленно возрастает.

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

• По двум направлениям: 0.002079.

Реализация алгоритма

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

Примеры библиотек, содержащих реализации алгоритма ортогонализации Грама-Шмидта:

1. C#
2. Matlab
3. Wolfram
4. NAG Fortran

Реализация для пакета Mathematica

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

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}}]

Реализация для Maxima

Данный скрипт, предназначенный для пакета Maxima, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в квадратных скобках. Количество векторов и их координат могут быть произвольными.

В данном случае для примера взяты векторы [-2,1,0], [-2,0,1], [-0.5,-1,1].

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

Реализация для 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/Процесс_Грама_―_Шмидта Процесс Грама ― Шмидта