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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1. [math]\displaystyle{ \mathbf{b}_1 = \mathbf{a}_1 }[/math]

2. [math]\displaystyle{ \mathbf{b}_i = \mathbf{a}_i - \sum_{j=1}^{i-1} \mathbf{proj}_{\mathbf{b}_{i-1}} \mathbf{a}_i }[/math] [math]\displaystyle{ (i = 2 \cdots N) }[/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 \frac{n(n-1)}{2}n }[/math] - вычитаний;

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

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

Исходя из вышесказанного, получаем итоговую сложность последовательного алгоритма [math]\displaystyle{ O(n^3) }[/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].

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

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

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

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

На основе графа зависимости по данным, представленном на Рис. 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{ N*3 }[/math] умножения, [math]\displaystyle{ (N-1)*2 }[/math] сложения и 1 деление. В цепочке зависимости данных присутствует [math]\displaystyle{ N-1 }[/math] проекция. Значит всего [math]\displaystyle{ N*N*3 }[/math] умножения, [math]\displaystyle{ N*(N-1)*2 }[/math] сложения и [math]\displaystyle{ N*1 }[/math] деления. Получаем итоговую сложность [math]\displaystyle{ O(N^{2}) }[/math].

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

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

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