Участник:Pandvik/Ортогонализация Грама - Шмидта
![]() | Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и ASA. |
Авторы описания алгоритма: Павлов Андрей (разделы [1.7], [1.8], [2.4], [2.7], [3] - поиск и изучение литературы), Филимонов Владимир (разделы [1.1]-[1.6], [1.9], [1.10], [2.7], [3] - поиск и изучение литературы).
Ортогонализация Грамма-Шмидта [1] | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(N^2M)[/math] |
Объём входных данных | [math]NM[/math] |
Объём выходных данных | [math]NM[/math] |
Содержание
- 1 ЧАСТЬ. Свойства и структура алгоритмов
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 ЧАСТЬ. Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 ЧАСТЬ. Свойства и структура алгоритмов
1.1 Общее описание алгоритма
В конечномерном евклидовом пространстве существует ортонормированный базис.
Для доказательства этого факта требуется находить и строить такие базисы. Построить ортонормированный базис можно, отталкиваясь от некоторого исходного базиса, при помощи алгоритма, который называют процессом ортогонализации Грама — Шмидта[1].
Процесс ортогонализации Грама-Шмидта используется для квадратных невырожденных матриц, которые преобразуются, либо уже преобразованы, к верхнему(нижнему) треугольному виду.
Процесс ортогонализации Грама-Шмидта нашёл применение в оптимизации оценивания параметров моделей управления объектом[2], в протоколах безопасности[3], в обработке сигналов[4], в вычислении локальных минимумов целочисленных решёток и многом другом.
Обычно, процесс ортогонализации используется как промежуточный шаг в других алгоритмах для уменьшения количества вычислений.
1.2 Математическое описание алгоритма
Исходные данные: квадратная матрица A с линейно независимыми векторами [math]\mathbf{a}_1,...,\mathbf{a}_N | \mathbf{a}_i \in \R^{M}[/math].
Определяется оператор проекции [math]\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b[/math], где [math]\left \langle a,b \right \rangle[/math] - скалярное произведение векторов a и b. Данный оператор используется для проецирования вектора a коллинеарно вектору b.
[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]\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]\vdots[/math]
[math]\mathbf{b}_N = \mathbf{a}_N - \sum_{j=1}^{N-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N[/math]
На основе каждого вектора [math]\mathbf{b}_j (j = 1 \cdots N)[/math] может быть получен нормированный вектор [math]\mathbf{e}_j = \frac{\mathbf{b}_j}{||\mathbf{b}_j||}[/math].
Результаты процесса ортогонализации Грама-Шмидта: [math]\mathbf{b}_1\cdots\mathbf{b}_N[/math] - система ортогональных векторов либо система ортонормированных векторов [math]\mathbf{e}_1\cdots\mathbf{e}_N[/math].
Подробное описание можно найти в [5] [6].
1.3 Вычислительное ядро алгоритма
Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов [math]\sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N[/math] для каждого [math]b_i[/math].
Подробное описание общей сложности алгоритма представлено в разделе #Последовательная сложность алгоритма
1.4 Макроструктура алгоритма
1. первый ортогональный вектор [math]\mathbf{b}_1[/math] задаётся с помощью координат вектора [math]\mathbf{a}_1[/math];
2. второй ортогональный вектор получаем с помощью вычитания из [math]\mathbf{a}_2[/math] проекции полученного ортогонального вектора [math]\mathbf{b}_1[/math] на [math]\mathbf{a}_2[/math];
[math]\vdots[/math]
N. каждый новый ортогональный вектор получаем с помощью вычитания из [math]\mathbf{a}_N[/math] суммы проекций полученных предыдущих векторов [math]\mathbf{b}_1 \cdots \mathbf{b}_{N-1}[/math] на [math]\mathbf{a}_N[/math], где [math]\mathbf{a}_1 \cdots \mathbf{a}_{N}[/math] - векторы исходной матрицы [math]A[/math], а [math]\mathbf{b}_1 \cdots \mathbf{b}_{N}[/math] - получаемые ортогональные векторы матрицы [math]B[/math].
Как можно заметить из #Математическое описание алгоритма, в алгоритме преобладают векторные вычисления. Ключевым этапом процесса ортогонализации Грама-Шмидта является вычисление большого числа проекций для каждого вектора, каждый из которых состоит из трёх операций сложения и одной операции деления. В свою очередь, можно сделать вывод, что основную часть представленного метода составляют вычисления проекций векторов и сумм этих проекций [math]\sum_{j=1}^{N-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N[/math], где [math]\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b[/math].
1.5 Схема реализации последовательного алгоритма
Пусть дано множество векторов [math]{a_1, a_2, \cdots , a_N}[/math]. Тогда:
Шаг 1: [math]\mathbf{b}_1 = \mathbf{a}_1[/math]
Для вычисления остальных [math]\mathbf{b}_i [/math] следует последовательно выполнять Шаг 2 для всех [math] i [/math] таких что при [math] 2 \le i \le N [/math].
Шаг 2: [math]\mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{j}} \mathbf{a}_i[/math], где для каждого [math]\mathbf{proj}_b a[/math] выполняется [math]\frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle}b[/math].
1.6 Последовательная сложность алгоритма
Оценим сложность последовательной реализации алгоритма.
Для вычисления одной операции проекции требуется:
[math]{\color{Blue}\bullet} \quad M \cdot 3[/math] - умножений;
[math]{\color{Blue}\bullet} \quad (M-1) \cdot 2[/math] - сложений;
[math]{\color{Blue}\bullet} \quad 1 [/math] - деление;
Общее число операций, необходимых для вычисления алгоритма следующее:
[math]{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot M \cdot 3[/math] - умножений;
[math]{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot (M-1) \cdot 3 [/math] - сложений и вычитаний;
[math]{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} [/math] - делений.
Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма [math]O(N^2 \cdot M)[/math].
1.7 Информационный граф
Опишем информационный граф алгоритма. В разделе #Схематичное описание информационного графа будет представлена схема с сокрытием подробностей реализации функции проекции. Такое ограничение позволяет продемонстрировать основные зависимости, не перегружая картинку лишними элементами. В разделе #Уточнение информационного графа будет раскрыта внутренняя структура операции проекции.
1.7.1 Схематичное описание информационного графа
Для вычисления [math]b_i[/math] требуется найти [math]proj_{b_{j}}a_{i}[/math] для всех [math]j \in [1, i][/math]. Следовательно для полного вычисления вектора [math]b_{i}[/math] требуется знать все [math]b_{j}[/math] с меньшим индексом. Такая зависимость по данным очень не удачна для параллелизма. Однако, если разбить процесс вычисления [math]b_{i}[/math] на несколько этапов, соответствующих функциям проекции ([math]proj_{b_{j}}a_{i}[/math]), то это позволит производить некоторые предварительные вычисления для [math]b_{j}[/math] до момента, когда станут известны все предшествующие ей [math]b_{i}[/math].
Далее в этом разделе будут представлены несколько вариантов информационного графа. Общим моментом является выделение двойными фиолетовыми кругами данных, поступающих на вход, а двойными зелёными кругами - те данные, которые будут получены на выходе алгоритма. В предложенных графах, остриё стрелок указывает на элемент, которому требуются некоторые данные. То есть исходящие стрелки обозначают результат операции, а входящие - аргументы операции.
На Рис. 1 изображена зависимость каждого из этапов от предыдущих вычислений. Первая строка овалов содержит все проекции, зависящие от [math]b_{1}[/math], вторая строка - все проекции, зависящие от [math]b_{2}[/math], [math]\;\cdots\;[/math], самая верхняя строка содержит проекцию, которая зависит от [math]b_{N-1}[/math]. Широкие красные стрелки обозначают суммирование элементов с сохранением результата в переменную, находящуюся на острие стрелки. Следует уточнить, что в данном варианте графа указаны только зависимости по данным вычисляемым непосредственно в процессе работы алгоритма. Так, например, опущена зависимость проекций от векторов [math]a_{i}[/math], что позволяет сделать рисунок более наглядным.
Рис. 2 показывает зависимость по данным немного в другом формате. Каждая строка представляет собой набор данных, которые требуются для вычисления [math]b_{i}[/math] из первого столбца. Второй и последующие столбцы группируют проекции, зависящие от одной из [math]b_{i}[/math] и, с помощью стрелки, показывают эту зависимость. В данном варианте схемы зависимости по данным опущены входные данные, так как они изначально доступны для каждой вычислительной единицы.
1.7.2 Уточнение информационного графа
В этом разделе будет представлено уточнение информационной структуры графа изображенного на Рис. 1 с помощью раскрытия внутренней структуры операции проекции. Полученный граф несёт в себе новую информацию, которую можно использовать для улучшения алгоритма, но становится менее понятным.
На Рис. 3 изображен граф, аналогичный графу Рис. 1, но функция [math]proj[/math] изображена более подробно. Оранжевым овалом обозначаются скалярные произведения двух векторов. Голубой треугольник - принимает на вход два числа (являющиеся результатами скалярных произведений) и вычисляет их частное. Зелёным шестиугольником обозначается произведение вектора [math]a_i[/math] на число. Фиолетовый ромб складывает два числа, от которых имеется зависимость. Зелёными стрелками будут выделены потоки наиболее значимые для алгоритма потоки данных. А именно - распределение зависимостей от вычисленных векторов [math]b_i[/math]. Стрелки, показывающие зависимость от начальных векторов окрашены в серый цвет. Для всех стрелок ромбом обозначен исток стрелки, точкой обозначено ветвление стрелки.
По сравнению с Рис. 1, более детальное описание информационной модели позволяет заметить повторяющиеся построчно скалярные произведения [math](b_i, b_i)[/math]. Зависимость от такого скалярного произведения изображена красным цветом на Рис. 3. Это замечание позволяет уменьшить суммарное число операций, но увеличивает объём данных передаваемых между вычислительными элементами.
1.8 Ресурс параллелизма алгоритма
На основе графа зависимости по данным, представленном на Рис. 1 или Рис. 2 раздела #Информационный граф можно выделить следующий ресурс для параллелизации. После вычисления [math]b_{i}[/math] можно начинать параллельно вычислять все [math]proj_{b_i}a_{j}[/math]. На Рис. 2 блоки, которые можно вычислять параллельно представлены столбцами.
Кроме того, при очень больших размерностях можно также распараллелить вычисление проекций. Так как в функции проекции участвует пара векторов [math]a[/math] и [math]b[/math], а сама операция проекции определяется как [math]\mathbf{proj}_b a = \frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle }b[/math], где [math]\left \langle a,b \right \rangle[/math] - скалярное произведение векторов [math]a[/math] и [math]b[/math], то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.
1.8.1 Оценка ускорения параллельной реализации алгоритма
Рассмотрим случай, когда присутствует возможность запуска бесконечного числа потоков. Попробуем оценить сложность вычисления ортогональных векторов, при условии, что параллельно будут запускаться все проекции использующие одинаковую [math]b_i[/math]. В этом случае, останется критический путь зависимости по данным: [math]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]{\color{Blue}\bullet} \quad M \cdot 3[/math] умножения в одной проекции;
[math]{\color{Blue}\bullet} \quad (M-1) \cdot 2[/math] сложения/вычитаний;
[math]{\color{Blue}\bullet} \quad 1[/math] деление;
В цепочке зависимости данных присутствует [math]N-1[/math] проекция, значит всего:
[math]{\color{Blue}\bullet} \quad N \cdot M \cdot 3[/math] - умножений;
[math]{\color{Blue}\bullet} \quad N \cdot (M-1) \cdot 2[/math] - сложений/вычитаний;
[math]{\color{Blue}\bullet} \quad N \cdot 1[/math] - делений.
Получаем итоговую сложность [math]O(N \cdot M)[/math].
Если дополнительно рассматривать возможность параллелизации векторных операций (скалярного произведения и умножения вектора на число), критический путь можно сократить до следующих операций:
[math]{\color{Blue}\bullet} \quad 3 \cdot N [/math] умножений;
[math]{\color{Blue}\bullet} \quad N [/math] делений;
[math]{\color{Blue}\bullet} \quad N \cdot logM [/math] сложений;
Следовательно оценка сложности параллельного алгоритма станет равна [math]O(N \cdot logM)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: прямоугольная матрица [math]A^{M \times N}[/math] (элементы [math]\mathbf{a}_{ij}[/math]).
Ограничение: матрица [math]A[/math] должна быть невырожденной.
Объём входных данных: [math] N \cdot N[/math].
Выходные данные: матрица [math]B[/math] (система ортогональных векторов [math]\mathbf{b}_{1} \cdots \mathbf{b}_{N}[/math]).
Объём выходных данных: [math]N \cdot N[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов : [math]\frac{N^2 \cdot M}{N \cdot log M} = \frac{N \cdot M}{log M}[/math].
Отношение числа операций к суммарному объему входных и выходных данных - линейное.
Представленный алгоритм является численно неустойчивым в силу ошибок округления, которые накапливаются в суммах проекций, так как они могут привести к неортогональности всей системы. Чтобы исправить численную неустойчивость требуется для каждого искомого вектора [math]\mathbf{b}_j[/math] делать переортагонализацию вектора путём последовательного вычитания из [math]\mathbf{a}_j[/math] компенент по направлению [math]\mathbf{b}_j[/math]. В результате внесённых изменений из формулы [math]\mathbf{b}_j = \mathbf{a}_j - \sum_{i=1}^{j-1} \mathbf{proj}_{\mathbf{b}_i} \mathbf{a}_j[/math] получаем вычисления следующего вида:
[math]\mathbf{a}_j^{(1)}[/math] = [math]\mathbf{a}_j[/math] - [math]\mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_j[/math]
[math]\mathbf{a}_j^{(2)}[/math] = [math]\mathbf{a}_j^{(1)}[/math] - [math]\mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_j^{(1)}[/math]
[math]\vdots[/math]
[math]\mathbf{a}_j^{(j-2)}[/math] = [math]\mathbf{a}_j^{(j-3)}[/math] - [math]\mathbf{proj}_{\mathbf{b}_{j-2}} \mathbf{a}_j^{(j-3)}[/math]
[math]\mathbf{b}_j[/math] = [math]\mathbf{a}_j^{(j-2)}[/math] - [math]\mathbf{proj}_{\mathbf{b}_{j-1}} \mathbf{a}_j^{(j-2)}[/math] .
Существуют несколько реализаций рассматриваемого алгоритма. Статья [7] описывает три версии данного алгоритма: строковый классический алгоритм Г-Ш (Грама - Шмидта), столбцовый классический алгоритм Г-Ш[8] и модифицированный алгоритм Г-Ш. Также существует ортогонализация с использованием алгоритма SVD (примером является реализация данного метода в библиотеке SciPy - scipy.linalg.orth(A) [2]). Модифицированный алгоритм Грама-Шмидта[9] является вычислительно более устойчивым, то есть позволяет уменьшить погрешность вычислений связанную с ошибкой выполнения операций с плавающей точкой.
Процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт [math] 0 [/math] (нулевой вектор) на шаге [math] j [/math], если [math]a_{j}[/math] является линейной комбинацией векторов [math]a_{1}, \ldots ,a_{j-1}[/math].
2 ЧАСТЬ. Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
2.4.1 Масштабируемость алгоритма
Теоретические оценки предела масштабируемости алгоритма были рассмотрены в разделе #Ресурс параллелизма алгоритма.
2.4.2 Масштабируемость реализации алгоритма
Проведём исследование масштабируемости параллельной реализации согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов"[10] Суперкомпьютерного комплекса Московского университета.
Набор и границы значений изменяемых параметров запуска реализации алгоритма:
- число процессоров [8 : 64] с шагом 4;
- размер матрицы [2000 : 7000] с шагом 500-1000.
В результате проведённых экспериментов был получен следующий диапазон эффективности распараллеливания алгоритма:
- минимальная эффективность реализации 0.58%;
- максимальная эффективность реализации 9.48%.
На следующих рисунках приведены графики времени работы и эффективности выбранной реализации алгоритма ортогонализации Грама-Шмидта в зависимости от изменяемых параметров запуска.
развернутьТаблица зависимости времени работы реализации от размерности задачи и числа вычислительных узлов |
---|
На Рис. 4 изображен график зависимости времени решения задачи от размерности входных данных и числа вычислительных узлов. При фиксированном значении размерности задачи, график зависимости времени работы реализации от числа узлов имеет вид вогнутой функции (пример для размерности матрицы 5000 изображен на Рис.6). То есть при увеличении числа узлов время работы данной реализации начинает увеличиваться после некоторого порогового значения, зависящего от размерности решаемой задачи. Есть предположение, что этот эффект является следствием большого объема данных, который передаётся между вычислительными узлами.
Рис. 5 изображает эффективность распараллеливания данной реализации в зависимости от размерности задачи и числа вычислительных узлов. По этому графику можно сделать вывод, что максимальная эффективность (9.5%) была достигнута на задаче размерности 4500 на 8 и 16 вычислительных узлах. Так же наблюдается снижение эффективности распараллеливания при увеличении числа вычислительных узлов.
На основе результатов тестирования можно определить число вычислительных узлов, на которых время работы рассматриваемой реализации является минимальным:
Размерность задачи: | 2000 | 2500 | 3000 | 3500 | 4000 | 4500 | 5000 | 5500 | 6000 | 7000 |
---|---|---|---|---|---|---|---|---|---|---|
Число узлов, на которых время работы данной реализации является минимальным: | 16 | 16 | 24 | 24 | 32 | 32 | 40 | 40 | 40 | 56 |
2.4.3 Особенности тестирования реализации алгоритма
Для исследования возможностей параллелизма алгоритма ортогонализации Грама-Шмидта была написана собственная реализация [3] с использованием технологии MPI. В данной реализации задачи вычисления векторов [math]b_{j}[/math] равномерно распределяются по всем доступным вычислительным узлам. Как только на одном из узлов будет полностью вычислен некоторый [math]b_{j}[/math], происходит рассылка результата вычисления на все остальные узлы.
Тестирование данной реализации производилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса МГУ им. Ломоносова.
развернутьТехнические подробности |
---|
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Пример реализации базового алгоритма на языке Си с использованием технологий OpenMP и MPI представлен в статье [11].
Во всех библиотеках предназначенных для математических вычислений реализованы алгоритмы ортогонализации. В наиболее популярных:
[math]{\color{Blue}\bullet} \quad [/math] Mathematica - Orthogonalize[{v1,v2,…}] [5]
[math]{\color{Blue}\bullet} \quad [/math] Maxima - gramschmidt(x) [6]
[math]{\color{Blue}\bullet} \quad [/math] MatLab - orth(A) [7].
3 Литература
- ↑ А.А. Самарский, А.В. Гулин, Численные методы, Москва: «НАУКА», 1989
- ↑ А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления, Томск: Томский политехнический университет, 2006
- ↑ А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности, Санкт-Петербург: НИУ ИТМО, 2016
- ↑ Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов, Москва: МГУ ВМК, 2009
- ↑ А.Н. Канатников, А.П.Крищенко, Линейная Алгебра, Москва: МГТУ им. Н.Э. Баумана, 2002
- ↑ Jim Lambers, Gram-Schmidt Orthogonalization, The University of Southern Mississippi, Department of Mathematics, 2014
- ↑ Å. Björck, Numerics of Gram-Schmidt orthogonalization, Linkiiping University, Department of Mathematics, 1994
- ↑ Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization, Massachusetts Institute of Technology, 2006
- ↑ С.А. Белов, Н.Ю. Золотых, Численные Методы Линейной Алгебры, Нижний Новгород, 2008
- ↑ Воеводин Вл. В., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. В., Практика суперкомпьютера «Ломоносов», Москва: «Открытые системы», 2012
- ↑ Tirana University, Department of Mathematics, Genci Berati, Gram–Schmidt Process in Different Parallel Platforms (Control Flow versus Data Flow), International Journal of Advanced Research in Artificial Intelligence, 2015