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

Авторы описания алгоритма: Павлов Андрей, Филимонов Владимир.

Эта работа ждет рассмотрения преподавателем Дата последней правки страницы: 15.10.2016 Авторы этой статьи считают, что задание выполнено.

Ортогонализация Грамма-Шмидта [1]
Последовательный алгоритм
Последовательная сложность	$O(N^2M)$
Объём входных данных	$NM$
Объём выходных данных	$NM$

1 ЧАСТЬ. Свойства и структура алгоритмов

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

В конечномерном евклидовом пространстве существует ортонормированный базис.

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

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

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

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

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

Исходные данные: квадратная матрица A с линейно независимыми векторами $\mathbf{a}_1,...,\mathbf{a}_N | \mathbf{a}_i \in \R^{M}$.

Определяется оператор проекции $\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b$, где $\left \langle a,b \right \rangle$ - скалярное произведение векторов a и b. Данный оператор используется для проецирования вектора a коллинеарно вектору 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$

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

$\vdots$

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

На основе каждого вектора $\mathbf{b}_j (j = 1 \cdots N)$ может быть получен нормированный вектор $\mathbf{e}_j = \frac{\mathbf{b}_j}{||\mathbf{b}_j||}$.

Результаты процесса ортогонализации Грама-Шмидта: $\mathbf{b}_1\cdots\mathbf{b}_N$ - система ортогональных векторов либо система ортонормированных векторов $\mathbf{e}_1\cdots\mathbf{e}_N$.

Подробное описание можно найти в .

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

Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов $\sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N$ для каждого $b_i$.

Подробное описание общей сложности алгоритма представлено в разделе #Последовательная сложность алгоритма

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

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

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

Пусть дано множество векторов ${a_1, a_2, \cdots , a_N}$. Тогда:

Шаг 1: $\mathbf{b}_1 = \mathbf{a}_1$

Для вычисления остальных $\mathbf{b}_i$ следует последовательно выполнять Шаг 2 для всех $i$ таких что при $2 \le i \le N$.

Шаг 2: $\mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{j}} \mathbf{a}_i$, где для каждого $\mathbf{proj}_b a$ выполняется $\frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle}b$.

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

Оценим сложность последовательной реализации алгоритма.

Для вычисления одной операции проекции требуется:

${\color{Blue}\bullet} \quad M \cdot 3$ - умножений;

${\color{Blue}\bullet} \quad (M-1) \cdot 2$ - сложений;

${\color{Blue}\bullet} \quad 1$ - деление;

Общее число операций, необходимых для вычисления алгоритма следующее:

