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

Эта работа прошла предварительную проверку Дата последней правки страницы: 02.11.2016 Данная работа соответствует формальным критериям. Проверено ASA.

Ортогонализация Грама-Шмидта
Последовательный алгоритм
Последовательная сложность	$O\left(k^2n\right)$
Объём входных данных	$kn$
Объём выходных данных	$kn$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(k\sqrt{n})$
Ширина ярусно-параллельной формы	$O(kn)$

Основные авторы описания: А.В.Кибанов, Т.З.Аджиева.

Все пункты данной статьи обсуждались и создавались совместно, доля ответственности равноправна.

1 Свойства и структура алгоритма

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

Процесс ортогонализации Грама–Шмидта — алгоритм построения для данной линейно независимой системы векторов евклидова или эрмитова пространства $V$ ортогональной системы ненулевых векторов, при котором каждый вектор построенной системы линейно выражается через векторы данной системы. Процесс ортогонализации Грама-Шмидта может быть применен к любой, в том числе и линейно-зависимой системе элементов евклидова пространства. Если ортогонализуемая система линейно-зависима, то на некотором шаге (возможно несколько раз) мы получим нулевой элемент, после отбрасывания которого можно продолжить процесс ортогонализации, однако в данной статье этот случай мы не рассматриваем.

Исторически это первый алгоритм, описывающий процесс ортогонализации. Традиционно его связывают с именами Йоргена Педерсена Грама и Эрхарда Шмидта. Йорган Педерсен Грам был датским актуарием, который в неявном виде представил суть процесса ортогонализации в 1883 году. Нельзя не сказать о том, что метод ортогонализации в несколько ином, модифицированом виде, ранее использовал Пьер-Симон Лаплас при нахождении массы Сатурна и Юпитера из систем нормальных уравнений и расчете распределения ошибок при этих вычислениях. Ясно, что Й. Грам не знал о том, что Лаплас ранее использовал этот метод. Эрхард Шмидт был учеником великих математиков Германа Шварца и Давида Гильберта. Алгоритм ортогонализации был опубликован Э. Шмидтом в 1907 году в его исследовании интегральных уравнений, которое в свою очередь дало основание для развития теории гильбертовых пространств. Шмидт использовал процесс ортогонализации применительно к геометрии гильбертова пространства решений интегральных уравнений отмечал, что процесс ортогонализации принципиально такой же, какой прежде использовал Й. Грам. В отечественной литературе иногда этот процесс называют ортогонализацией Сонина–Шмидта. Это название связано с тем, что разработанный Сониным метод ортогонализации системы функций для частного случая по существу не отличается от метода ортогонализации системы функций, предложенного ранее Шмидтом.

Классический процесс ортогонализации Грама–Шмидта применяется в таких совменных областях науки, как анализ изображений, цифровая обработка сигналов, нейронные сети.

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

Bходные данные: $k$ линейно независимых векторов $\mathbf{a_1},\mathbf{a_2},\ldots,\mathbf{a_k}$ длины $n$ $\left(\alpha_{ij}\right.$, $j=1,2, \ldots,n$, — координаты вектора $\left.\mathbf{a_i}\right)$ .

Вычисляемые данные:

$k$ ортогональных векторов $\mathbf{b_1},\mathbf{b_2},...,\mathbf{b_k}$ длины $n$, причем $\mathbf{b_1}=\mathbf{a_1}$, либо

$k$ ортонормированных векторов $\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}$ длины $n$, причем $\mathbf{q_i}=\frac{\mathbf{b_i}}{|\mathbf{b_i}|}$, $i=1,2, \ldots,k$.

Формулы процесса ортогонализации:

$\begin{align} \mathbf{b_{1}} & =\mathbf{a_{1}}, \\ \mathbf{b_{2}}& =\mathbf{a_{2}}-\mathbf{b_{1}} proj_{\mathbf{b_1}}\mathbf{a_{2}},\\ & ...\\ \mathbf{b_{i}} & =\mathbf{ a_{i}}-\sum\limits_{j=1}^{i-1}\mathbf{b_{i}} proj_{\mathbf{b_j}}\mathbf{a_{i}}.\\ & ...\\ \mathbf{b_{k}} & =\mathbf{ a_{k}}-\sum\limits_{j=1}^{k-1}\mathbf{b_{k}} proj_{\mathbf{b_j}}\mathbf{a_{k}}.\\ \end{align}$

