Участник:Pandvik/Ортогонализация Грама - Шмидта
Авторы описания алгоритма: Павлов Андрей (разделы [1.7], [1.8], [2.7], [3] - поиск и изучение литературы), Филимонов Владимир (разделы [1.1]-[1.6], [1.9], [1.10], [2.7], [3] - поиск и изучение литературы).
| Ортогонализация Грамма-Шмидта [1] | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(N^2M) }[/math] |
| Объём входных данных | [math]\displaystyle{ NM }[/math] |
| Объём выходных данных | [math]\displaystyle{ NM }[/math] |
ЧАСТЬ. Свойства и структура алгоритмов
Общее описание алгоритма
В конечномерном евклидовом пространстве существует ортонормированный базис.
Для доказательства этого факта требуется находить и строить такие базисы. Построить ортонормированный базис можно, отталкиваясь от некоторого исходного базиса, при помощи алгоритма, который называют процессом ортогонализации Грама — Шмидта[1].
Процесс ортогонализации Грама-Шмидта используется для квадратных невырожденных матриц, которые преобразуются, либо уже преобразованы, к верхнему(нижнему) треугольному виду.
Процесс ортогонализации Грама-Шмидта нашёл применение в оптимизации оценивания параметров моделей управления объектом[2], в протоколах безопасности[3], в обработке сигналов[4], в вычислении локальных минимумов целочисленных решёток и многом другом.
Обычно, процесс ортогонализации используется как промежуточный шаг в других алгоритмах для уменьшения количества вычислений.
Математическое описание алгоритма
Исходные данные: квадратная матрица A с линейно независимыми векторами [math]\displaystyle{ \mathbf{a}_1,...,\mathbf{a}_N | \mathbf{a}_i \in \R^{M} }[/math].
Определяется оператор проекции [math]\displaystyle{ \mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b }[/math], где [math]\displaystyle{ \left \langle a,b \right \rangle }[/math] - скалярное произведение векторов a и b. Данный оператор используется для проецирования вектора a коллинеарно вектору b.
[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{ \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 }[/math]
[math]\displaystyle{ \vdots }[/math]
[math]\displaystyle{ \mathbf{b}_N = \mathbf{a}_N - \sum_{j=1}^{N-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N }[/math]
На основе каждого вектора [math]\displaystyle{ \mathbf{b}_j (j = 1 \cdots N) }[/math] может быть получен нормированный вектор [math]\displaystyle{ \mathbf{e}_j = \frac{\mathbf{b}_j}{||\mathbf{b}_j||} }[/math].
Результаты процесса ортогонализации Грама-Шмидта: [math]\displaystyle{ \mathbf{b}_1\cdots\mathbf{b}_N }[/math] - система ортогональных векторов либо система ортонормированных векторов [math]\displaystyle{ \mathbf{e}_1\cdots\mathbf{e}_N }[/math].
Подробное описание можно найти в [5] [6].
Вычислительное ядро алгоритма
Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов [math]\displaystyle{ \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N }[/math] для каждого [math]\displaystyle{ b_i }[/math].
Подробное описание общей сложности алгоритма представлено в разделе #Последовательная сложность алгоритма
Макроструктура алгоритма
Как можно заметить из #Математическое описание алгоритма, в алгоритме преобладают векторные вычисления. Ключевым этапом процесса ортогонализации Грама-Шмидта является вычисление большого числа проекций для каждого вектора, каждый из которых состоит из трёх операций сложения и одной операции деления. В свою очередь, можно сделать вывод, что основную часть представленного метода составляют вычисления проекций векторов и сумм этих проекций [math]\displaystyle{ \sum_{j=1}^{N-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N }[/math], где [math]\displaystyle{ \mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b }[/math].
Схема реализации последовательного алгоритма
Пусть дано множество векторов [math]\displaystyle{ {a_1, a_2, \cdots , a_N} }[/math]. Тогда:
Шаг 1: [math]\displaystyle{ \mathbf{b}_1 = \mathbf{a}_1 }[/math]
Для вычисления остальных [math]\displaystyle{ \mathbf{b}_i }[/math] следует последовательно выполнять Шаг 2 для всех [math]\displaystyle{ i }[/math] таких что при [math]\displaystyle{ 2 \le i \le N }[/math].
Шаг 2: [math]\displaystyle{ \mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{j}} \mathbf{a}_i }[/math], где для каждого [math]\displaystyle{ \mathbf{proj}_b a }[/math] выполняется [math]\displaystyle{ \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle}b }[/math].
Последовательная сложность алгоритма
Оценим сложность последовательной реализации алгоритма.
Для вычисления одной операции проекции требуется:
[math]\displaystyle{ {\color{Blue}\bullet} \quad M \cdot 3 }[/math] - умножений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad (M-1) \cdot 2 }[/math] - сложений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad 1 }[/math] - деление;
Общее число операций, необходимых для вычисления алгоритма следующее:
[math]\displaystyle{ {\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot M \cdot 3 }[/math] - умножений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot (M-1) \cdot 3 }[/math] - сложений и вычитаний;
[math]\displaystyle{ {\color{Blue}\bullet} \quad \frac{N(N-1)}{2} }[/math] - делений.
Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма [math]\displaystyle{ O(N^2 \cdot M) }[/math].
Информационный граф
Опишем информационный граф алгоритма. В разделе #Схематичное описание информационного графа будет представлена схема с сокрытием подробностей реализации функции проекции. Такое ограничение позволяет продемонстрировать основные зависимости, не перегружая картинку лишними элементами. В разделе #Уточнение информационного графа будет раскрыта внутренняя структура операции проекции.
Схематичное описание информационного графа
Для вычисления [math]\displaystyle{ b_i }[/math] требуется найти [math]\displaystyle{ proj_{b_{j}}a_{i} }[/math] для всех [math]\displaystyle{ j \in [1, i] }[/math]. Следовательно для полного вычисления вектора [math]\displaystyle{ b_{i} }[/math] требуется знать все [math]\displaystyle{ b_{j} }[/math] с меньшим индексом. Такая зависимость по данным очень не удачна для параллелизма. Однако, если разбить процесс вычисления [math]\displaystyle{ b_{i} }[/math] на несколько этапов, соответствующих функциям проекции ([math]\displaystyle{ proj_{b_{j}}a_{i} }[/math]), то это позволит производить некоторые предварительные вычисления для [math]\displaystyle{ b_{j} }[/math] до момента, когда станут известны все предшествующие ей [math]\displaystyle{ b_{i} }[/math].
Далее в этом разделе будут представлены несколько вариантов информационного графа. Общим моментом является выделение двойными фиолетовыми кругами данных, поступающих на вход, а двойными зелёными кругами - те данные, которые будут получены на выходе алгоритма. В предложенных графах, остриё стрелок указывает на элемент, которому требуются некоторые данные. То есть исходящие стрелки обозначают результат операции, а входящие - аргументы операции.
На Рис. 1 изображена зависимость каждого из этапов от предыдущих вычислений. Первая строка овалов содержит все проекции, зависящие от [math]\displaystyle{ b_{1} }[/math], вторая строка - все проекции, зависящие от [math]\displaystyle{ b_{2} }[/math], [math]\displaystyle{ \;\cdots\; }[/math], самая верхняя строка содержит проекцию, которая зависит от [math]\displaystyle{ b_{N-1} }[/math]. Широкие красные стрелки обозначают суммирование элементов с сохранением результата в переменную, находящуюся на острие стрелки. Следует уточнить, что в данном варианте графа указаны только зависимости по данным вычисляемым непосредственно в процессе работы алгоритма. Так, например, опущена зависимость проекций от векторов [math]\displaystyle{ a_{i} }[/math], что позволяет сделать рисунок более наглядным.
Рис. 2 показывает зависимость по данным немного в другом формате. Каждая строка представляет собой набор данных, которые требуются для вычисления [math]\displaystyle{ b_{i} }[/math] из первого столбца. Второй и последующие столбцы группируют проекции, зависящие от одной из [math]\displaystyle{ b_{i} }[/math] и, с помощью стрелки, показывают эту зависимость. В данном варианте схемы зависимости по данным опущены входные данные, так как они изначально доступны для каждой вычислительной единицы.
Уточнение информационного графа
В этом разделе будет представлено уточнение информационной структуры графа изображенного на Рис. 1 с помощью раскрытия внутренней структуры операции проекции. Полученный граф несёт в себе новую информацию, которую можно использовать для улучшения алгоритма, но становится менее понятным.
На Рис. 3 изображен граф, аналогичный графу Рис. 1, но функция [math]\displaystyle{ proj }[/math] изображена более подробно. Оранжевым овалом обозначаются скалярные произведения двух векторов. Голубой треугольник - принимает на вход два числа (являющиеся результатами скалярных произведений) и вычисляет их частное. Зелёным шестиугольником обозначается произведение вектора [math]\displaystyle{ a_i }[/math] на число. Фиолетовый ромб складывает два числа, от которых имеется зависимость.
По сравнению с Рис. 1, более детальное описание информационной модели позволяет заметить повторяющиеся построчно скалярные произведения [math]\displaystyle{ (b_i, b_i) }[/math]. Это замечание позволяет уменьшить суммарное число операций, но увеличивает связность графа.
Ресурс параллелизма алгоритма
На основе графа зависимости по данным, представленном на Рис. 1 или Рис. 2 раздела #Информационный граф можно выделить следующий ресурс для параллелизации. После вычисления [math]\displaystyle{ b_{i} }[/math] можно начинать параллельно вычислять все [math]\displaystyle{ proj_{b_i}a_{j} }[/math]. На Рис. 2 блоки, которые можно вычислять параллельно представлены столбцами.
Кроме того, при очень больших размерностях можно также распараллелить вычисление проекций. Так как в функции проекции участвует пара векторов [math]\displaystyle{ a }[/math] и [math]\displaystyle{ 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{ \left \langle a,b \right \rangle }[/math] - скалярное произведение векторов [math]\displaystyle{ a }[/math] и [math]\displaystyle{ b }[/math], то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.
Оценка ускорения параллельной реализации алгоритма
Рассмотрим случай, когда присутствует возможность запуска бесконечного числа потоков. Попробуем оценить сложность вычисления ортогональных векторов, при условии, что параллельно будут запускаться все проекции использующие одинаковую [math]\displaystyle{ b_i }[/math]. В этом случае, останется критический путь зависимости по данным: [math]\displaystyle{ 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 }[/math]. Следовательно для оценки сложность распараллеленной версии алгоритма достаточно почитать сложность представленной цепочки вычислений.
[math]\displaystyle{ {\color{Blue}\bullet} \quad M \cdot 3 }[/math] умножения в одной проекции;
[math]\displaystyle{ {\color{Blue}\bullet} \quad (M-1) \cdot 2 }[/math] сложения;
[math]\displaystyle{ {\color{Blue}\bullet} \quad 1 }[/math] деление;
В цепочке зависимости данных присутствует [math]\displaystyle{ N-1 }[/math] проекция, значит всего:
[math]\displaystyle{ {\color{Blue}\bullet} \quad N \cdot M \cdot 3 }[/math] - умножений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad N \cdot (M-1) \cdot 2 }[/math] - сложений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad N \cdot 1 }[/math] - делений.
Получаем итоговую сложность [math]\displaystyle{ O(N \cdot M) }[/math].
Если дополнительно рассматривать возможность параллелизации векторных операций (скалярного произведения и умножения вектора на число), критический путь можно сократить до следующих операций:
[math]\displaystyle{ {\color{Blue}\bullet} \quad 3 \cdot N }[/math] умножений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad N }[/math] делений;
[math]\displaystyle{ {\color{Blue}\bullet} \quad N \cdot logM }[/math] сложений;
Следовательно оценка сложности параллельного алгоритма станет равна [math]\displaystyle{ O(N \cdot logM) }[/math].
Входные и выходные данные алгоритма
Входные данные: прямоугольная матрица [math]\displaystyle{ A^{M \times N} }[/math] (элементы [math]\displaystyle{ \mathbf{a}_{ij} }[/math]).
Ограничение: матрица [math]\displaystyle{ A }[/math] должна быть невырожденной.
Объём входных данных: [math]\displaystyle{ N \cdot N }[/math].
Выходные данные: матрица [math]\displaystyle{ B }[/math] (система ортогональных векторов [math]\displaystyle{ \mathbf{b}_{1} \cdots \mathbf{b}_{N} }[/math]).
Объём выходных данных: [math]\displaystyle{ N \cdot N }[/math].
Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов : [math]\displaystyle{ \frac{N^2 \cdot M}{N \cdot log M} = \frac{N \cdot M}{log M} }[/math].
Отношение числа операций к суммарному объему входных и выходных данных - линейное.
Существуют несколько реализаций рассматриваемого алгоритма. Статья [7] описывает три версии данного алгоритма: строковый классический алгоритм Г-Ш (Грама - Шмидта), столбцовый классический алгоритм Г-Ш[8] и модифицированный алгоритм Г-Ш. Модифицированный алгоритм Грама-Шмидта[9] является вычислительно более устойчивым [10], то есть позволяет уменьшить погрешность вычислений связанную с ошибкой выполнения операций с плавающей точкой.
Процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт [math]\displaystyle{ 0 }[/math] (нулевой вектор) на шаге [math]\displaystyle{ j }[/math], если [math]\displaystyle{ a_{j} }[/math] является линейной комбинацией векторов [math]\displaystyle{ a_{1}, \ldots ,a_{j-1} }[/math].
ЧАСТЬ. Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
TODO: (до 15 ноября)
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Пример реализации базового алгоритма на языке Си с использованием технологий OpenMP и MPI представлен в статье [11].
Во всех библиотеках предназначенных для математических вычислений реализованы алгоритмы ортогонализации. В наиболее популярных:
[math]\displaystyle{ {\color{Blue}\bullet} \quad }[/math] Mathematica - Orthogonalize[{v1,v2,…}] [2]
[math]\displaystyle{ {\color{Blue}\bullet} \quad }[/math] Maxima - gramschmidt(x) [3]
[math]\displaystyle{ {\color{Blue}\bullet} \quad }[/math] MatLab - orth(A) [4]
[math]\displaystyle{ {\color{Blue}\bullet} \quad }[/math] SciPy - scipy.linalg.orth(A) [5] - Ортогонализация с использованием алгоритма SVD.
Литература
- ↑ А.А. Самарский, А.В. Гулин, Численные методы
- ↑ А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления
- ↑ А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности
- ↑ МГУ ВМК, Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов
- ↑ А.Н. Канатников, А.П.Крищенко, Линейная Алгебра
- ↑ The University of Southern Mississippi, Department of Mathematics, Jim Lambers, Gram-Schmidt Orthogonalization
- ↑ Å. Björck, Numerics of Gram-Schmidt orthogonalization.
- ↑ Massachusetts Institute of Technology, Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization
- ↑ С.А. Белов, Н.Ю. Золотых, Численные Методы Линейной Алгебры
- ↑ Massachusetts Institute of Technology, Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization
- ↑ Tirana University, Department of Mathematics, Genci Berati, Gram–Schmidt Process in Different Parallel Platforms (Control Flow versus Data Flow)