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

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

Основные авторы описания: Белов Н. А. (пункты 1.3, 1.4, 1.5, 1.6, 1.9, 1.10, 2.4, 2.7, 3), Богомолов Д. А. (пункты 1.1, 1.2, 1.6, 1.7, 1.8, 2.4, 3).

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

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

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

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

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

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

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

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

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

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

[math]\mathbf{b}_1 = \mathbf{a}_1[/math]

[math]\mathbf{b}_2 = \mathbf{a}_2 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_2 [/math]

[math]\mathbf{b}_3 = \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_3 - \mathbf{proj}_{\mathbf{b}_2} \mathbf{a}_3 [/math]

[math]\vdots[/math]

[math]\mathbf{b}_m = \mathbf{a}_m - \sum_{j=1}^{m-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_m[/math]

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

[math]\beta_{1j} = \alpha_{1j}[/math]

[math]\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} [/math]

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

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

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

Из математического описания алгоритма и описания ядра алгоритма в предыдущих разделах следует, что основную часть данного метода составляют операции вычисления проекций векторов, а также их сумм. Таким образом, в алгоритме можно выделить следующие макрооперации: скалярное произведение векторов в количестве [math](m-1)m[/math], а также множественные (всего [math]m-1[/math]) вычисления сумм [math]\mathbf{a}_i - \sum_{j=1}^{i-1} \frac{\left \langle a_i,b_j \right \rangle}{\left \langle b_j,b_j \right \rangle }b_j[/math].

Макроструктура может быть представлена в следующем виде:

1. Вычисление проекций [math]\mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_i \forall j: j=\overline{1, i - 1}[/math] для текущего шага [math]i[/math].
2. Вычисление суммы [math]\mathbf{a}_i - \sum_{j=1}^{i-1} \frac{\left \langle a_i,b_j \right \rangle}{\left \langle b_j,b_j \right \rangle }b_j[/math] для текущего шага [math]i[/math].
3. Переход на следующий шаг.

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

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

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

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

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

Матрица [math]Q[/math] с ортонормированным базисом [math]A[/math].

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

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 Последовательная сложность алгоритма

Рассмотрим [math]m[/math] векторов длины [math]n[/math].

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

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

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

Граф состоит из элементов 2 видов: прямоугольников и овалов. Прямоугольниками обозначаются данные, причем красному цвету соответствуют входные данные, а зеленому ― выходные. Овалы означают операции над переданными данными. Голубым цветом обозначены более сложно организованные операции, а синим ― менее. Для нахождения [math]\mathbf{b}_i[/math] вектора требуются значения всех векторов с номерами, меньшими [math]i[/math], однако, расчеты проекций можно начинать, не дожидаясь вычисления всех векторов, а производить сразу, как только соответствующий ортогональный вектор будет рассчитан. Эта зависимость и отражена на графе ниже.

Рисунок 1. Информационный граф алгоритма Грама-Шмидта.

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

Для построения базиса методом Грама–Шмидта в параллельном случае, имея в распоряжении [math]m[/math] и более потоков, необходимо выполнить следующие ярусы:

1. [math]m-1[/math] ярус вычисления проекций;
2. [math]m-1[/math] ярус вычитания векторов.

Эти ярусы чередуются. Количество операций проекции и вычитания векторов уменьшается на единицу в каждом следующем ярусе. Вычисление проекции вектора требует 2 скалярных произведения, 1 деление и [math]n[/math] произведений, то есть:

1. [math]2n-2[/math] ([math]+[/math]);
2. [math]3n[/math] ([math]\times[/math]);
3. [math]1[/math] ([math]\div[/math]).

Если считать, что вычисление проекции происходит последовательно, то сложность по высоте [math]O(n)[/math], по ширине [math]O(1)[/math] (то есть одна операция проекции в ярусе представляется набором ярусов с простыми последовательными операциями). В самом первом ярусе находится [math]m-1[/math] проекций и это наибольшее количество операций проекций, тогда сложность по высоте будет [math]O(n)[/math], по ширине ― [math]O(m)[/math].

Вычисление вычитания векторов требует [math]n[/math] вычитаний, и если считать, что операция происходит последовательно, то сложность по высоте [math]O(n)[/math], по ширине [math]O(1)[/math]. Во втором ярусе находится [math]m-1[/math] вычитаний и это наибольшее количество операций вычитаний, тогда сложность по высоте будет [math]O(n)[/math], по ширине ― [math]O(m)[/math].

Из вышеизложенного следует, что при классификации по ширине ЯПФ алгоритм имеет сложность [math]O(m)[/math], а по высоте ― [math]O(nm)[/math].

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

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

Объем входных данных. [math]mn[/math]

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

Объем выходных данных. [math]mn[/math]

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