Здесь $proj_{\mathbf{b_j}}\mathbf{a_{i}}$, для $j=1,...,k-1$ — это число, равное по величине проекции вектора $\mathbf{a_{j}}$ на ось, проходящую через вектор $\mathbf{b_j}$.

Формула для ее вычисления, полученная из определения скалярного произведения: $proj_{\mathbf{b_j}}\mathbf{a_{i}}=\dfrac{(\mathbf{ai},\mathbf{bj})}{\mathbf{\left \|\mathbf{b_j}\right \|}}$, для $j=1,...,i-1$

Произведение длин $\mathbf{|b_1|},\mathbf{|b_2|},\ldots ,\mathbf{|b_k|}$ равно объему параллелепипеда, построенного на векторах системы $\left\{\mathbf{a_i}\right\}$, $i=1,2,\ldots ,k$, как на ребрах.

Явное выражение векторов $\mathbf{b_i}$ для $i=1,...,k$ через $\mathbf{a_1},\mathbf{a_2},...,\mathbf{a_k}$ дает формула

$\mathbf{b_i}= \begin{vmatrix} (a_1,a_1) & \cdots & (a_1,a_i-1) &a_1 \\ (a_{2},a_1) & \cdots & (a_2,a_i-1) &a_2 \\ \cdots & \cdots & \cdots& \cdots \\ (a_{i},a_1) & \cdots & (a_i,a_i-1) &a_i \\ \end{vmatrix}$. (В правой части этого равенства определитель следует формально разложить по последнему столбцу).

Иногда полученные векторы нормируются сразу после их нахождения и находится система ортонормированных векторов $\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}$. В этом случае знаменатель в формуле для вычисления проекции $|\mathbf{q_i}|=1$ для $i=1,...,k-1$, что существенно упрощает вычисления алгоритма. В настоящей статье применяя процесс ортогонализации к системе линейно независимых векторов, мы так и будем поступать.

Если в описанном алгоритме в формулах для вычисления векторов $\mathbf{b_i}\,$, для $i=1,...,k$ процесса ортогонализации заменить $proj_{\mathbf{b_j}}\mathbf{a_{i}}$, для $j=1,...,i-1$ на $proj_{\mathbf{b_j}}\mathbf{b_{i}}$, для $j=1,...,i-1$,алгоритм, полученный таким образом называется модифицированным алгоритмом Грама–Шмидта.

Чаще всего, процесс Грама–Шмидта приводит к значительному нарушению ортогональности. Что касается модифицированного процесса Грама–Шмидта, он проявляет себя достаточно неусточиво, хоть и чуть менее, чем классический процесс.

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

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

• $\dfrac{k(k+1)}{2}$ скалярных произведений векторов,

включая нахождения скалярных произведений одинаковых векторов дальнейшего нахождения их длин;

• $\dfrac{(k-1)k}{2}$ умножений векторов на число после вычисления нового базисного вектора на итерации.

1.4 Макроструктура алгоритма

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

Таким образом, основновные макрооперации алгоритма:

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

Теперь опишем макроструктуру процесса ортогонализации. На верхнем уровне она представляет из себя последовательное вычисление $k$ ортонормированных векторов $\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_i},\ldots,\mathbf{q_k}$.

1. При $i=1$:

$\mathbf{q_{1}}= \frac{\mathbf{a_{1}}}{len(\mathbf{a_1})}$.

2. При $i=2,...,k$

$\mathbf{b_{i}}=\mathbf{a_{i}}-\sum\limits_{k=1}^{i-1}(\mathbf{ai},\mathbf{q_k})\mathbf{q_{k}}$

$\mathbf{q_{i}}=\frac{\mathbf{b_{i}}}{len(\mathbf{b_i})}$

Здесь $len(\cdot)$ — обозначение длины вектора, вычисляемой в $\mathbb{R}^n$ стандартным образом, как арифметический квадратный корень из скалярного произведения.

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

Входные данные: $k$ линейно независимых векторов $\mathbf{a_1},\mathbf{a_2},\ldots,\mathbf{a_k}$ длины $n$ $\left(\alpha_{ij}\right.$, где $j=1,2, \ldots,n$ — координаты вектора $\mathbf{a_i}$, для $\left.i=1,2, \ldots,k\right)$.

Вычисляемые данные:

