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

Файл:Symbol confirmed.svg

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и ASA.

Ортогонализация Грама-Шмидта
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(nm^2) }[/math]
Объём входных данных	[math]\displaystyle{ nm }[/math]
Объём выходных данных	[math]\displaystyle{ nm }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(nm) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ 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).

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \vdots }[/math]

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

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

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

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

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

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

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

Из математического описания алгоритма и описания ядра алгоритма в предыдущих разделах следует, что основную часть данного метода составляют операции вычисления проекций векторов, а также их сумм. Таким образом, в алгоритме можно выделить следующие макрооперации: скалярное произведение векторов в количестве [math]\displaystyle{ (m-1)m }[/math], а также множественные (всего [math]\displaystyle{ m-1 }[/math]) вычисления сумм [math]\displaystyle{ \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]\displaystyle{ \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_i \forall j: j=\overline{1, i - 1} }[/math] для текущего шага [math]\displaystyle{ i }[/math].
2. Вычисление суммы [math]\displaystyle{ \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]\displaystyle{ i }[/math].
3. Переход на следующий шаг.

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

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

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

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

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

Матрица [math]\displaystyle{ Q }[/math] с ортонормированным базисом [math]\displaystyle{ 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1. Отношение последовательной и параллельной сложности в предположении доступности неограниченного числа необходимых ресурсов — [math]\displaystyle{ m }[/math]. При этом необходимо отметить, что отношение числа операций к суммарному объему входных и выходных данных (вычислительная мощность алгоритма) — линейное.
2. Алгоритм почти полностью детерминирован, то есть гарантирована единственность результата выполнения. Это означает, что возможно накопление ошибок округления, то есть векторы выходной матрицы часто не точно ортогональны. Из-за потери ортогональности в процессе вычислений классический процесс Грама-Шмидта называют численно неустойчивым. Однако, процесс Грама-Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Такой алгоритм называется модифицированным процессом Грама-Шмидта. Рассмотрим его отличие подробнее.
   [math]\displaystyle{ \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]\displaystyle{ \mathbf{b}_j }[/math] в классическом процессе, а формула (2) — в модифицированном.
   Разница между ними заключается в том, от каких векторов вычисляются компоненты: от [math]\displaystyle{ \mathbf{a}_j }[/math] в классическом процессе или от результата предыдущего вычитания, то есть от [math]\displaystyle{ \mathbf{a}_j^{(k)} }[/math] в модифицированном процессе. Таким образом, [math]\displaystyle{ \mathbf{a}_j^{(k)} }[/math] представляет собой [math]\displaystyle{ \mathbf{a}_j }[/math], из которого удалены компоненты в направлениях [math]\displaystyle{ \mathbf{b}_1,\,\;\ldots,\;\mathbf{b}_{k} }[/math]. Компонента вектора [math]\displaystyle{ \mathbf{a}_j }[/math] в направлении [math]\displaystyle{ \mathbf{b}_{k+1} }[/math] при этом удалении не затрагивается и поэтому она в [math]\displaystyle{ \mathbf{a}_j^{(k)} }[/math] такая же, как в [math]\displaystyle{ \mathbf{a}_j }[/math].
3. Процесс Грама-Шмидта может применяться к бесконечной последовательности линейно независимых векторов, а также к линейно зависимым векторам. В этом случае он выдаёт [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]. Если это может случиться, то для сохранения ортогональности выходных векторов и для предотвращения деления на ноль при ортонормировании алгоритм должен делать проверку на нулевые векторы и отбрасывать их. Количество векторов, выдаваемых алгоритмом, будет равно размерности подпространства, порождённого векторами (то есть количеству линейно независимых векторов, которые можно выделить среди исходных векторов).
4. Процесс Грама-Шмидта может быть истолкован как разложение невырожденной квадратной матрицы в произведение ортогональной (или унитарной матрицы в случае эрмитова пространства) и верхнетреугольной матрицы с положительными диагональными элементами — QR-разложение.

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

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

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

Основная идея реализации: нулевой вычислительный узел является управляющим, остальные ― подчиненными. Подчиненные узлы запрашивают "задание" у управляющего узла ― номер вектора, который требуется расчитать. Когда подчиненный узел заканчивает расчет, то пересылает готовый вектор на управляющий узел и запрашивает следующее задание. Если заданий больше нет, то подчиненный узел заканчивает свою работу. Каждый узел имеет список уже готовых векторов. Если при расчете требуется вектор, которого нет в списке, то он запрашивается у управляющего узла.

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

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

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

Использованные модули и компиляторы
module add slurm module add openmpi/1.8.4-gcc module add gcc/5.2.0 Для самой программы был использован компилятор mpicxx, для генератора и тестера ― g++.

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

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

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

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

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

Рисунок 4. График эффективности распараллеливания алгоритма в зависимости от числа процессоров и размерности матриц.

Минимальная эффективность распараллеливания 0.12%. Максимальная ― 95%.

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

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

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

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

Пакет Mathematica

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

Пакет Maxima

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

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

Пакет MATLAB

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

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

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