Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 103 промежуточные версии 5 участников)
Строка 1: Строка 1:
−
Авторы описания алгоритма: Павлов Андрей, Филимонов Владимир.
+
{{Assignment|ASA|Frolov}}
 +
 
 +
Авторы описания алгоритма: [[Участник:Pandvik|Павлов Андрей]] (разделы [1.7], [1.8], [2.4], [2.7], [3] - поиск и изучение литературы),  [[Участник:Vvf63|Филимонов Владимир]] (разделы [1.1]-[1.6], [1.9], [1.10], [2.7], [3] - поиск и изучение литературы).
 +
{{algorithm
 +
| name              = Ортогонализация Грамма-Шмидта [https://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%BE%D1%86%D0%B5%D1%81%D1%81_%D0%93%D1%80%D0%B0%D0%BC%D0%B0_%E2%80%95_%D0%A8%D0%BC%D0%B8%D0%B4%D1%82%D0%B0]
 +
| serial_complexity = <math>O(N^2M)</math>
 +
| input_data        = <math>NM</math>
 +
| output_data      = <math>NM</math>
 +
}}
  
 
= ЧАСТЬ. Свойства и структура алгоритмов =
 
= ЧАСТЬ. Свойства и структура алгоритмов =
Строка 6: Строка 14:
 
В конечномерном евклидовом пространстве существует ортонормированный базис.
 
В конечномерном евклидовом пространстве существует ортонормированный базис.
  
−
Для доказательства этого факта требуется находить и строить такие базисы. Построить ортонормированный базис можно, отталкиваясь от некоторого исходного базиса, при помощи алгоритма, который называют процессом ортогонализации Грама — Шмидта.
+
Для доказательства этого факта требуется находить и строить такие базисы. Построить ортонормированный базис можно, отталкиваясь от некоторого исходного базиса, при помощи алгоритма, который называют процессом ортогонализации Грама — Шмидта<ref>А.А. Самарский, А.В. Гулин, Численные методы, Москва: «НАУКА», 1989</ref>.
  
 
Процесс ортогонализации Грама-Шмидта используется для квадратных невырожденных матриц, которые преобразуются, либо уже преобразованы, к верхнему(нижнему) треугольному виду.
 
Процесс ортогонализации Грама-Шмидта используется для квадратных невырожденных матриц, которые преобразуются, либо уже преобразованы, к верхнему(нижнему) треугольному виду.
  
−
Процесс ортогонализации Грама-Шмидта нашёл применение в оптимизации оценивания параметров моделей управления объектом, в протоколах безопасности, в обработке сигналов, в вычислении локальных минимумов целочисленных решёток и многом другом.
+
Процесс ортогонализации Грама-Шмидта нашёл применение в оптимизации оценивания параметров моделей управления объектом<ref>А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления, Томск: Томский политехнический университет, 2006</ref>, в протоколах безопасности<ref>А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности, Санкт-Петербург: НИУ ИТМО, 2016</ref>, в обработке сигналов<ref>Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов, Москва: МГУ ВМК, 2009</ref>, в вычислении локальных минимумов целочисленных решёток и многом другом.
  
 
Обычно, процесс ортогонализации используется как промежуточный шаг в других алгоритмах для уменьшения количества вычислений.
 
Обычно, процесс ортогонализации используется как промежуточный шаг в других алгоритмах для уменьшения количества вычислений.
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
−
Исходные данные: квадратная матрица A с линейно независимыми векторами <math>\mathbf{a}_1,...,\mathbf{a}_N</math>.
+
Исходные данные: квадратная матрица 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{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.
Строка 23: Строка 31:
 
<math>\mathbf{b}_2 = \mathbf{a}_2 - \mathbf{proj}_{\mathbf{b}_1} \mathbf{a}_2 </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}_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>\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>\vdots</math>
Строка 34: Строка 42:
  
 
Результаты процесса ортогонализации Грама-Шмидта: <math>\mathbf{b}_1\cdots\mathbf{b}_N</math> - система ортогональных векторов либо система ортонормированных векторов <math>\mathbf{e}_1\cdots\mathbf{e}_N</math>.
 
Результаты процесса ортогонализации Грама-Шмидта: <math>\mathbf{b}_1\cdots\mathbf{b}_N</math> - система ортогональных векторов либо система ортонормированных векторов <math>\mathbf{e}_1\cdots\mathbf{e}_N</math>.
 +
 +
Подробное описание можно найти в <ref>А.Н. Канатников, А.П.Крищенко, Линейная Алгебра, Москва: МГТУ им. Н.Э. Баумана, 2002</ref> <ref>Jim Lambers, Gram-Schmidt Orthogonalization, The University of Southern Mississippi, Department of Mathematics, 2014</ref>.
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
  
−
Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов <math>\sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N</math>, на основе которой получается <math>b_i</math>.
+
Основная вычислительная сложность данного алгоритма заключается в вычислении суммы проекций векторов <math>\sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_j} \mathbf{a}_N</math> для каждого <math>b_i</math>.
  
 
Подробное описание общей сложности алгоритма представлено в разделе [[#Последовательная сложность алгоритма]]
 
Подробное описание общей сложности алгоритма представлено в разделе [[#Последовательная сложность алгоритма]]
  
 
== Макроструктура алгоритма ==
 
== Макроструктура алгоритма ==
 +
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. <math>\mathbf{b}_1 = \mathbf{a}_1</math>
 
  
−
2. <math>\mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{j}} \mathbf{a}_i</math> при <math>i = 2 \cdots N</math>, где для каждого <math>\mathbf{proj}_b a</math> выполняется <math>\frac{\left \langle a,b \right \rangle}{\left \langle b,b \right \rangle}b</math>.
+
Пусть дано множество векторов <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>.
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
−
Последовательный алгоритм требует:
+
Оценим сложность последовательной реализации алгоритма.
 +
 
 +
Для вычисления одной операции проекции требуется:
 +
 
 +
<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}n</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}n</math> - делений;
+
<math>{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} \cdot (M-1) \cdot 3 </math> - сложений и вычитаний;
  
−
<math>{\color{Blue}\bullet} \quad n \cdot 2 \cdot (n-1)(n-2)</math> - умножений.
+
<math>{\color{Blue}\bullet} \quad \frac{N(N-1)}{2} </math> - делений.
  
−
Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма <math>O(n^3)</math>.
+
Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма <math>O(N^2 \cdot M)</math>.
  
 
== Информационный граф ==
 
== Информационный граф ==
−
Опишем информационный граф алгоритма.
+
 
 +
Опишем информационный граф алгоритма. В разделе [[#Схематичное описание информационного графа]] будет представлена схема с сокрытием подробностей реализации функции проекции. Такое ограничение позволяет продемонстрировать основные зависимости, не перегружая картинку лишними элементами. В разделе [[#Уточнение информационного графа]] будет раскрыта внутренняя структура операции проекции.
 +
 
 +
=== Схематичное описание информационного графа ===
  
 
Для вычисления <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>.
 
Для вычисления <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>.
+
Далее в этом разделе будут представлены несколько вариантов информационного графа. Общим моментом является выделение двойными фиолетовыми кругами данных, поступающих на вход, а двойными зелёными кругами - те данные, которые будут получены на выходе алгоритма. В предложенных графах, остриё стрелок указывает на элемент, которому требуются некоторые данные. То есть исходящие стрелки обозначают результат операции, а входящие - аргументы операции.
 +
 
 +
На Рис. 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> и, с помощью стрелки, показывают эту зависимость.  
+
Рис. 2 показывает зависимость по данным немного в другом формате. Каждая строка представляет собой набор данных, которые требуются для вычисления <math>b_{i}</math> из первого столбца. Второй и последующие столбцы группируют проекции, зависящие от одной из <math>b_{i}</math> и, с помощью стрелки, показывают эту зависимость. В данном варианте схемы зависимости по данным опущены входные данные, так как они изначально доступны для каждой вычислительной единицы.
  
−
[[Файл:Diagram_Gram-Schmidt_DataFlow.png|center|frame|Рис. 1. Информационный граф алгоритма Грама-Шмидта]]
+
[[Файл:Diagram_Gram-Schmidt_DataFlow.png|100px|center|frame|Рис. 1. Информационный граф алгоритма Грама-Шмидта]]
−
[[Файл:Diagram_Gram-Schmidt_DataFlow_1.png|center|frame|Рис. 2. Информационный граф алгоритма Грама-Шмидта]]
+
[[Файл:Diagram_Gram-Schmidt_DataFlow_1.png|100px|center|frame|Рис. 2. Информационный граф алгоритма Грама-Шмидта]]
  
 +
=== Уточнение информационного графа ===
 +
В этом разделе будет представлено уточнение информационной структуры графа изображенного на Рис. 1 с помощью раскрытия внутренней структуры операции проекции. Полученный граф несёт в себе новую информацию, которую можно использовать для улучшения алгоритма, но становится менее понятным.
  
 +
На Рис. 3 изображен граф, аналогичный графу Рис. 1, но функция <math>proj</math> изображена более подробно. Оранжевым овалом обозначаются скалярные произведения двух векторов. Голубой треугольник - принимает на вход два числа (являющиеся результатами скалярных произведений) и вычисляет их частное. Зелёным шестиугольником обозначается произведение вектора <math>a_i</math> на число. Фиолетовый ромб складывает два числа, от которых имеется зависимость. Зелёными стрелками будут выделены потоки наиболее значимые для алгоритма потоки данных. А именно - распределение зависимостей от вычисленных векторов <math>b_i</math>. Стрелки, показывающие зависимость от начальных векторов окрашены в серый цвет. Для всех стрелок ромбом обозначен исток стрелки, точкой обозначено ветвление стрелки.
 +
 +
По сравнению с Рис. 1, более детальное описание информационной модели позволяет заметить повторяющиеся построчно скалярные произведения <math>(b_i, b_i)</math>. Зависимость от такого скалярного произведения изображена красным цветом на Рис. 3. Это замечание позволяет уменьшить суммарное число операций, но увеличивает объём данных передаваемых между вычислительными элементами.
 +
 +
[[Файл:Diagram_Gram-Schmidt_DataFlow_proj.png|center|100px|frame|Рис. 3. Уточненный информационный граф алгоритма Грама-Шмидта]]
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
  
−
На основе графа зависимости по данным, представленном на Рис. 1 или Рис. 2 раздела [[#Информационный граф]] можно выделить следующий ресурс параллелизма. После вычисления <math>b_{i}</math> можно начинать параллельно вычислять все <math>proj_{b_i}a_{j}</math>. На Рис. 2 блоки, которые можно вычислять параллельно представлены столбцами.
+
На основе графа зависимости по данным, представленном на Рис. 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>, то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.
 
Кроме того, при очень больших размерностях можно также распараллелить вычисление проекций. Так как в функции проекции участвует пара векторов <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>, то с помощью использования векторных операций можно получить некоторый выигрыш в производительности.
Строка 83: Строка 130:
 
Рассмотрим случай, когда присутствует возможность запуска бесконечного числа потоков. Попробуем оценить сложность вычисления ортогональных векторов, при условии, что параллельно будут запускаться все проекции использующие одинаковую <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>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 N \cdot 3</math> умножения в одной проекции;
+
<math>{\color{Blue}\bullet} \quad M \cdot 3</math> умножения в одной проекции;
  
−
<math>{\color{Blue}\bullet} \quad (N-1) \cdot 2</math> сложения и 1 деление.
+
<math>{\color{Blue}\bullet} \quad (M-1) \cdot 2</math> сложения/вычитаний;
 +
 
 +
<math>{\color{Blue}\bullet} \quad 1</math> деление;
  
 
В цепочке зависимости данных присутствует <math>N-1</math> проекция, значит всего:
 
В цепочке зависимости данных присутствует <math>N-1</math> проекция, значит всего:
  
−
<math>{\color{Blue}\bullet} \quad N \cdot N \cdot 3</math> - умножений;
+
<math>{\color{Blue}\bullet} \quad N \cdot M \cdot 3</math> - умножений;
  
−
<math>{\color{Blue}\bullet} \quad N \cdot (N-1) \cdot 2</math> - сложений;
+
<math>{\color{Blue}\bullet} \quad N \cdot (M-1) \cdot 2</math> - сложений/вычитаний;
  
 
<math>{\color{Blue}\bullet} \quad N \cdot 1</math> - делений.
 
<math>{\color{Blue}\bullet} \quad N \cdot 1</math> - делений.
  
−
Получаем итоговую сложность <math>O(N^{2})</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>.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
−
Входные данные: квадратная матрица A (элементы <math>\mathbf{a}_{ij}</math>). Матрица может быть приведена к верхнему(нижнему) треугольному виду.
+
Входные данные: прямоугольная матрица <math>A^{M \times N}</math> (элементы <math>\mathbf{a}_{ij}</math>).
  
−
Ограничение: матрица A должна быть невырожденной.
+
Ограничение: матрица <math>A</math> должна быть невырожденной.
  
−
Объём входных данных: <math>\frac{n(n+1)}{2}</math> (после преобразования матрицы к треугольному виду достаточно хранить только ненулевые элементы).
+
Объём входных данных: <math> N \cdot N</math>.
  
−
Выходные данные: матрица B (система ортогональных векторов <math>\mathbf{b}_{1} \cdots \mathbf{b}_{N}</math>).
+
Выходные данные: матрица <math>B</math> (система ортогональных векторов <math>\mathbf{b}_{1} \cdots \mathbf{b}_{N}</math>).
  
−
Объём выходных данных: <math>\frac{n(n+1)}{2}</math> (в силу треугольного вида матрицы).
+
Объём выходных данных: <math>N \cdot N</math>.
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
−
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов является линейной(последовательная реализация - <math>O(n^3)</math>, параллельная реализация - <math>O(n^2)</math>).
+
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов : <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> .
 +
 
 +
Существуют несколько реализаций рассматриваемого алгоритма. Статья <ref>Å. Björck, Numerics of Gram-Schmidt orthogonalization, Linkiiping University, Department of Mathematics, 1994</ref> описывает три версии данного алгоритма: строковый классический алгоритм Г-Ш (Грама - Шмидта), столбцовый классический алгоритм Г-Ш<ref>Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization, Massachusetts Institute of Technology, 2006</ref> и модифицированный алгоритм Г-Ш. Также существует ортогонализация с использованием алгоритма SVD (примером является реализация данного метода в библиотеке SciPy - scipy.linalg.orth(A) [http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.orth.html]). Модифицированный алгоритм Грама-Шмидта<ref>С.А. Белов, Н.Ю. Золотых, Численные Методы Линейной Алгебры, Нижний Новгород, 2008</ref> является вычислительно более устойчивым, то есть позволяет уменьшить погрешность вычислений связанную с ошибкой выполнения операций с плавающей точкой.
 +
 
 +
Процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт <math> 0 </math>  (нулевой вектор) на шаге <math> j </math>, если <math>a_{j}</math> является линейной комбинацией векторов <math>a_{1}, \ldots ,a_{j-1}</math>.
  
 
= ЧАСТЬ. Программная реализация алгоритма =
 
= ЧАСТЬ. Программная реализация алгоритма =
Строка 124: Строка 197:
  
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
 +
==== Масштабируемость алгоритма ====
 +
Теоретические оценки предела масштабируемости алгоритма были рассмотрены в разделе [[#Ресурс параллелизма алгоритма]].
 +
==== Масштабируемость реализации алгоритма ====
 +
Проведём исследование масштабируемости параллельной реализации согласно [[Scalability methodology|методике]]. Исследование проводилось на суперкомпьютере "Ломоносов"<ref name="Lom">Воеводин Вл. В., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. В., Практика суперкомпьютера «Ломоносов», Москва: «Открытые системы», 2012</ref> [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 +
 +
Набор и границы значений изменяемых [[Глоссарий#Параметры запуска|параметров запуска]] реализации алгоритма:
 +
 +
* число процессоров [8 : 64] с шагом 4;
 +
* размер матрицы [2000 : 7000] с шагом 500-1000.
 +
 +
В результате проведённых экспериментов был получен следующий диапазон [[Глоссарий#Эффективность распараллеливания|эффективности распараллеливания]] алгоритма:
 +
 +
* минимальная эффективность реализации 0.58%;
 +
* максимальная эффективность реализации 9.48%.
 +
 +
На следующих рисунках приведены графики времени работы и [[Глоссарий#Эффективность распараллеливания|эффективности]] выбранной реализации алгоритма ортогонализации Грама-Шмидта в зависимости от изменяемых параметров запуска.
 +
[[Файл:Gram-Schmidt stat.png|center|100px|frame|Рис. 4. График времени работы реализации алгоритма Грама-Шмидта (с использованием технологии MPI)]]
 +
 +
[[Файл:Gram-Schmidt stat eff.png|center|100px|frame|Рис. 5. График эффективности распараллеливания реализации алгоритма Грама-Шмидта относительно последовательной реализации]]
 +
 +
{| class="wikitable mw-collapsible mw-collapsed"
 +
! Таблица зависимости времени работы реализации от размерности задачи и числа вычислительных узлов
 +
|-
 +
|
 +
''' Время(сек.) | N | M | CPUS '''
 +
  8.214269    2000 2000 8
 +
  5.535392    2000 2000 16
 +
  7.121008    2000 2000 24
 +
  10.995030  2000 2000 32
 +
  17.709424  2000 2000 40
 +
 +
  15.969897  2500 2500 8
 +
  9.438978    2500 2500 16
 +
  9.658329    2500 2500 24
 +
  12.944879  2500 2500 32
 +
  17.741019  2500 2500 40
 +
  25.414062  2500 2500 48
 +
 +
  27.463480  3000 3000 8
 +
  15.192174  3000 3000 16
 +
  13.546558  3000 3000 24
 +
  14.893935  3000 3000 32
 +
  19.726213  3000 3000 40
 +
  27.127187  3000 3000 48
 +
  35.446089  3000 3000 56
 +
 +
  43.528986  3500 3500 8
 +
  23.267714  3500 3500 16
 +
  18.832487  3500 3500 24
 +
  20.074890  3500 3500 32
 +
  22.985693  3500 3500 40
 +
  27.376960  3500 3500 48
 +
  37.711549  3500 3500 56
 +
  48.896100  3500 3500 64
 +
 +
  64.889743  4000 4000 8
 +
  34.014309  4000 4000 16
 +
  25.912094  4000 4000 24
 +
  24.668768  4000 4000 32
 +
  27.262719  4000 4000 40
 +
  32.295430  4000 4000 48
 +
  41.851132  4000 4000 56
 +
  48.027403  4000 4000 64
 +
 +
  92.075673  4500 4500 8
 +
  47.766134  4500 4500 16
 +
  35.229498  4500 4500 24
 +
  31.972291  4500 4500 32
 +
  33.277924  4500 4500 40
 +
  38.653318  4500 4500 48
 +
  45.888597  4500 4500 56
 +
  54.089960  4500 4500 64
 +
 +
  126.197428  5000 5000 8
 +
  64.872589  5000 5000 16
 +
  46.852009  5000 5000 24
 +
  40.500379  5000 5000 32
 +
  39.817172  5000 5000 40
 +
  42.630904  5000 5000 48
 +
  48.521401  5000 5000 56
 +
  55.565563  5000 5000 64
 +
 +
  85.984754  5500 5500 16
 +
  60.770627  5500 5500 24
 +
  50.246163  5500 5500 32
 +
  47.708869  5500 5500 40
 +
  49.754522  5500 5500 48
 +
  52.557597  5500 5500 56
 +
  59.170178  5500 5500 64
 +
 +
  111.130398  6000 6000 16
 +
  77.988458  6000 6000 24
 +
  63.733449  6000 6000 32
 +
  57.554799  6000 6000 40
 +
  57.613826  6000 6000 48
 +
  62.380190  6000 6000 56
 +
  68.328533  6000 6000 64
 +
 +
  121.137472  7000 7000 24
 +
  96.815082  7000 7000 32
 +
  84.331911  7000 7000 40
 +
  79.391282  7000 7000 48
 +
  78.346921  7000 7000 56
 +
  84.036373  7000 7000 64
 +
 +
 +
|}
 +
 +
На Рис. 4 изображен график зависимости времени решения задачи от размерности входных данных и числа вычислительных узлов. При фиксированном значении размерности задачи, график зависимости времени работы реализации от числа узлов имеет вид вогнутой функции (пример для размерности матрицы 5000 изображен на Рис.6). То есть при увеличении числа узлов время работы данной реализации начинает увеличиваться после некоторого порогового значения, зависящего от размерности решаемой задачи. Есть предположение, что этот эффект является следствием большого объема данных, который передаётся между вычислительными узлами.
 +
 +
Рис. 5 изображает эффективность распараллеливания данной реализации в зависимости от размерности задачи и числа вычислительных узлов. По этому графику можно сделать вывод, что максимальная эффективность (9.5%) была достигнута на задаче размерности 4500 на 8 и 16 вычислительных узлах. Так же наблюдается снижение эффективности распараллеливания при увеличении числа вычислительных узлов.
 +
 +
На основе результатов тестирования можно определить число вычислительных узлов, на которых время работы рассматриваемой реализации является минимальным:
 +
 +
{| class="wikitable"
 +
|-
 +
! Размерность задачи:
 +
! 2000
 +
! 2500
 +
! 3000
 +
! 3500
 +
! 4000
 +
! 4500
 +
! 5000
 +
! 5500
 +
! 6000
 +
! 7000
 +
|-
 +
| Число узлов, на которых время работы данной реализации является минимальным:
 +
| 16
 +
| 16
 +
| 24
 +
| 24
 +
| 32
 +
| 32
 +
| 40
 +
| 40
 +
| 40
 +
| 56
 +
|}
 +
 +
[[Файл:Gram-Schmidt stat5000.png|center|100px|frame|Рис. 6. График времени работы параллельной реализации алгоритма Грама-Шмидта при размерности исходной матрицы 5000x5000]]
 +
 +
==== Особенности тестирования реализации алгоритма ====
 +
 +
Для исследования возможностей параллелизма алгоритма ортогонализации Грама-Шмидта была написана собственная реализация [http://bitbucket.org/pandvik/gram_schmidt_process] с использованием технологии MPI. В данной реализации задачи вычисления векторов <math>b_{j}</math> равномерно распределяются по всем доступным вычислительным узлам. Как только на одном из узлов будет полностью вычислен некоторый <math>b_{j}</math>, происходит рассылка результата вычисления на все остальные узлы.
 +
 +
Тестирование данной реализации производилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster/lomonosov.html Суперкомпьютерного комплекса МГУ им. Ломоносова].
 +
 +
{| class="wikitable mw-collapsible mw-collapsed"
 +
! Технические подробности
 +
|-
 +
|
 +
Компилятор: openmpi/1.8.4-gcc
 +
 +
Планировщик задач: slurm/2.5.6
 +
 +
Очереди: test, regular4
 +
 +
Команда постановки скрипта в очередь: sbatch -n <CPUS> --time=20 ompi a.out
 +
 +
Тестовые данные (матрицы различной размерности) были сгенерированы с помощью скрипта Gram_Schmidt_process/scripts/matrixGenerator.py из репозитория [http://bitbucket.org/pandvik/gram_schmidt_process]. Генерировались случайные квадратные матрицы выбранной размерности.
 +
|}
  
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==
Строка 130: Строка 366:
  
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
−
Пример реализации базового алгоритма с использованием технологий OpenMP и MPI представлен в статье <ref>Tirana University, Department of Mathematics, Genci Berati, Gram–Schmidt Process in Different Parallel Platforms (Control Flow versus Data Flow)</ref>.
+
Пример реализации базового алгоритма на языке Си с использованием технологий OpenMP и MPI представлен в статье <ref>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</ref>.
 +
 
 +
Во всех библиотеках предназначенных для математических вычислений реализованы алгоритмы ортогонализации. В наиболее популярных:
 +
 
 +
<math>{\color{Blue}\bullet} \quad </math> Mathematica - Orthogonalize[{v1,v2,…}] [https://reference.wolfram.com/language/tutorial/VectorOperations.html]
 +
 
 +
<math>{\color{Blue}\bullet} \quad </math> Maxima - gramschmidt(x) [http://maxima.sourceforge.net/docs/manual/maxima_23.html]
  
−
Существует реализация процесса отртогонализации Грама-Шмидта, которая является вычислительно более устойчивой, такая реализации имеет название Модифицированный процесс Грама-Шмидта.
+
<math>{\color{Blue}\bullet} \quad </math> MatLab - orth(A) [https://www.mathworks.com/help/symbolic/orth.html].
  
 
= Литература =
 
= Литература =

Текущая версия на 09:41, 15 декабря 2016

Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу
Данное задание было проверено и зачтено.
Проверено 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], в протоколах безопасности[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. Информационный граф алгоритма Грама-Шмидта
Рис. 2. Информационный граф алгоритма Грама-Шмидта

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. Это замечание позволяет уменьшить суммарное число операций, но увеличивает объём данных передаваемых между вычислительными элементами.

Рис. 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. График времени работы реализации алгоритма Грама-Шмидта (с использованием технологии MPI)
Рис. 5. График эффективности распараллеливания реализации алгоритма Грама-Шмидта относительно последовательной реализации
Таблица зависимости времени работы реализации от размерности задачи и числа вычислительных узлов

Время(сек.) | N | M | CPUS 

 8.214269    2000 2000 8
 5.535392    2000 2000 16
 7.121008    2000 2000 24
 10.995030   2000 2000 32
 17.709424   2000 2000 40

 15.969897   2500 2500 8
 9.438978    2500 2500 16
 9.658329    2500 2500 24
 12.944879   2500 2500 32
 17.741019   2500 2500 40
 25.414062   2500 2500 48

 27.463480   3000 3000 8
 15.192174   3000 3000 16
 13.546558   3000 3000 24
 14.893935   3000 3000 32
 19.726213   3000 3000 40
 27.127187   3000 3000 48
 35.446089   3000 3000 56

 43.528986   3500 3500 8
 23.267714   3500 3500 16
 18.832487   3500 3500 24
 20.074890   3500 3500 32
 22.985693   3500 3500 40
 27.376960   3500 3500 48
 37.711549   3500 3500 56
 48.896100   3500 3500 64

 64.889743   4000 4000 8
 34.014309   4000 4000 16
 25.912094   4000 4000 24
 24.668768   4000 4000 32
 27.262719   4000 4000 40
 32.295430   4000 4000 48
 41.851132   4000 4000 56
 48.027403   4000 4000 64

 92.075673   4500 4500 8
 47.766134   4500 4500 16
 35.229498   4500 4500 24
 31.972291   4500 4500 32
 33.277924   4500 4500 40
 38.653318   4500 4500 48
 45.888597   4500 4500 56
 54.089960   4500 4500 64

 126.197428  5000 5000 8
 64.872589   5000 5000 16
 46.852009   5000 5000 24
 40.500379   5000 5000 32
 39.817172   5000 5000 40
 42.630904   5000 5000 48
 48.521401   5000 5000 56
 55.565563   5000 5000 64

 85.984754   5500 5500 16
 60.770627   5500 5500 24
 50.246163   5500 5500 32
 47.708869   5500 5500 40
 49.754522   5500 5500 48
 52.557597   5500 5500 56
 59.170178   5500 5500 64

 111.130398  6000 6000 16
 77.988458   6000 6000 24
 63.733449   6000 6000 32
 57.554799   6000 6000 40
 57.613826   6000 6000 48
 62.380190   6000 6000 56
 68.328533   6000 6000 64

 121.137472  7000 7000 24
 96.815082   7000 7000 32
 84.331911   7000 7000 40
 79.391282   7000 7000 48
 78.346921   7000 7000 56
 84.036373   7000 7000 64


На Рис. 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
Рис. 6. График времени работы параллельной реализации алгоритма Грама-Шмидта при размерности исходной матрицы 5000x5000

2.4.3 Особенности тестирования реализации алгоритма

Для исследования возможностей параллелизма алгоритма ортогонализации Грама-Шмидта была написана собственная реализация [3] с использованием технологии MPI. В данной реализации задачи вычисления векторов [math]b_{j}[/math] равномерно распределяются по всем доступным вычислительным узлам. Как только на одном из узлов будет полностью вычислен некоторый [math]b_{j}[/math], происходит рассылка результата вычисления на все остальные узлы.

Тестирование данной реализации производилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса МГУ им. Ломоносова.

Технические подробности

Компилятор: openmpi/1.8.4-gcc

Планировщик задач: slurm/2.5.6

Очереди: test, regular4

Команда постановки скрипта в очередь: sbatch -n <CPUS> --time=20 ompi a.out

Тестовые данные (матрицы различной размерности) были сгенерированы с помощью скрипта Gram_Schmidt_process/scripts/matrixGenerator.py из репозитория [4]. Генерировались случайные квадратные матрицы выбранной размерности.

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 Литература

  1. ↑ А.А. Самарский, А.В. Гулин, Численные методы, Москва: «НАУКА», 1989
  2. ↑ А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления, Томск: Томский политехнический университет, 2006
  3. ↑ А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности, Санкт-Петербург: НИУ ИТМО, 2016
  4. ↑ Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов, Москва: МГУ ВМК, 2009
  5. ↑ А.Н. Канатников, А.П.Крищенко, Линейная Алгебра, Москва: МГТУ им. Н.Э. Баумана, 2002
  6. ↑ Jim Lambers, Gram-Schmidt Orthogonalization, The University of Southern Mississippi, Department of Mathematics, 2014
  7. ↑ Å. Björck, Numerics of Gram-Schmidt orthogonalization, Linkiiping University, Department of Mathematics, 1994
  8. ↑ Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization, Massachusetts Institute of Technology, 2006
  9. ↑ С.А. Белов, Н.Ю. Золотых, Численные Методы Линейной Алгебры, Нижний Новгород, 2008
  10. ↑ Воеводин Вл. В., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. В., Практика суперкомпьютера «Ломоносов», Москва: «Открытые системы», 2012
  11. ↑ 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
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Pandvik/Ортогонализация_Грама_-_Шмидта&oldid=20967