$k$ ортогональных векторов $\mathbf{b_1},\mathbf{b_2},...,\mathbf{b_k}$ длины $n$ $\left(\beta_{ij}\right.$, где $j=1,2, \ldots,n$ — координаты вектора $\mathbf{b_i}$, для $\left.i=1,2, \ldots,k\right)$, причем $\mathbf{b_1}=\mathbf{a_1}$, либо

$k$ ортонормированных векторов $\mathbf{q_1},\mathbf{q_2},\ldots,\mathbf{q_k}$ длины $n$ ($\gamma_{ij}$, $j=1,2, \ldots,n$ — координаты вектора $\mathbf{q_i}$, для $\left.i=1,2, \ldots,k\right)$, причем $\mathbf{q_1}=\frac{\mathbf{b_i}}{|\mathbf{b_i}|}$, для $i=1,2, \ldots,k$.

Последовательность исполнения метода:

1. $len(\mathbf{a_1})=\sqrt{\alpha_{11}^2+\alpha_{12}^2+\ldots+\alpha_{1n}^2}$.

2. При $i$ от $1$ до $n$:

$\gamma_{1j}= \frac{\alpha_{1j}}{len(\mathbf{a_1})}$.

3. При $i$ от $2$ до $k$:

при $j$ от $1$ до $n$:

$\beta_{ij}=\alpha_{ij}-\sum\limits_{k=1}^{j-1}\left(\sum\limits_{m=1}^n \alpha_{im}\gamma_{km}\right)\gamma_{kj}$

$len(\mathbf{b_i})=\sqrt{\beta_{i1}^2+\beta_{i2}^2+\ldots+\beta_{in}^2}$

при $j$ от $1$ до $n$:

$\gamma_{ij}= \frac{\beta_{ij}}{len(\mathbf{b_{i}})}$.

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

Для ортогонализации $n$ линейно независимых векторов длины $n$ в последовательном варианте с помощью классического метода ортогонализации Грама-Шмидта потребуется:

• $k$ вычислений квадратного корня
• $\frac{k(k-1)(n-1)}{2}$ сложений,
• $\frac{k(k-1)(n-1)}{2}$ вычитаний,
• $k^2 n$ умножений,
• $k$ делений.

Умножения, сложения и вычитания составляют основную часть алгоритма.

При классификации по последовательной сложности, таким образом, метода ортогонализации Грама-Шмидта с кубической сложностью.

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

Представим граф алгоритма с помощью текстового описания, а также с помощью изображения.

Все вершины графа можно разбить на пять групп, каждой из которых соответствует вид исполняемой операции. Каждая вершина расположена в точках с целочисленными координатами, в пространствах с числом измерений от $1$ до $3$. Дадим описание каждой из групп.

Первая группа — вершины, которым соответствует операция скалярного умножения вектора на себя. На графе они располагаются в одномерной области. $i$ изменяется от $1$ до $k$.

Значениями аргументов для этих функций являются:

• при $i = 1$ — входные данные $\{a_{1p}\}, p = 1..n$;
• при $i \gt 1$ — результат выполнения функции в вершине пятой группы с координатами $\{i-1,1,p\}, p = 1..n$.

Результат этих операций является промежуточным данным алгоритма.

Вторая группа — вершины, которым соответствует операция SQRT извлечения квадратного корня. Они расположены в одномерной области, $i$ изменяется от $1$ до $k$.

Значением аргумента для этих функций является результат выполнения функции в вершине первой группы, с координатой $i$.

Результат является промежуточным данным алгоритма.

Третья группа – вершины, которым соответствует операция деления. Располагаются в двумерной области, $i$ изменяется от $1$ до $k$, $p$ изменяется от $1$ до $n$.

Значениями аргументов для этих функций являются:

• делимое:
  • при $i = 1$ — входные данные $a_{1p}$;
  • при $i \gt 1$ — результат выполнения функции в вершине пятой группы с координатами $i-1,1,p$;
• делитель - результат выполнения функции в вершине второй группы с координатой $i$.

Результат является выходным данным алгоритма $q_{ip}$.

Четвертая группа — вершины, которым соответствует операция скалярного умножения векторов $\mathbf{a}\cdot\mathbf{b}$. Располагаются в двумерной области, $i$ изменяется от $1$ до $k-1$, $j$ изменяется от $1$ до $k-i$.

Значениями аргументов для этих функций являются:

• $\mathbf{a}$ — входные данные $\{a_{1+j,p}\},p = 1..n$;
• $\mathbf{b}$ — результат выполнения функции в вершине третьей группы с координатой $\{i,p\}, p = 1..n$.