1. Отношение последовательной и параллельной сложности в предположении доступности неограниченного числа необходимых ресурсов — [math]m[/math]. При этом необходимо отметить, что отношение числа операций к суммарному объему входных и выходных данных (вычислительная мощность алгоритма) — линейное.
2. Алгоритм почти полностью детерминирован, то есть гарантирована единственность результата выполнения. Это означает, что возможно накопление ошибок округления, то есть векторы выходной матрицы часто не точно ортогональны. Из-за потери ортогональности в процессе вычислений классический процесс Грама-Шмидта называют численно неустойчивым. Однако, процесс Грама-Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Такой алгоритм называется модифицированным процессом Грама-Шмидта. Рассмотрим его отличие подробнее.
   [math] \begin{array}{l} \mathbf{b}_j\\ \\ \mathbf{b}_j \end{array} \begin{array}{c} = \\ \\ = \end{array} \underbrace{ \underbrace{ \underbrace{ \begin{array}{c} \mathbf{a}_j \\ \\ \mathbf{a}_j \end{array} \begin{array}{c} - \\ \\ - \end{array} \begin{array}{c} \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j \\ || \\ \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j \end{array} }_{\mathbf{a}_j^{(1)}} \begin{array}{c} - \\ \\ - \end{array} \begin{array}{c} \mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_j \\ || \\ \mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_j^{(1)} \end{array} }_{\mathbf{a}_j^{(2)}} \begin{array}{c} - \\ \\ - \end{array} \begin{array}{c} \mathbf{proj}_{\mathbf{b}_3}\,\mathbf{a}_j \\ || \\ \mathbf{proj}_{\mathbf{b}_3}\,\mathbf{a}_j^{(2)} \end{array} }_{\mathbf{a}_j^{(3)}} \begin{array}{ccc} - & \ldots & - \\ \\ - & \ldots & - \end{array} \begin{array}{cc} \mathbf{proj}_{\mathbf{b}_{j-1}}\,\mathbf{a}_j & (1) \\ || & \\ \mathbf{proj}_{\mathbf{b}_{j-1}}\,\mathbf{a}_j^{(j-2)} & (2) \end{array} [/math]
   Формула (1) показывает вычисление [math]\mathbf{b}_j[/math] в классическом процессе, а формула (2) — в модифицированном.
   Разница между ними заключается в том, от каких векторов вычисляются компоненты: от [math]\mathbf{a}_j[/math] в классическом процессе или от результата предыдущего вычитания, то есть от [math]\mathbf{a}_j^{(k)}[/math] в модифицированном процессе. Таким образом, [math]\mathbf{a}_j^{(k)}[/math] представляет собой [math]\mathbf{a}_j[/math], из которого удалены компоненты в направлениях [math]\mathbf{b}_1,\,\;\ldots,\;\mathbf{b}_{k}[/math]. Компонента вектора [math]\mathbf{a}_j[/math] в направлении [math]\mathbf{b}_{k+1}[/math] при этом удалении не затрагивается и поэтому она в [math]\mathbf{a}_j^{(k)}[/math] такая же, как в [math]\mathbf{a}_j[/math].
3. Процесс Грама-Шмидта может применяться к бесконечной последовательности линейно независимых векторов, а также к линейно зависимым векторам. В этом случае он выдаёт [math]\mathbf{0}[/math] (нулевой вектор) на шаге [math]j[/math], если [math]\mathbf{a}_j[/math] является линейной комбинацией векторов [math]\mathbf{a}_1,\;\ldots,\;\mathbf{a}_{j-1}[/math]. Если это может случиться, то для сохранения ортогональности выходных векторов и для предотвращения деления на ноль при ортонормировании алгоритм должен делать проверку на нулевые векторы и отбрасывать их. Количество векторов, выдаваемых алгоритмом, будет равно размерности подпространства, порождённого векторами (то есть количеству линейно независимых векторов, которые можно выделить среди исходных векторов).
4. Процесс Грама-Шмидта может быть истолкован как разложение невырожденной квадратной матрицы в произведение ортогональной (или унитарной матрицы в случае эрмитова пространства) и верхнетреугольной матрицы с положительными диагональными элементами — QR-разложение.

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

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

В рамках статьи была написана своя параллельная реализация процесса ортогонализации Грама-Шмидта с использованием OpenMPI. Реализация вместе с данными и вспомогательными скриптами доступна в репозитории GramSchmidt.

Испытания реализованного алгоритма проводились на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета со следующими значениями параметров:

• размер матрицы по одному измерению: 100, 500 и далее с шагом 500 до 10000 элементов,
• число процессов: 1, 2, 4, 8 и далее с шагом 8 до 128. Тестирование с использованием 1, 2, 4 и 8 процессов не использовалось для данных, где размер превышал 5000.

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

Рисунок 2. График производительности реализации алгоритма в зависимости от числа процессоров и размерности матриц.

Размер матрицы:	100	500	1000	1500	2000	2500	3000	3500	4000	4500	5000	5500	6000	6500	7000	7500	8000	8500	9000	9500	10000
Число узлов, на которых время работы данной реализации является минимальным:	100	500	1000	1500	2000	2500	3000	3500	4000	4500	5000	5500	6000	6500	7000	7500	8000	8500	72	80	80

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

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

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

2.7.1 Пакет Mathematica

Приведем пример скрипта, предназначенного для пакета Mathematica, который проводит процесс ортогонализации Грама-Шмидта над векторами, заданными в фигурных скобках последней строки. Количество векторов и их координат могут быть произвольными. В данном случае для примера взяты векторы [math]\{-2,\;1,\;0\}[/math], [math]\{-2,\;0,\;1\}[/math], [math]\{-0{.}5,\;-1,\;1\}[/math].

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, который проводит процесс ортогонализации Грама-Шмидта над векторами, заданными в фигурных скобках последней строки. Количество векторов и их координат могут быть произвольными. В данном случае для примера взяты векторы [math]\{-2,\;1,\;0\}[/math], [math]\{-2,\;0,\;1\}[/math], [math]\{-0{.}5,\;-1,\;1\}[/math].

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

2.7.3 Пакет MATLAB

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

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

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