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

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

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

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

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

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

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

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

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

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

Исходные данные: квадратная матрица A с линейно независимыми векторами [math]\mathbf{a}_1,...,\mathbf{a}_N[/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].

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

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

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

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

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].

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

Последовательный алгоритм требует:

[math]{\color{Blue}\bullet} \quad \frac{n(n-1)}{2}n[/math] - вычитаний;

[math]{\color{Blue}\bullet} \quad \frac{n(n-1)}{2}n[/math] - делений;

[math]{\color{Blue}\bullet} \quad n \cdot 2 \cdot (n-1)(n-2)[/math] - умножений.

Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма [math]O(n^3)[/math].

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

Опишем информационный граф алгоритма.

Для вычисления [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].

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

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

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

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 N \cdot 3[/math] умножения в одной проекции;

[math]{\color{Blue}\bullet} \quad (N-1) \cdot 2[/math] сложения и 1 деление.

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

[math]{\color{Blue}\bullet} \quad N \cdot N \cdot 3[/math] - умножений;

[math]{\color{Blue}\bullet} \quad N \cdot (N-1) \cdot 2[/math] - сложений;

[math]{\color{Blue}\bullet} \quad N \cdot 1[/math] - делений.

Получаем итоговую сложность [math]O(N^{2})[/math].

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

Входные данные: квадратная матрица A (элементы [math]\mathbf{a}_{ij}[/math]). Матрица может быть приведена к верхнему(нижнему) треугольному виду.

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

Объём входных данных: [math]\frac{n(n+1)}{2}[/math] (после преобразования матрицы к треугольному виду достаточно хранить только ненулевые элементы).

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

Объём выходных данных: [math]\frac{n(n+1)}{2}[/math] (в силу треугольного вида матрицы).

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов является линейной(последовательная реализация - [math]O(n^3)[/math], параллельная реализация - [math]O(n^2)[/math]).

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

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

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

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

Пример реализации базового алгоритма с использованием технологий OpenMP и MPI представлен в статье [7].

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

3 Литература

</references>

1. А.А. Самарский, А.В. Гулин, Численные методы.

2. А.Е. Карелин, А.А. Светлаков, Использование ортогонализации Грама-Шмидта для повышения экономичности многоточечных алгоритмов рекуррентного оценивания параметров моделей объектов управления.

3. The University of Southern Mississippi, Department of Mathematics, Jim Lambers, Gram-Schmidt Orthogonalization.

4. Massachusetts Institute of Technology, Per-Olof Persson, Introduction to Numerical Methods: Gram-Schmidt Orthogonalization.

5. МГУ ВМК, Ф.Ф. Дедус, Л.И. Куликова, А.Н. Панкратов, Р.К. Тетуев, Классические ортогональные базисы в задачах аналитического описания и обработки информационных сигналов.

6. А.В. Пискова, А. А. Менщиков, А.Г. Коробейников, Использование ортогонализации Грама-Шмидта в алгоритме приведения базиса решётки для протоколов безопасности.

7. Tirana University, Department of Mathematics, Genci Berati, Gram–Schmidt Process in Different Parallel Platforms (Control Flow versus Data Flow).

8. The University of Iowa, Department of Computer Science, S. Oliveira, L. Borges, M. Holzrichter, T. Soma, Analysis of Different Partitioning Schemes for Parallel Gram-Schmidt Algorithms.

9. Islamic Azad University (IAU), Mechanical and Aerospace Engineering Department, Science and Research Branch, S. Roholah Ghodsi, Bahman Mehri, Mohammad Taeibi-Rahni, A Parallel implementation of Gram-Schmidt Algorithm for Dense Linear System of Equations.