Результат является промежуточным данным алгоритма.

Пятая группа – вершины которые соответствуют операции $a - b c$. Располагаются в трехмерной области, $i$ изменяется от $1$ до $k-1$, $j$ изменяется от $1$ до $k-i$, $p$ изменяется от $1$ до $n$.

Значениями аргументов для этих функций являются:

• $a$:
  • при $i = 1$ — входные данные $a_{1+j,p}$;
  • при $i \gt 1$ — результат выполнения функции в вершине пятой группы с координатами $i-1,j+1,p$;
• $b$ — результат выполнения функции в вершине четвертой группы, с координатами $i,j$;
• $c$ — результат выполнения функции в вершине третьей группы с координатой $i,p$.

Результат является промежуточным данным алгоритма.

Описанный граф для случая $k = 3$, $n = 3$ изображен на рисунке. Каждая группа вершин отличается цветом и буквенным обозначением. "S" — первая группа, "SQ" — вторая группа, "Div" — третья группа, "Dot" — четвертая группа, "f" — пятая группа. Квадратные метки синего и розового цветов обозначают передачу входных и выходных данных соответственно.

Рисунок 1. Граф алгоритма с отображением входных и выходных данных

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

Для построения ортонормированного базиса методом Грама–Шмидта необходимо выполнить следующие ярусы операций:

• $k$ ярусов вычисления скалярных произведений векторов на самих себя (по одной операций на каждом, сложность операции по высоте и по ширине ЯПФ — корень квадратный);
• $k$ ярусов вычисления квадратного корня (по одной операции на каждом);
• $k$ ярусов деления (по $n$ операций на каждом);
• $k-1$ ярус вычисления скалярных произведений пары векторов (от $1$ до $(k-1)$ операций на каждом,сложность операции по высоте и по ширине ЯПФ — корень квадратный);
• $k-1$ ярус умножения/вычитания чисел (от $n$ до $(k-1)n$ операций на каждом).

В этом случае сложность алгоритма при классификации по высоте ЯПФ — $O(k\sqrt{n})$, при классификации по ширине ЯПФ — $O(kn)$.

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

Входные данные: матрица $A$ с элементами $\alpha_{ij}$, где $i$ принимает значения $1,\ldots,k$, $j$ принимает значения $1,\ldots,n$. Здесь $k\leqslant n$. По строкам этой матрицы записаны $k$ линейно-независимых векторов некоторого пространства размерности $n$.

Объем входных данных: $kn$.

Выходные данные: матрица $Q$ с элементами $\gamma_{ij}$. Строки этой матрицы также задают $k$ векторов единичной длины размерности $n$, причем всякое скалярное произведение двух различных векторов этой системы равно $0$.

Объем выходных данных: $kn$.

Замечание: Алгоритм может работать с вырожденной матрицей, он выдаёт нулевую строку на шаге $j$, если строка $j$ матрицы $A$ является линейной комбинацией предыдущих строк. В этом случае для корректного продолжения работы алгоритма необходима дополнительная проверка на равенство строки нулю и пропуск соответствующих строк. Количество построенных векторов будет равняться размерности пространства, порождаемого строками матрицы.

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

Отношение последовательной и параллельной сложности алгоритма в случае неограниченных ресурсов равно $\dfrac{k^2n}{k\sqrt{n}}=k\sqrt{n}$. Вычислительная мощность алгоритма линейная.

Алгоритм почти полностью детерминирован — единственность результата выполнения гарантирована, однако возможно накопление ошибок округления.

Алгоритм является численно неустойчивым — ошибки округления могут привести к неортогональности полученных векторов.

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

2 Программная реализация алгоритма

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

Для классического метода ортогонализации Грама–Шмидта существуют реализации в математических программных пакетах, таких как

• Maxima,
• Mathematica,
• Matlab.

Также метод реализован в библиотеке SpectralPython для языка Python. Для устранения проблемы ошибок численных округлений иногда используется повторное применение метода к уже построенному ортогональному базису.

Поскольку метод ортогонализации Грама–Шмидта гарантирует получение $j$ ортогональных векторов к концу шага $j$, а не в самом конце работы алгоритма, используется в итерационных методах, таких как метод Арнольди. Для поддержки работы этого метода реализация ортогонализации Грама–Шмидта включена в библиотеку ARPACK.