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

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

Основные авторы описания: Белов Н. А., Богомолов Д. А..

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

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

Ортогонализация ― алгоритм построения для данной линейно независимой системы векторов евклидова или эрмитова пространства V ортогональной системы ненулевых векторов, порождающих то же самое подпространство в V.

Ортогонализация Грама–Шмидта (процесс Грама–Шмидта) ― наиболее известный алгоритм ортогонализации. Назван в честь Йоргена Педерсена Грама и Эрхарда Шмидта, однако ранее уже появлялся в работах Лапласа и Коши. Является частным случаем разложения Ивасавы, так как может быть представлен как разложение невырожденной квадратной матрицы в произведение ортогональной (или унитарной матрицы в случае эрмитова пространства) и верхнетреугольной матрицы с положительными диагональными элементами.

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

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

Входные данные. $m$ линейно независимых векторов $\mathbf{a}_1,...,\mathbf{a}_m$ c размерностью пространства $n$, записанных в матрице $A$ с элементами $\alpha_{ij}$.

Выходные данные. $m$ ортогональных векторов $\mathbf{b}_1,...,\mathbf{b}_m$ c размерностью пространства $n$, записанных в матрице $B$ с элементами $\beta_{ij}$.

Определим оператор проекции (проецирует вектор $\mathbf{a}$ коллинеарно вектору $\mathbf{b}$) $\mathbf{proj_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b}$, где $\mathbf{\left \langle a,b \right \rangle}$ ― скалярное произведение векторов $\mathbf{a}$ и $\mathbf{b}$.

Тогда ортогональные векторы вычисляются следующим образом:

$\mathbf{b}_1 = \mathbf{a}_1$

$\mathbf{b}_2 = \mathbf{a}_2 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_2$

$\mathbf{b}_3 = \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_2} \mathbf{a}_3$

$\vdots$

$\mathbf{b}_m = \mathbf{a}_m - \sum_{j=1}^{m-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_m$

Или поэлементно:

$\beta_{1j} = \alpha_{1j}$

$\beta_{ij} = \alpha_{ij} - (\sum_{k=1}^{i-1} \mathbf{proj}_{\mathbf{b}_k} \mathbf{a}_i)_j = \alpha_{ij} - \sum_{k=1}^{i-1}\frac{\sum_{s=1}^{n}\alpha_{is}\beta_{ks}}{\sum_{s=1}^{n}\beta_{ks}\beta_{ks}}\beta_{kj}$

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

Основное время работы алгоритма приходится на вычисление сумм $\sum_{i=1}^{j-1} \mathbf{proj}_{\mathbf{b}_i} \mathbf{a}_k \forall j: j=\overline{1, m}$.

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

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

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

Рассмотрим алгоритм, написанный в виде функции для математической системы MATLAB.

Входные данные.

Матрица $A$ с линейно-независимыми векторами в столбцах.

Выходные данные.

Матрица $Q$ с ортонормированным базисом $A$.

Использование.

Q = gramschmidt( A );

Реализация.

function Q = gramschmidt( A )
% Число векторов.
n = size( A, 2 );
 
% Инициализируем выходную матрицу.
Q = zeros( n );
 
% Преобразуем каждый вектор в базисный:
% (1) j-ый базисный вектор будет ортогонален каждому из предыдущих 1..j-1 векторов;
% (2) будет единичной длины.
for j = 1 : n
    % Выбираем j-ый вектор
    u = A( :, j );
    
    % Специальный случай для j = 1: просто нормируем этот вектор и помещаем
    % первым в найденный базис, так как нет векторов, относительно которых
    % нужно было бы делать его ортогональным.

    % Удаляем из оригинального вектора "u" все компоненты натянутые на базис
    % из векторов 1..j-1, их вклад будет удален.
    % => j-ый вектор ортогонален остальным.
    % Предыдущие базисные вектора были найдены на предыдуших шагах.
    % => соблюдаем принцип ортогональности, но не ортонормированности.
    for i = 1 : j - 1
        u = u - proj( Q(:,i), A(:,j) );
    end
     
    % Нормируем вектор
    Q(:,j) = u ./ norm( u );
end
end
 
% Проекция вектора "a" коллинеарно вектору "e".
function p = proj( e, a )
p = (e' * a) / (e' * e) .* e;
end

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

Рассмотрим $m$ векторов длины $n$.

1. Скалярное произведение векторов требует $n-1$ сложение и $n$ произведений.
2. Вычитание проекции вектора требует $2$ скалярных произведения, $1$ деление, $n$ умножений и $n$ сложений, то есть:
   1. $3n-2$ ($+$);
   2. $3n$ ($\times$);
   3. $1$ ($\div$).
3. Вычисление $i$-го вектора требует $i-1$ вычитаний проекций, то есть:
   1. $(3n-2)(i-1)$ ($+$);
   2. $3n(i-1)$ ($\times$);
   3. $i-1$ ($\div$).
4. Мы вычисляем вектора от $i=1$ до $m$, поэтому множители $(i-1)$ выражаются треугольным числом $(m-1)\frac{m}{2}$.
   1. $(3n-2)(m-1)\frac{m}{2}$ ($+$);
   2. $3n(m-1)\frac{m}{2}$ ($\times$);
   3. $(m-1)\frac{m}{2}$ ($\div$).

Таким образом, cложность алгоритма $O(nm^2)$.

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

Входные данные. Матрица $A$ с элементами $\alpha_{ij}$, где $i = \overline{1, m}$ (количество линейно независимых векторов) и $j = \overline{1, n}$ (размерность пространства), при этом $m \leqslant n$.

Объем входных данных. $mn$

Выходные данные. Матрица $B$ с элементами $\beta_{ij}$, где $i = \overline{1, m}$ и $j = \overline{1, n}$, тогда $b_{1}, \dots, b_{m}$ — система ортогональных векторов.

Объем выходных данных. $mn$

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

Отношение последовательной и параллельной сложности -- ?

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

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

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

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

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

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

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

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

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

Срок сдачи продлен до 15 ноября.

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

Данный алгоритм реализован во множестве библиотек и математических программных пакетах.

1. Библиотека Epicycle (C#).
2. Библиотека Spectral Python (SPy).
3. Библиотека vect (Haskell).

2.7.1 Пакет 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}}]

2.7.2 Пакет 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);

2.7.3 Пакет MATLAB

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

В данном случае для примера взяты векторы $\{-2,\;1,\;0\}$, $\{-2,\;0,\;1\}$, $\{-0{.}5,\;-1,\;1\}$.

[x y]=gschmidt([-2 1 0; -2 0 1; -0.5 -1 1])