${\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot M \cdot 3$ - умножений;

${\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot (M-1) \cdot 3$ - сложений и вычитаний;

${\color{Blue}\bullet} \quad \frac{N(N-1)}{2}$ - делений.

Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма $O(N^2 \cdot M)$.

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

Опишем информационный граф алгоритма. В разделе #Схематичное описание информационного графа будет представлена схема с сокрытием подробностей реализации функции проекции. Такое ограничение позволяет продемонстрировать основные зависимости, не перегружая картинку лишними элементами. В разделе #Уточнение информационного графа будет раскрыта внутренняя структура операции проекции.

1.7.1 Схематичное описание информационного графа

Для вычисления $b_i$ требуется найти $proj_{b_{j}}a_{i}$ для всех $j \in [1, i]$. Следовательно для полного вычисления вектора $b_{i}$ требуется знать все $b_{j}$ с меньшим индексом. Такая зависимость по данным очень не удачна для параллелизма. Однако, если разбить процесс вычисления $b_{i}$ на несколько этапов, соответствующих функциям проекции ($proj_{b_{j}}a_{i}$), то это позволит производить некоторые предварительные вычисления для $b_{j}$ до момента, когда станут известны все предшествующие ей $b_{i}$.

На Рис. 1 изображена зависимость каждого из этапов от предыдущих вычислений. Первая строка овалов содержит все проекции, зависящие от $b_{1}$, вторая строка - все проекции, зависящие от $b_{2}$, $\;\cdots\;$, самая верхняя строка содержит проекцию, которая зависит от $b_{N-1}$.

Рис. 2 показывает зависимость по данным немного в другом формате. Каждая строка представляет собой набор данных, которые требуются для вычисления $b_{i}$ из первого столбца. Второй и последующие столбцы группируют проекции, зависящие от одной из $b_{i}$ и, с помощью стрелки, показывают эту зависимость.

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

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

1.7.2 Уточнение информационного графа

В этом разделе будет представлено уточнение информационной структуры графа изображенного на Рис. 1 с помощью раскрытия внутренней структуры операции проекции. Полученный граф несёт в себе новую информацию, которую можно использовать для улучшения алгоритма, но становится менее понятным.

На Рис. 3 изображен граф, аналогичный графу Рис. 1, но функция $proj$ изображена более подробно. Оранжевым овалом обозначаются скалярные произведения двух векторов. Голубой треугольник - принимает на вход два числа (являющиеся результатами скалярных произведений) и вычисляет их частное. Зелёным шестиугольником обозначается произведение вектора $a_i$ на число. Фиолетовый ромб складывает два числа, от которых имеется зависимость.

По сравнению с Рис. 1, более детальное описание информационной модели позволяет заметить повторяющиеся построчно скалярные произведения $(b_i, b_i)$. Это замечание позволяет уменьшить суммарное число операций, но увеличивает связность графа.

Рис. 3. Уточненный информационный граф алгоритма Грама-Шмидта

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

На основе графа зависимости по данным, представленном на Рис. 1 или Рис. 2 раздела #Информационный граф можно выделить следующий ресурс для параллелизации. После вычисления $b_{i}$ можно начинать параллельно вычислять все $proj_{b_i}a_{j}$. На Рис. 2 блоки, которые можно вычислять параллельно представлены столбцами.

Кроме того, при очень больших размерностях можно также распараллелить вычисление проекций. Так как в функции проекции участвует пара векторов $a$ и $b$, а сама операция проекции определяется как $\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b$, где $\left \langle a,b \right \rangle$ - скалярное произведение векторов $a$ и $b$, то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.

1.8.1 Оценка ускорения параллельной реализации алгоритма

Рассмотрим случай, когда присутствует возможность запуска бесконечного числа потоков. Попробуем оценить сложность вычисления ортогональных векторов, при условии, что параллельно будут запускаться все проекции использующие одинаковую $b_i$. В этом случае, останется критический путь зависимости по данным: $b_1 \leftarrow proj_{b_2}a_{2} \leftarrow b_2 \leftarrow proj_{b_3}a_{3} \leftarrow \ldots \leftarrow proj_{b_{N-1}}a_{N-1} \leftarrow b_N$. Следовательно для оценки сложность распараллеленной версии алгоритма достаточно почитать сложность представленной цепочки вычислений.

${\color{Blue}\bullet} \quad M \cdot 3$ умножения в одной проекции;

${\color{Blue}\bullet} \quad (M-1) \cdot 2$ сложения;

${\color{Blue}\bullet} \quad 1$ деление;

В цепочке зависимости данных присутствует $N-1$ проекция, значит всего:

${\color{Blue}\bullet} \quad N \cdot M \cdot 3$ - умножений;

${\color{Blue}\bullet} \quad N \cdot (M-1) \cdot 2$ - сложений;

${\color{Blue}\bullet} \quad N \cdot 1$ - делений.

Получаем итоговую сложность $O(N \cdot M)$.

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

${\color{Blue}\bullet} \quad 3 \cdot N$ умножений;

${\color{Blue}\bullet} \quad N$ делений;

${\color{Blue}\bullet} \quad N \cdot logM$ сложений;

Следовательно оценка сложности параллельного алгоритма станет равна $O(N \cdot logM)$.

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

Входные данные: прямоугольная матрица $A^{M \times N}$ (элементы $\mathbf{a}_{ij}$).

Ограничение: матрица $A$ должна быть невырожденной.

Объём входных данных: $N \cdot N$.

Выходные данные: матрица $B$ (система ортогональных векторов $\mathbf{b}_{1} \cdots \mathbf{b}_{N}$).

Объём выходных данных: $N \cdot N$.

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов : $\frac{N^2 \cdot M}{N \cdot log M} = \frac{N \cdot M}{log M}$.

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

Существуют несколько реализаций рассматриваемого алгоритма. Статья описывает три версии данного алгоритма: строковый классический алгоритм Г-Ш (Грама - Шмидта), столбцовый классический алгоритм Г-Ш и модифицированный алгоритм Г-Ш. Модифицированный алгоритм Грама-Шмидта является вычислительно более устойчивым , то есть позволяет уменьшить погрешность вычислений связанную с ошибкой выполнения операций с плавающей точкой.

Процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт $0$ (нулевой вектор) на шаге $j$, если $a_{j}$ является линейной комбинацией векторов $a_{1}, \ldots ,a_{j-1}$.

2 ЧАСТЬ. Программная реализация алгоритма

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

TODO: (до 15 ноября)

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

Пример реализации базового алгоритма на языке Си с использованием технологий OpenMP и MPI представлен в статье .

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

${\color{Blue}\bullet} \quad$ Mathematica - Orthogonalize[{v1,v2,…}] [2]

${\color{Blue}\bullet} \quad$ Maxima - gramschmidt(x) [3]

${\color{Blue}\bullet} \quad$ MatLab - orth(A) [4]

${\color{Blue}\bullet} \quad$ SciPy - scipy.linalg.orth(A) [5] - Ортогонализация с использованием алгоритма SVD.