Участник:DVIN/Ортогонализация Грамма-Шмидта: различия между версиями
DVIN (обсуждение | вклад) (Новая страница: «''' Ортогонализация Грамма-Шмидта ''' == Свойства и структура алгоритмов == === Общее описани…») |
ASA (обсуждение | вклад) |
||
(не показано 37 промежуточных версий 3 участников) | |||
Строка 1: | Строка 1: | ||
− | + | {{Assignment|ASA}} | |
+ | |||
+ | {{algorithm | ||
+ | | name = Ортогонализация Грама-Шмидта | ||
+ | | serial_complexity = <math>O(N^2k)</math> | ||
+ | | pf_height = <math>O(N)</math> | ||
+ | | pf_width = <math>O(Nk)</math> | ||
+ | | input_data = <math>Nk</math> | ||
+ | | output_data = <math>Nk</math> | ||
+ | }} | ||
+ | |||
+ | Основные авторы описания: [[Участник:DVIN|Инжелевская Дарья Валерьевна]] | ||
== Свойства и структура алгоритмов == | == Свойства и структура алгоритмов == | ||
=== Общее описание алгоритма === | === Общее описание алгоритма === | ||
+ | Создатели алгоритма Эрхард Шмидт и Йорген Педерсен Грам. | ||
+ | |||
+ | '''Эрхард Шмидт''' (14 января 1876 — 6 декабря 1959) — немецкий математик, с 1917 года профессор Берлинского университета. В 1946—58 первый директор Института математики АН ГДР. Основные труды по теории функций, интегральным уравнениям, функциональному анализу. Определил и изучил геометрически гильбертово пространство, используя аналогии с геометрией Евклида. Занимался квадратичными формами. Известен оператор Гильберта-Шмидта и процесс Грама ― Шмидта. | ||
+ | |||
+ | '''Йорген Педерсен Грам''' (27 июня 1850 - 29 апреля 1916) датский математик, родился в герцогстве Шлезвиг, Дания и умер в Копенгагене, Дания. Основные его работы "On series expansions determined by the methods of least squares", "Investigations of the number of primes less than a given number". Математический метод, который носит его имя, процесс Грама-Шмидта, впервые был опубликован в 1883. Теорема Грама и Определитель Грама также назван в его честь. | ||
+ | |||
+ | Процесс Грама-Шмидта — это один из алгоритмов, в которых на основе множества линейно независимых векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math> строится множество ортогональных векторов <math>{\displaystyle \mathbf {b}_{1},\;\ldots ,\;\mathbf {b} _{N}} </math> или ортонормированных векторов <math>{\displaystyle \mathbf {e} _{1},\;\ldots ,\;\mathbf {e}_{N}} </math>, причём так, что каждый вектор <math>{\displaystyle \mathbf {b} _{j}} </math> или <math>{\displaystyle \mathbf {e} _{j}}</math> может быть выражен линейной комбинацией векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\; \mathbf {a} _{j}}</math>. | ||
+ | |||
+ | === Математическое описание алгоритма === | ||
+ | ==== Классический процесс Грама — Шмидта ==== | ||
+ | |||
+ | ===== Алгоритм ===== | ||
+ | Пусть имеются линейно независимые векторы <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N</math>. | ||
+ | |||
+ | Определим оператор проекции следующим образом: <math>\mathbf{proj}_{\mathbf{b}}\,\mathbf{a} = {\langle \mathbf{a}, \mathbf{b} \rangle \over \langle \mathbf{b}, \mathbf{b}\rangle} \mathbf{b} ,</math> | ||
+ | |||
+ | где <math>\langle \mathbf{a}, \mathbf{b} \rangle</math> — скалярное произведение векторов <math>\mathbf{a}</math> и <math>\mathbf{b}</math>. | ||
+ | |||
+ | Скалярное произведение для двух векторов <math>\mathbf{ a= [a_1, a_2, ...,a_k]}</math> и <math>\mathbf{ b= [b_1, b_2, ..., b_k]}</math> в '''''k'''''-мерном действительном пространстве определяется как: | ||
+ | :<math>\langle \mathbf{a}, \mathbf{b} \rangle=\sum_{i=1}^k a_ib_i=a_1b_1+a_2b_2+\cdots+ a_kb_k</math>. | ||
+ | |||
+ | Этот оператор проецирует вектор <math>\mathbf{a}</math> коллинеарно вектору <math>\mathbf{b}</math>. | ||
+ | |||
+ | Ортогональность векторов <math>\mathbf{a}</math> и <math>\mathbf{b}</math> достигается на шаге (2). | ||
+ | |||
+ | Классический процесс Грама — Шмидта выполняется следующим образом: | ||
+ | |||
+ | : <math> | ||
+ | {\begin{array}{lclr} | ||
+ | {\mathbf {b}}_{1}&=&{\mathbf {a}}_{1}&(1)\\ | ||
+ | {\mathbf {b}}_{2}&=&{\mathbf {a}}_{2}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{2}&(2)\\ | ||
+ | {\mathbf {b}}_{3}&=&{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{2}}}\,{\mathbf {a}}_{3}&(3)\\ | ||
+ | {\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}&(4)\\ | ||
+ | &\vdots &&\\{\mathbf {b}}_{N}&=&{\mathbf {a}}_{N}-\displaystyle \sum _{{j=1}}^{{N-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{N}&(N) | ||
+ | \end{array}} | ||
+ | </math> | ||
+ | |||
+ | |||
+ | На основе каждого вектора <math>\mathbf{b}_j \;(j = 1 \ldots N)</math> может быть получен нормированный вектор: <math>\mathbf{e}_j = {\mathbf{b}_j\over \| \mathbf{b}_j \|}</math> (у нормированного вектора направление будет таким же, как у исходного, а длина — единичной). Причем берется норма согласованная со скалярным произведением <math>\| x \| = \sqrt{\langle x, x \rangle}</math> | ||
+ | |||
+ | Результаты процесса Грама — Шмидта: | ||
+ | |||
+ | <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math> — система ортогональных векторов либо | ||
+ | |||
+ | <math>\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N</math> — система ортонормированных векторов. | ||
+ | |||
+ | Вычисление <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math> носит название ортогонализации Грама — Шмидта, а <math>\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N</math> — ортонормализации Грама — Шмидта. | ||
+ | |||
+ | ===== Доказательство ===== | ||
+ | Докажем ортогональность векторов <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N</math>. | ||
+ | |||
+ | Для этого вычислим скалярное произведение <math>\langle \mathbf{b}_1, \mathbf{b}_2 \rangle</math>, подставив в него формулу (2). Мы получим ноль. Равенство нулю скалярного произведения векторов означает, что эти векторы ортогональны. Затем вычислим скалярное произведение <math>\langle \mathbf{b}_1, \mathbf{b}_3 \rangle</math>, используя результат для <math>\langle \mathbf{b}_1, \mathbf{b}_2 \rangle</math> и формулу (3). Мы снова получим ноль, то есть векторы <math>\mathbf{b}_1</math> и <math>\mathbf{b}_3</math> ортогональны. Общее доказательство выполняется методом математической индукции. | ||
+ | |||
+ | ===== Геометрическая интерпретация — вариант 1 ===== | ||
+ | [[Файл:Gram-Schmidt-step2.svg|right|frame|Рис. 1. Второй шаг процесса Грама — Шмидта]] | ||
+ | Рассмотрим формулу (2) — второй шаг алгоритма. | ||
+ | Её геометрическое представление изображено на рис. 1: | ||
+ | |||
+ | 1 — получение проекции вектора <math>\mathbf{a}_2</math> на <math>\mathbf{b}_1</math>; | ||
+ | |||
+ | 2 — вычисление <math>\mathbf{a}_2-\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_2</math>, то есть [[перпендикуляр]]а, которым выполняется проецирование конца <math>\mathbf{a}_2</math> на <math>\mathbf{b}_1</math>. | ||
+ | |||
+ | Этот перпендикуляр — вычисляемый в формуле (2) вектор <math>\mathbf{b}_2</math>; | ||
+ | |||
+ | 3 — перемещение полученного на шаге 2 вектора <math>\mathbf{b}_2</math> в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (2). | ||
+ | |||
+ | На рисунке видно, что вектор <math>\mathbf{b}_2</math> ортогонален вектору <math>\mathbf{b}_1</math>, так как <math>\mathbf{b}_2</math> является перпендикуляром, по которому <math>\mathbf{a}_2</math> проецируется на <math>\mathbf{b}_1</math>. | ||
+ | |||
+ | Рассмотрим формулу (3) — третий шаг алгоритма — в следующем варианте: | ||
+ | : <math> | ||
+ | \begin{array}{lcr} | ||
+ | \mathbf{b}_3 = \mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3+\mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3) & & (6) \\ | ||
+ | \end{array} | ||
+ | </math> | ||
+ | |||
+ | Её геометрическое представление изображено на рис. 2: | ||
+ | [[Файл:Gram-schmidt-step3base1.png|center|frame|Рис. 2. Третий шаг процесса Грама — Шмидта]] | ||
+ | |||
+ | 1 — получение проекции вектора <math>\mathbf{a}_3</math> на <math>\mathbf{b}_1</math>; | ||
+ | |||
+ | 2 — получение проекции вектора <math>\mathbf{a}_3</math> на <math>\mathbf{b}_2</math>; | ||
+ | |||
+ | 3 — вычисление суммы <math>\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3</math>, то есть проекции вектора <math>\mathbf{a}_3</math> на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>. Эта плоскость закрашена на рисунке серым цветом; | ||
+ | |||
+ | 4 — вычисление <math>\mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3)</math>, то есть перпендикуляра, которым выполняется проецирование конца <math>\mathbf{a}_3</math> на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>. Этот перпендикуляр — вычисляемый в формуле (6) вектор <math>\mathbf{b}_3</math>; | ||
+ | |||
+ | 5 — перемещение полученного <math>\mathbf{b}_3</math> в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (6). | ||
+ | |||
+ | На рисунке видно, что вектор <math>\mathbf{b}_3</math> ортогонален векторам <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>, так как <math>\mathbf{b}_3</math> является перпендикуляром, по которому <math>\mathbf{a}_3</math> проецируется на плоскость, образуемую векторами <math>\mathbf{b}_1</math> и <math>\mathbf{b}_2</math>. | ||
+ | |||
+ | Таким образом, в процессе Грама — Шмидта для вычисления <math>\mathbf{b}_j</math> выполняется проецирование <math>\mathbf{a}_j</math> ортогонально на [[гиперплоскость]]?, формируемую векторами <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}</math>. Вектор <math>\mathbf{b}_j</math> затем вычисляется как разность между <math>\mathbf{a}_j</math> и его проекцией. То есть <math>\mathbf{b}_j</math> — это перпендикуляр? от конца <math>\mathbf{a}_j</math> к гиперплоскости?, формируемой векторами <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}</math>. Поэтому <math>\mathbf{b}_j</math> ортогонален векторам, образующим эту гиперплоскость?. | ||
+ | |||
+ | ===== Численная неустойчивость ===== | ||
+ | При вычислении на ЭВМ по формулам (1) — (5) векторы <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{N-1}</math> часто не точно ортогональны из-за [[Ошибка квантования|ошибок округления]]. Из-за потери ортогональности в процессе вычислений классический процесс Грама — Шмидта называют [[Численная устойчивость|численно неустойчивым]]. | ||
+ | |||
+ | ==== Модифицированный процесс Грама — Шмидта ==== | ||
+ | Процесс Грама — Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Вместо вычисления <math>\mathbf{b}_j</math> как: <math> \mathbf{b}_j = \mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_j - \ldots - \mathbf{proj}_{\mathbf{b}_{j-1}}\,\mathbf{a}_j \; (7) </math> этот вектор вычисляется следующим образом: | ||
+ | : <math> \begin{array}{lclr} | ||
+ | \mathbf{a}_j^{(1)} & = &\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j & (8) \\ | ||
+ | \mathbf{a}_j^{(2)} & = &\mathbf{a}_j^{(1)} - \mathbf{proj}_{\mathbf{b}_2} \, \mathbf{a}_j^{(1)} & (9) \\& \vdots & & \\ | ||
+ | \mathbf{a}_j^{(j-2)} & = & \mathbf{a}_j^{(j-3)} - \mathbf{proj}_{\mathbf{b}_{j-2}} \, \mathbf{a}_j^{(j-3)} & (10) \\ | ||
+ | \mathbf{b}_j & = & \mathbf{a}_j^{(j-2)} - \mathbf{proj}_{\mathbf{b}_{j-1}} \, \mathbf{a}_j^{(j-2)} & (11) | ||
+ | \end{array} </math> | ||
+ | |||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
+ | |||
+ | Вычислительное ядро последовательной версии метода ортогонализации Грамма-Шмидта можно составить из множественных (всего их <math>\frac{N(N-1)}{2}</math>) вычислений проекций : <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> | ||
+ | |||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
+ | Данный алгоритм использует в качестве составных частей другие алгоритмы. В дальнейшем имеет смысл описывать алгоритм не в максимально детализированном виде (т.е. на уровне арифметических операций), а давать только его макроструктуру, то есть описывает структуру и состав макроопераций. Типичной макрооперацией, часто встречающиеся в алгоритме является оператор проекции векторов. Как записано и в описании ядра алгоритма, основную часть метода ортогонализации Грамма-Шмидта составляют множественные (всего их <math>\frac{N(N-1)}{2}</math>) вычисления оператора проекции. | ||
+ | |||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
+ | |||
+ | Следующий алгоритм реализует нормализацию Грамма-Шмидта. Векторы <math>v_1,...,v_k</math> заменяются набором ортонормированных векторов, которые имеют ту же линейную оболочку. | ||
+ | [[File:Gram-Schmidt ortho.png|none|thumb]] | ||
+ | |||
+ | Вычислительная сложность этого <math>2Nk^2</math> операции с плавающей точкой, где N - размерность векторов. | ||
+ | |||
+ | Последовательность исполнения метода следующая: | ||
+ | |||
+ | 1. <math>\mathbf {b}_{1}=\mathbf {a}_{1}</math> | ||
+ | |||
+ | 2. Далее для всех векторов <math>\mathbf {b}_{i}</math> для <math>i=2 ... N</math> производим вычисление по следующей формуле: | ||
+ | <math>{\mathbf {b}}_{i}={\mathbf {a}}_{i}-\displaystyle \sum _{{j=1}}^{{i-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}</math>. | ||
+ | |||
+ | В ней на каждом шаге <math>i</math> по очереди вычисляются все <math>{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}</math> для <math>j=1 ... i-1</math> | ||
+ | |||
+ | Пример реализации на Python. Функция работает для произвольного количества векторов любой размерности. При этом если количество векторов <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N</math> больше их размерности или они линейно зависимы то функция возвращает максимально возможное число линейно независимых векторов <math>\mathbf{b}_1,\;\ldots,\;\mathbf{b}_n</math>, а остальные векторы <math>\mathbf{b}_{n+1},\;\ldots,\;\mathbf{b}_N</math> нулевые. | ||
+ | <source lang="python"> | ||
+ | import random | ||
+ | |||
+ | def GramSchmidt(*a): | ||
+ | k=len(a[0]) | ||
+ | N=len(a); | ||
+ | b = [[0] * k for i in range(N)] | ||
+ | b[0]=a[0] | ||
+ | for i in range(1,N): | ||
+ | sum=a[i] | ||
+ | for j in range(0,i): | ||
+ | scolar_ab=0 | ||
+ | scolar_bb=0 | ||
+ | proj=[i for i in range(k)] | ||
+ | for n in range(k): | ||
+ | scolar_ab+=b[j][n]*a[i][n] | ||
+ | scolar_bb+=b[j][n]*b[j][n] | ||
+ | for n in range(k): | ||
+ | proj[n]=(scolar_ab/scolar_bb)*b[j][n] | ||
+ | for n in range(k): | ||
+ | sum[n]-=proj[n] | ||
+ | b[i]=sum | ||
+ | return b; | ||
+ | |||
+ | l1=[random.randrange(0,10) for i in range(3)] | ||
+ | l2=[random.randrange(0,10) for i in range(3)] | ||
+ | l3=[random.randrange(0,10) for i in range(3)] | ||
+ | print(l1,l2,l3) | ||
+ | print(GramSchmidt(l1,l2,l3)) | ||
+ | </source> | ||
+ | |||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
+ | |||
+ | Для построения ортогонального набора векторов в последовательном варианте требуется: | ||
+ | |||
+ | При k=N: | ||
+ | |||
+ | * <math>\frac{N(N-1)}{2}</math> делений, | ||
+ | |||
+ | * <math>\frac{N(N-1)(3N-1)}{2}</math> сложений (вычитаний), | ||
+ | |||
+ | * <math>\frac{3N^2 (N-1)}{2}</math> умножений. | ||
+ | |||
+ | При k≠N | ||
+ | |||
+ | * <math>\frac{N(N-1)}{2}</math> делений, | ||
+ | |||
+ | * <math>\frac{N(N-1)(3k-1)}{2}</math> сложений (вычитаний), | ||
+ | |||
+ | * <math>\frac{3Nk (N-1)}{2}</math> умножений. | ||
+ | |||
+ | При классификации по последовательной сложности, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам с кубической сложностью. | ||
+ | |||
=== Информационный граф === | === Информационный граф === | ||
+ | Опишем граф алгоритма как аналитически, так и в виде рисунка. | ||
+ | |||
+ | Граф алгоритма состоит из двух групп вершин, расположенных в целочисленных узлах двух областей. | ||
+ | |||
+ | '''Первая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i} = {\langle \mathbf{a_i}, \mathbf{b_j} \rangle \over \langle \mathbf{b_j}, \mathbf{b_j}\rangle} \mathbf{b_j} ,</math> . | ||
+ | |||
+ | Естественно введённые координаты области таковы: | ||
+ | |||
+ | i — меняется в диапазоне от 2 до N, принимая все целочисленные значения; | ||
+ | |||
+ | j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения. | ||
+ | |||
+ | Аргументы операции следующие: | ||
+ | |||
+ | <math>a_i</math>: | ||
+ | элементы входных данных, а именно <math>a_i</math>; | ||
+ | |||
+ | <math>b_j</math>: | ||
+ | результат срабатывания операции, соответствующей вершине из второй группы, с координатой j. | ||
+ | |||
+ | '''Вторая''' группа вершин расположена в двухмерной области, соответствующая ей операция <math>a - b</math>. | ||
+ | |||
+ | Естественно введённые координаты области таковы: | ||
+ | |||
+ | i — меняется в диапазоне от 2 до N, принимая все целочисленные значения; | ||
+ | |||
+ | j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения. | ||
+ | |||
+ | Аргументы операции следующие: | ||
+ | |||
+ | a: | ||
+ | |||
+ | j=1 - входные данные <math>a_j</math> | ||
+ | |||
+ | j>1 - результат срабатывания операции, соответствующей вершине из второй группы, с координатой j-1 | ||
+ | |||
+ | b: | ||
+ | |||
+ | результат срабатывания операции, соответствующей вершине из первой группы, с координатой j | ||
+ | [[Файл:picture00001.jpg|мини|upright=2|centre|frame|Рис. 5. Граф алгоритма с отображением входных и выходных данных. Proj - вычисление оператора проекции, F - операция a-b, In - входные данные, Out - результаты.]] | ||
+ | |||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
+ | |||
+ | В параллельном варианте, в отличие от последовательного, вычисления <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> можно выполнять параллельно. То есть после каждого вычисления <math>b_j</math> можно запускать параллельное вычисление <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> для всех <math>i=j+1..N</math>. Кроме того можно параллельно производить вычитание соответствующего <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> для всех <math>i=j+1..N</math>. | ||
+ | |||
+ | Для ортогонализации в параллельном варианте требуется последовательно выполнить следующие ярусы: | ||
+ | * <math>N</math> ярусов в каждом из которых будет выполнятся, | ||
+ | * по <math>N - 1</math> операторов проекций (в каждом из ярусов линейное количество проекций, в зависимости от яруса — от <math>N-1</math> до <math>1</math>), | ||
+ | * по <math>N - 1</math> сложений (в каждом из ярусов линейное количество сложений, в зависимости от яруса — от <math>N-1</math> до <math>1</math>), | ||
+ | |||
+ | При классификации по высоте (количество ярусов в ЯПФ ) ЯПФ, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам со сложностью <math>O(N)</math>. | ||
+ | |||
+ | При классификации по ширине (максимальное количество вершин в ярусе) ЯПФ его сложность будет <math>O(N^2)</math>. | ||
+ | |||
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
+ | |||
+ | '''Входные данные''': множества линейно независимых векторов <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math>, каждый из которых описывается <math>\mathbf{ a_i= [a^i_1, a^i_2, ..., a^i_k]}</math> . | ||
+ | |||
+ | Дополнительные ограничения: | ||
+ | * вектора <math>{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}</math> линейно независимы, поэтому <math>k \geqslant N</math>. | ||
+ | |||
+ | '''Объём входных данных''': <math>Nk</math> | ||
+ | |||
+ | '''Выходные данные''': множество ортогональных векторов <math>{\displaystyle \mathbf {b} _{1},\;\ldots ,\;\mathbf {b} _{N}} </math>, каждый из которых описывается <math>\mathbf{ b_i= [b^i_1, b^i_2, ..., b^i_k]}</math> . | ||
+ | |||
+ | '''Объём выходных данных''': <math>Nk</math> | ||
+ | |||
=== Свойства алгоритма === | === Свойства алгоритма === | ||
+ | |||
+ | Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение кубической к квадратичной). | ||
+ | |||
+ | При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, также линейна. | ||
+ | |||
+ | Алгоритм почти полностью детерминирован — единственность результата выполнения гарантирована, однако возможно накопление ошибок округления при использовании классического процесса Грама-Шмидта. При этом модифицированный алгоритм дает более хороший результат. | ||
+ | |||
+ | Алгоритм является численно неустойчивым — ошибки округления могут привести к неортогональности полученных векторов. | ||
+ | |||
+ | Процесс Грама — Шмидта может применяться также к бесконечной последовательности линейно независимых векторов. | ||
+ | |||
+ | Кроме того, процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт <math>\mathbf{0}</math> (нулевой вектор) на шаге <math>j</math>, если <math>\mathbf{a}_j</math> является линейной комбинацией векторов <math>\mathbf{a}_1,\;\ldots,\;\mathbf{a}_{j-1}</math>. Если это может случиться, то для сохранения ортогональности выходных векторов и для предотвращения деления на ноль при ортонормировании алгоритм должен делать проверку на нулевые векторы и отбрасывать их. Количество векторов, выдаваемых алгоритмом, будет равно размерности подпространства, порождённого векторами (т.е. количеству линейно независимых векторов, которые можно выделить среди исходных векторов). | ||
+ | |||
+ | Процесс Грама ― Шмидта может быть истолкован как разложение невырожденной квадратной матрицы в произведение ортогональной (или унитарной матрицы в случае эрмитова пространства) и верхнетреугольной матрицы с положительными диагональными элементами ― QR-разложение, что есть частный случай разложения Ивасавы. | ||
== Программная реализация алгоритма == | == Программная реализация алгоритма == | ||
=== Особенности реализации последовательного алгоритма === | === Особенности реализации последовательного алгоритма === | ||
=== Локальность данных и вычислений === | === Локальность данных и вычислений === | ||
− | === Возможные способы и особенности параллельной реализации алгоритма === | + | === Возможные способы и особенности параллельной реализации алгоритма=== |
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
+ | Проведём исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. Исследовалась собственная параллельная реализация алгоритма, написанная на языке C++ с использованием стандарта MPI. | ||
+ | |||
+ | Версия MPI - intel opemmpi/1.8.4-icc | ||
+ | |||
+ | Модель используемого CPU - Intel Xeon X5570 2.93GHz | ||
+ | |||
+ | Набор и границы значений изменяемых реализации алгоритма: | ||
+ | |||
+ | * число процессоров [10:100] с шагом 10; | ||
+ | * размерность векторов [50 : 500] с шагом 50. | ||
+ | |||
+ | На следующем рисунке приведен график времени работы алгоритма в зависимости от числа процессов и размерности векторов. | ||
+ | [[File:gram-schmidt.jpg|мини|upright=3|centre|framecentre|Рис. 6 Время работы алгоритма в зависимости от числа процессоров и размерности задачи]] | ||
+ | |||
+ | Так как алгоритм не поддается полному распараллеливанию, то при увеличении числа процессоров не происходит значительного ускорения. Это связано с тем, что параллельно можно вычислить лишь <math>b_j</math>. То есть после каждого вычисления <math>b_j</math> можно запускать параллельное вычисление <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> для всех <math>i=j+1..N</math>. Кроме того можно параллельно производить вычитание соответствующего <math>\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}</math> для всех <math>i=j+1..N</math>. Но после этого все процессоры должны синхронизоваться и дождаться когда получат следующий вектор <math>b_j</math>. | ||
+ | |||
+ | Проведём еще одно исследование исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта для меньшего числа процессоров. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. Исследовалась параллельная реализация алгоритма, написанная с использованием стандарта MPI. | ||
+ | |||
+ | Набор и границы значений изменяемых реализации алгоритма: | ||
+ | |||
+ | * число процессоров [1:10] с шагом 1; | ||
+ | * размерность векторов [20 : 400] с шагом 20. | ||
+ | На графике 1 приведено время работы программы для разных значений размерности матрицы и числа процессоров. | ||
+ | [[File:GramTime.jpg|мини|upright=3|centre|framecentre|1|Рис. 7 Время работы программы в зависимости от числа процессоров и размерности задачи]] | ||
+ | На графике 2 приведена производительность работы программы для разных значений размерности матрицы и числа процессоров. | ||
+ | [[File:GramPerf.jpg|мини|upright=3|centre|framecentre|2|Рис. 8 Производительность программы в зависимости от количество процессоров и размерности задачи]] | ||
+ | На графике 3 приведена эффективность распараллеливания работы программы для разных значений размерности матрицы и числа процессоров. | ||
+ | [[File:GramParEffi.jpg|мини|upright=3|centre|framecentre|3|Рис. 9 Эффективность распараллеливания алгоритма]] | ||
+ | Оценим значения масштабируемости данной реализации: | ||
+ | Общей оценкой масштабируемости по числу процессов будем называть среднее значение вклада каждого элемента рассматриваемой области значений в общую оценку по числу процессов по всем элементам рассматриваемой области параметров запуска. | ||
+ | * По числу процессоров: -0.009265. С увеличением числа процессоров наблюдается слабое снижение эффективности работы параллельной реализации алгоритма. | ||
+ | Общей оценкой масштабируемости по размеру решаемой задачи будем называть среднее значение вклада каждого элемента рассматриваемой области значений в общую оценку по размеру решаемой задачи по всем элементам рассматриваемой области параметров запуска. | ||
+ | * По размеру задачи: 0.00002555. При увеличении размерности задачи эффективность алгоритма очень медленно возрастает. | ||
+ | Общей оценкой масштабируемости будем называть средний размер вклада каждого элемента рассматриваемой области значений в общую оценку по всем элементам рассматриваемой области параметров запуска. | ||
+ | * По двум направлениям: 0.002079. | ||
+ | [https://bitbucket.org/d_2007_90/gram-schmidt-process/src Реализация алгоритма] | ||
+ | |||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
=== Выводы для классов архитектур === | === Выводы для классов архитектур === | ||
=== Существующие реализации алгоритма === | === Существующие реализации алгоритма === | ||
+ | |||
+ | Примеры библиотек, содержащих реализации алгоритма ортогонализации Грама-Шмидта: | ||
+ | # [http://numerics.mathdotnet.com/Matrix.html C#] | ||
+ | # [http://www.mathworks.com/help/symbolic/mupad_ref/linalg-orthog.html?searchHighlight=orthogonalization Matlab] | ||
+ | # [https://reference.wolfram.com/language/ref/Orthogonalize.html Wolfram] | ||
+ | # [http://www.nag.co.uk/numeric/Fl/manual/pdf/F05/f05aaf.pdf NAG Fortran] | ||
+ | |||
+ | ====Реализация для пакета Mathematica==== | ||
+ | Данный скрипт, предназначенный для пакета Mathematica, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в фигурных скобках последней строки. Количество векторов и их координат могут быть произвольными. В данном случае для примера взяты векторы {-2,1,0}, {-2,0,1}, {-0.5,-1,1}. | ||
+ | |||
+ | <source lang='mathematica'> | ||
+ | Projection[v1_, v2_] := (v1.v2*v2)/v2.v2 | ||
+ | MultipleProjection[v1_, vecs_] := Plus @@ (Projection[v1, #1] &) /@ vecs | ||
+ | GramSchmidt[mat_] := Fold[Join[#1, {#2 - MultipleProjection[#2, #1]}] &, {}, mat] | ||
+ | GramSchmidt[{{-2, 1, 0}, {-2, 0, 1}, {-0.5, -1, 1}}] | ||
+ | </source> | ||
+ | ====Реализация для Maxima==== | ||
+ | Данный скрипт, предназначенный для пакета Maxima, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в квадратных скобках. Количество векторов и их координат могут быть произвольными. | ||
+ | |||
+ | В данном случае для примера взяты векторы [-2,1,0], [-2,0,1], [-0.5,-1,1]. | ||
+ | <source lang="tsql"> | ||
+ | load(eigen); | ||
+ | x: matrix ([-2,1,0],[-2,0,1],[-0.5,-1,1]); | ||
+ | y: gramschmidt(x); | ||
+ | </source> | ||
+ | ====Реализация для MATLAB==== | ||
+ | В MATLAB существует встроенная функция [Q,R]=qr(A) которая осуществляет ортогонализацию Грамма-Шмидта столбцов матрицы А. При этом столбцы матрицы предполагаются линейно независимыми. | ||
+ | |||
+ | [Q, R] = qr(A) возвращает матрицы Q с ортонормированными столбцами и обратимую верхнетреугольную матрицу R, так что A=Q*R | ||
== Литература == | == Литература == | ||
+ | * Беклемишев Д. В. Курс аналитической геометрии и линейной алгебры. — М.: Наука. | ||
+ | * Гельфанд И. М. Лекции по линейной алгебре. — 5-е. — М.: Добросвет, МЦНМО, 1998. — 319 с. — ISBN 5-7913-0015-8. | ||
+ | *Ильин В. А., Позняк Э. Г. Линейная алгебра. 6-е изд. — М.: Физматлит, 2010. — 280 с. — ISBN 978-5-9221-0481-4. | ||
+ | *Халмош П. Конечномерные векторные пространства = Finite-Dimensional Vector Spaces. — М.: Физматгиз, 1963. — 263 с. | ||
+ | *Фаддеев Д. К. Лекции по алгебре. — 5-е. — СПб.: Лань, 2007. — 416 с. | ||
+ | |||
+ | |||
+ | == Ссылки == | ||
+ | * https://ru.wikipedia.org/wiki/Процесс_Грама_―_Шмидта ''Процесс Грама ― Шмидта'' |
Текущая версия на 15:52, 23 ноября 2016
Эта работа прошла предварительную проверку Дата последней правки страницы: 23.11.2016 Данная работа соответствует формальным критериям. Проверено ASA. |
Ортогонализация Грама-Шмидта | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(N^2k)[/math] |
Объём входных данных | [math]Nk[/math] |
Объём выходных данных | [math]Nk[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(N)[/math] |
Ширина ярусно-параллельной формы | [math]O(Nk)[/math] |
Основные авторы описания: Инжелевская Дарья Валерьевна
Содержание
- 1 Свойства и структура алгоритмов
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
- 4 Ссылки
1 Свойства и структура алгоритмов
1.1 Общее описание алгоритма
Создатели алгоритма Эрхард Шмидт и Йорген Педерсен Грам.
Эрхард Шмидт (14 января 1876 — 6 декабря 1959) — немецкий математик, с 1917 года профессор Берлинского университета. В 1946—58 первый директор Института математики АН ГДР. Основные труды по теории функций, интегральным уравнениям, функциональному анализу. Определил и изучил геометрически гильбертово пространство, используя аналогии с геометрией Евклида. Занимался квадратичными формами. Известен оператор Гильберта-Шмидта и процесс Грама ― Шмидта.
Йорген Педерсен Грам (27 июня 1850 - 29 апреля 1916) датский математик, родился в герцогстве Шлезвиг, Дания и умер в Копенгагене, Дания. Основные его работы "On series expansions determined by the methods of least squares", "Investigations of the number of primes less than a given number". Математический метод, который носит его имя, процесс Грама-Шмидта, впервые был опубликован в 1883. Теорема Грама и Определитель Грама также назван в его честь.
Процесс Грама-Шмидта — это один из алгоритмов, в которых на основе множества линейно независимых векторов [math]{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}[/math] строится множество ортогональных векторов [math]{\displaystyle \mathbf {b}_{1},\;\ldots ,\;\mathbf {b} _{N}} [/math] или ортонормированных векторов [math]{\displaystyle \mathbf {e} _{1},\;\ldots ,\;\mathbf {e}_{N}} [/math], причём так, что каждый вектор [math]{\displaystyle \mathbf {b} _{j}} [/math] или [math]{\displaystyle \mathbf {e} _{j}}[/math] может быть выражен линейной комбинацией векторов [math]{\displaystyle \mathbf {a} _{1},\;\ldots ,\; \mathbf {a} _{j}}[/math].
1.2 Математическое описание алгоритма
1.2.1 Классический процесс Грама — Шмидта
1.2.1.1 Алгоритм
Пусть имеются линейно независимые векторы [math]\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N[/math].
Определим оператор проекции следующим образом: [math]\mathbf{proj}_{\mathbf{b}}\,\mathbf{a} = {\langle \mathbf{a}, \mathbf{b} \rangle \over \langle \mathbf{b}, \mathbf{b}\rangle} \mathbf{b} ,[/math]
где [math]\langle \mathbf{a}, \mathbf{b} \rangle[/math] — скалярное произведение векторов [math]\mathbf{a}[/math] и [math]\mathbf{b}[/math].
Скалярное произведение для двух векторов [math]\mathbf{ a= [a_1, a_2, ...,a_k]}[/math] и [math]\mathbf{ b= [b_1, b_2, ..., b_k]}[/math] в k-мерном действительном пространстве определяется как:
- [math]\langle \mathbf{a}, \mathbf{b} \rangle=\sum_{i=1}^k a_ib_i=a_1b_1+a_2b_2+\cdots+ a_kb_k[/math].
Этот оператор проецирует вектор [math]\mathbf{a}[/math] коллинеарно вектору [math]\mathbf{b}[/math].
Ортогональность векторов [math]\mathbf{a}[/math] и [math]\mathbf{b}[/math] достигается на шаге (2).
Классический процесс Грама — Шмидта выполняется следующим образом:
- [math] {\begin{array}{lclr} {\mathbf {b}}_{1}&=&{\mathbf {a}}_{1}&(1)\\ {\mathbf {b}}_{2}&=&{\mathbf {a}}_{2}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{2}&(2)\\ {\mathbf {b}}_{3}&=&{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{1}}}\,{\mathbf {a}}_{3}-{\mathbf {proj}}_{{{\mathbf {b}}_{2}}}\,{\mathbf {a}}_{3}&(3)\\ {\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}&(4)\\ &\vdots &&\\{\mathbf {b}}_{N}&=&{\mathbf {a}}_{N}-\displaystyle \sum _{{j=1}}^{{N-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{N}&(N) \end{array}} [/math]
На основе каждого вектора [math]\mathbf{b}_j \;(j = 1 \ldots N)[/math] может быть получен нормированный вектор: [math]\mathbf{e}_j = {\mathbf{b}_j\over \| \mathbf{b}_j \|}[/math] (у нормированного вектора направление будет таким же, как у исходного, а длина — единичной). Причем берется норма согласованная со скалярным произведением [math]\| x \| = \sqrt{\langle x, x \rangle}[/math]
Результаты процесса Грама — Шмидта:
[math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N[/math] — система ортогональных векторов либо
[math]\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N[/math] — система ортонормированных векторов.
Вычисление [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N[/math] носит название ортогонализации Грама — Шмидта, а [math]\mathbf{e}_1,\;\ldots,\;\mathbf{e}_N[/math] — ортонормализации Грама — Шмидта.
1.2.1.2 Доказательство
Докажем ортогональность векторов [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_N[/math].
Для этого вычислим скалярное произведение [math]\langle \mathbf{b}_1, \mathbf{b}_2 \rangle[/math], подставив в него формулу (2). Мы получим ноль. Равенство нулю скалярного произведения векторов означает, что эти векторы ортогональны. Затем вычислим скалярное произведение [math]\langle \mathbf{b}_1, \mathbf{b}_3 \rangle[/math], используя результат для [math]\langle \mathbf{b}_1, \mathbf{b}_2 \rangle[/math] и формулу (3). Мы снова получим ноль, то есть векторы [math]\mathbf{b}_1[/math] и [math]\mathbf{b}_3[/math] ортогональны. Общее доказательство выполняется методом математической индукции.
1.2.1.3 Геометрическая интерпретация — вариант 1
Рассмотрим формулу (2) — второй шаг алгоритма. Её геометрическое представление изображено на рис. 1:
1 — получение проекции вектора [math]\mathbf{a}_2[/math] на [math]\mathbf{b}_1[/math];
2 — вычисление [math]\mathbf{a}_2-\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_2[/math], то есть перпендикуляра, которым выполняется проецирование конца [math]\mathbf{a}_2[/math] на [math]\mathbf{b}_1[/math].
Этот перпендикуляр — вычисляемый в формуле (2) вектор [math]\mathbf{b}_2[/math];
3 — перемещение полученного на шаге 2 вектора [math]\mathbf{b}_2[/math] в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (2).
На рисунке видно, что вектор [math]\mathbf{b}_2[/math] ортогонален вектору [math]\mathbf{b}_1[/math], так как [math]\mathbf{b}_2[/math] является перпендикуляром, по которому [math]\mathbf{a}_2[/math] проецируется на [math]\mathbf{b}_1[/math].
Рассмотрим формулу (3) — третий шаг алгоритма — в следующем варианте:
- [math] \begin{array}{lcr} \mathbf{b}_3 = \mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3+\mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3) & & (6) \\ \end{array} [/math]
Её геометрическое представление изображено на рис. 2:
1 — получение проекции вектора [math]\mathbf{a}_3[/math] на [math]\mathbf{b}_1[/math];
2 — получение проекции вектора [math]\mathbf{a}_3[/math] на [math]\mathbf{b}_2[/math];
3 — вычисление суммы [math]\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3[/math], то есть проекции вектора [math]\mathbf{a}_3[/math] на плоскость, образуемую векторами [math]\mathbf{b}_1[/math] и [math]\mathbf{b}_2[/math]. Эта плоскость закрашена на рисунке серым цветом;
4 — вычисление [math]\mathbf{a}_3-(\mathbf{proj}_{\mathbf{b}_1}\mathbf{a}_3 + \mathbf{proj}_{\mathbf{b}_2}\mathbf{a}_3)[/math], то есть перпендикуляра, которым выполняется проецирование конца [math]\mathbf{a}_3[/math] на плоскость, образуемую векторами [math]\mathbf{b}_1[/math] и [math]\mathbf{b}_2[/math]. Этот перпендикуляр — вычисляемый в формуле (6) вектор [math]\mathbf{b}_3[/math];
5 — перемещение полученного [math]\mathbf{b}_3[/math] в начало координат. Это перемещение сделано на рисунке лишь для наглядности. Оно не является математическим действием и поэтому не отражается в формуле (6).
На рисунке видно, что вектор [math]\mathbf{b}_3[/math] ортогонален векторам [math]\mathbf{b}_1[/math] и [math]\mathbf{b}_2[/math], так как [math]\mathbf{b}_3[/math] является перпендикуляром, по которому [math]\mathbf{a}_3[/math] проецируется на плоскость, образуемую векторами [math]\mathbf{b}_1[/math] и [math]\mathbf{b}_2[/math].
Таким образом, в процессе Грама — Шмидта для вычисления [math]\mathbf{b}_j[/math] выполняется проецирование [math]\mathbf{a}_j[/math] ортогонально на гиперплоскость?, формируемую векторами [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}[/math]. Вектор [math]\mathbf{b}_j[/math] затем вычисляется как разность между [math]\mathbf{a}_j[/math] и его проекцией. То есть [math]\mathbf{b}_j[/math] — это перпендикуляр? от конца [math]\mathbf{a}_j[/math] к гиперплоскости?, формируемой векторами [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{j-1}[/math]. Поэтому [math]\mathbf{b}_j[/math] ортогонален векторам, образующим эту гиперплоскость?.
1.2.1.4 Численная неустойчивость
При вычислении на ЭВМ по формулам (1) — (5) векторы [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_{N-1}[/math] часто не точно ортогональны из-за ошибок округления. Из-за потери ортогональности в процессе вычислений классический процесс Грама — Шмидта называют численно неустойчивым.
1.2.2 Модифицированный процесс Грама — Шмидта
Процесс Грама — Шмидта может быть сделан более вычислительно устойчивым путём небольшой модификации. Вместо вычисления [math]\mathbf{b}_j[/math] как: [math] \mathbf{b}_j = \mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_2}\,\mathbf{a}_j - \ldots - \mathbf{proj}_{\mathbf{b}_{j-1}}\,\mathbf{a}_j \; (7) [/math] этот вектор вычисляется следующим образом:
- [math] \begin{array}{lclr} \mathbf{a}_j^{(1)} & = &\mathbf{a}_j - \mathbf{proj}_{\mathbf{b}_1}\,\mathbf{a}_j & (8) \\ \mathbf{a}_j^{(2)} & = &\mathbf{a}_j^{(1)} - \mathbf{proj}_{\mathbf{b}_2} \, \mathbf{a}_j^{(1)} & (9) \\& \vdots & & \\ \mathbf{a}_j^{(j-2)} & = & \mathbf{a}_j^{(j-3)} - \mathbf{proj}_{\mathbf{b}_{j-2}} \, \mathbf{a}_j^{(j-3)} & (10) \\ \mathbf{b}_j & = & \mathbf{a}_j^{(j-2)} - \mathbf{proj}_{\mathbf{b}_{j-1}} \, \mathbf{a}_j^{(j-2)} & (11) \end{array} [/math]
1.3 Вычислительное ядро алгоритма
Вычислительное ядро последовательной версии метода ортогонализации Грамма-Шмидта можно составить из множественных (всего их [math]\frac{N(N-1)}{2}[/math]) вычислений проекций : [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math]
1.4 Макроструктура алгоритма
Данный алгоритм использует в качестве составных частей другие алгоритмы. В дальнейшем имеет смысл описывать алгоритм не в максимально детализированном виде (т.е. на уровне арифметических операций), а давать только его макроструктуру, то есть описывает структуру и состав макроопераций. Типичной макрооперацией, часто встречающиеся в алгоритме является оператор проекции векторов. Как записано и в описании ядра алгоритма, основную часть метода ортогонализации Грамма-Шмидта составляют множественные (всего их [math]\frac{N(N-1)}{2}[/math]) вычисления оператора проекции.
1.5 Схема реализации последовательного алгоритма
Следующий алгоритм реализует нормализацию Грамма-Шмидта. Векторы [math]v_1,...,v_k[/math] заменяются набором ортонормированных векторов, которые имеют ту же линейную оболочку.
Вычислительная сложность этого [math]2Nk^2[/math] операции с плавающей точкой, где N - размерность векторов.
Последовательность исполнения метода следующая:
1. [math]\mathbf {b}_{1}=\mathbf {a}_{1}[/math]
2. Далее для всех векторов [math]\mathbf {b}_{i}[/math] для [math]i=2 ... N[/math] производим вычисление по следующей формуле: [math]{\mathbf {b}}_{i}={\mathbf {a}}_{i}-\displaystyle \sum _{{j=1}}^{{i-1}}{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}[/math].
В ней на каждом шаге [math]i[/math] по очереди вычисляются все [math]{\mathbf {proj}}_{{{\mathbf {b}}_{j}}}\,{\mathbf {a}}_{i}[/math] для [math]j=1 ... i-1[/math]
Пример реализации на Python. Функция работает для произвольного количества векторов любой размерности. При этом если количество векторов [math]\mathbf{a}_1,\;\ldots,\;\mathbf{a}_N[/math] больше их размерности или они линейно зависимы то функция возвращает максимально возможное число линейно независимых векторов [math]\mathbf{b}_1,\;\ldots,\;\mathbf{b}_n[/math], а остальные векторы [math]\mathbf{b}_{n+1},\;\ldots,\;\mathbf{b}_N[/math] нулевые.
import random
def GramSchmidt(*a):
k=len(a[0])
N=len(a);
b = [[0] * k for i in range(N)]
b[0]=a[0]
for i in range(1,N):
sum=a[i]
for j in range(0,i):
scolar_ab=0
scolar_bb=0
proj=[i for i in range(k)]
for n in range(k):
scolar_ab+=b[j][n]*a[i][n]
scolar_bb+=b[j][n]*b[j][n]
for n in range(k):
proj[n]=(scolar_ab/scolar_bb)*b[j][n]
for n in range(k):
sum[n]-=proj[n]
b[i]=sum
return b;
l1=[random.randrange(0,10) for i in range(3)]
l2=[random.randrange(0,10) for i in range(3)]
l3=[random.randrange(0,10) for i in range(3)]
print(l1,l2,l3)
print(GramSchmidt(l1,l2,l3))
1.6 Последовательная сложность алгоритма
Для построения ортогонального набора векторов в последовательном варианте требуется:
При k=N:
- [math]\frac{N(N-1)}{2}[/math] делений,
- [math]\frac{N(N-1)(3N-1)}{2}[/math] сложений (вычитаний),
- [math]\frac{3N^2 (N-1)}{2}[/math] умножений.
При k≠N
- [math]\frac{N(N-1)}{2}[/math] делений,
- [math]\frac{N(N-1)(3k-1)}{2}[/math] сложений (вычитаний),
- [math]\frac{3Nk (N-1)}{2}[/math] умножений.
При классификации по последовательной сложности, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
Опишем граф алгоритма как аналитически, так и в виде рисунка.
Граф алгоритма состоит из двух групп вершин, расположенных в целочисленных узлах двух областей.
Первая группа вершин расположена в двумерной области, соответствующая ей операция [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i} = {\langle \mathbf{a_i}, \mathbf{b_j} \rangle \over \langle \mathbf{b_j}, \mathbf{b_j}\rangle} \mathbf{b_j} ,[/math] .
Естественно введённые координаты области таковы:
i — меняется в диапазоне от 2 до N, принимая все целочисленные значения;
j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения.
Аргументы операции следующие:
[math]a_i[/math]: элементы входных данных, а именно [math]a_i[/math];
[math]b_j[/math]: результат срабатывания операции, соответствующей вершине из второй группы, с координатой j.
Вторая группа вершин расположена в двухмерной области, соответствующая ей операция [math]a - b[/math].
Естественно введённые координаты области таковы:
i — меняется в диапазоне от 2 до N, принимая все целочисленные значения;
j — меняется в диапазоне от 1 до i-1, принимая все целочисленные значения.
Аргументы операции следующие:
a:
j=1 - входные данные [math]a_j[/math]
j>1 - результат срабатывания операции, соответствующей вершине из второй группы, с координатой j-1
b:
результат срабатывания операции, соответствующей вершине из первой группы, с координатой j
1.8 Ресурс параллелизма алгоритма
В параллельном варианте, в отличие от последовательного, вычисления [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math] можно выполнять параллельно. То есть после каждого вычисления [math]b_j[/math] можно запускать параллельное вычисление [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math] для всех [math]i=j+1..N[/math]. Кроме того можно параллельно производить вычитание соответствующего [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math] для всех [math]i=j+1..N[/math].
Для ортогонализации в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]N[/math] ярусов в каждом из которых будет выполнятся,
- по [math]N - 1[/math] операторов проекций (в каждом из ярусов линейное количество проекций, в зависимости от яруса — от [math]N-1[/math] до [math]1[/math]),
- по [math]N - 1[/math] сложений (в каждом из ярусов линейное количество сложений, в зависимости от яруса — от [math]N-1[/math] до [math]1[/math]),
При классификации по высоте (количество ярусов в ЯПФ ) ЯПФ, таким образом, ортогонализация Грамма-Шмидта относится к алгоритмам со сложностью [math]O(N)[/math].
При классификации по ширине (максимальное количество вершин в ярусе) ЯПФ его сложность будет [math]O(N^2)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: множества линейно независимых векторов [math]{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}[/math], каждый из которых описывается [math]\mathbf{ a_i= [a^i_1, a^i_2, ..., a^i_k]}[/math] .
Дополнительные ограничения:
- вектора [math]{\displaystyle \mathbf {a} _{1},\;\ldots ,\;\mathbf {a} _{N}}[/math] линейно независимы, поэтому [math]k \geqslant N[/math].
Объём входных данных: [math]Nk[/math]
Выходные данные: множество ортогональных векторов [math]{\displaystyle \mathbf {b} _{1},\;\ldots ,\;\mathbf {b} _{N}} [/math], каждый из которых описывается [math]\mathbf{ b_i= [b^i_1, b^i_2, ..., b^i_k]}[/math] .
Объём выходных данных: [math]Nk[/math]
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение кубической к квадратичной).
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, также линейна.
Алгоритм почти полностью детерминирован — единственность результата выполнения гарантирована, однако возможно накопление ошибок округления при использовании классического процесса Грама-Шмидта. При этом модифицированный алгоритм дает более хороший результат.
Алгоритм является численно неустойчивым — ошибки округления могут привести к неортогональности полученных векторов.
Процесс Грама — Шмидта может применяться также к бесконечной последовательности линейно независимых векторов.
Кроме того, процесс Грама — Шмидта может применяться к линейно зависимым векторам. В этом случае он выдаёт [math]\mathbf{0}[/math] (нулевой вектор) на шаге [math]j[/math], если [math]\mathbf{a}_j[/math] является линейной комбинацией векторов [math]\mathbf{a}_1,\;\ldots,\;\mathbf{a}_{j-1}[/math]. Если это может случиться, то для сохранения ортогональности выходных векторов и для предотвращения деления на ноль при ортонормировании алгоритм должен делать проверку на нулевые векторы и отбрасывать их. Количество векторов, выдаваемых алгоритмом, будет равно размерности подпространства, порождённого векторами (т.е. количеству линейно независимых векторов, которые можно выделить среди исходных векторов).
Процесс Грама ― Шмидта может быть истолкован как разложение невырожденной квадратной матрицы в произведение ортогональной (или унитарной матрицы в случае эрмитова пространства) и верхнетреугольной матрицы с положительными диагональными элементами ― QR-разложение, что есть частный случай разложения Ивасавы.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Проведём исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Исследовалась собственная параллельная реализация алгоритма, написанная на языке C++ с использованием стандарта MPI.
Версия MPI - intel opemmpi/1.8.4-icc
Модель используемого CPU - Intel Xeon X5570 2.93GHz
Набор и границы значений изменяемых реализации алгоритма:
- число процессоров [10:100] с шагом 10;
- размерность векторов [50 : 500] с шагом 50.
На следующем рисунке приведен график времени работы алгоритма в зависимости от числа процессов и размерности векторов.
Так как алгоритм не поддается полному распараллеливанию, то при увеличении числа процессоров не происходит значительного ускорения. Это связано с тем, что параллельно можно вычислить лишь [math]b_j[/math]. То есть после каждого вычисления [math]b_j[/math] можно запускать параллельное вычисление [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math] для всех [math]i=j+1..N[/math]. Кроме того можно параллельно производить вычитание соответствующего [math]\mathbf{proj}_{\mathbf{b_j}}\,\mathbf{a_i}[/math] для всех [math]i=j+1..N[/math]. Но после этого все процессоры должны синхронизоваться и дождаться когда получат следующий вектор [math]b_j[/math].
Проведём еще одно исследование исследование масштабируемости параллельной реализации процесса ортогонализации Грама-Шмидта для меньшего числа процессоров. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Исследовалась параллельная реализация алгоритма, написанная с использованием стандарта MPI.
Набор и границы значений изменяемых реализации алгоритма:
- число процессоров [1:10] с шагом 1;
- размерность векторов [20 : 400] с шагом 20.
На графике 1 приведено время работы программы для разных значений размерности матрицы и числа процессоров.
На графике 2 приведена производительность работы программы для разных значений размерности матрицы и числа процессоров.
На графике 3 приведена эффективность распараллеливания работы программы для разных значений размерности матрицы и числа процессоров.
Оценим значения масштабируемости данной реализации: Общей оценкой масштабируемости по числу процессов будем называть среднее значение вклада каждого элемента рассматриваемой области значений в общую оценку по числу процессов по всем элементам рассматриваемой области параметров запуска.
- По числу процессоров: -0.009265. С увеличением числа процессоров наблюдается слабое снижение эффективности работы параллельной реализации алгоритма.
Общей оценкой масштабируемости по размеру решаемой задачи будем называть среднее значение вклада каждого элемента рассматриваемой области значений в общую оценку по размеру решаемой задачи по всем элементам рассматриваемой области параметров запуска.
- По размеру задачи: 0.00002555. При увеличении размерности задачи эффективность алгоритма очень медленно возрастает.
Общей оценкой масштабируемости будем называть средний размер вклада каждого элемента рассматриваемой области значений в общую оценку по всем элементам рассматриваемой области параметров запуска.
- По двум направлениям: 0.002079.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Примеры библиотек, содержащих реализации алгоритма ортогонализации Грама-Шмидта:
2.7.1 Реализация для пакета Mathematica
Данный скрипт, предназначенный для пакета Mathematica, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в фигурных скобках последней строки. Количество векторов и их координат могут быть произвольными. В данном случае для примера взяты векторы {-2,1,0}, {-2,0,1}, {-0.5,-1,1}.
Projection[v1_, v2_] := (v1.v2*v2)/v2.v2
MultipleProjection[v1_, vecs_] := Plus @@ (Projection[v1, #1] &) /@ vecs
GramSchmidt[mat_] := Fold[Join[#1, {#2 - MultipleProjection[#2, #1]}] &, {}, mat]
GramSchmidt[{{-2, 1, 0}, {-2, 0, 1}, {-0.5, -1, 1}}]
2.7.2 Реализация для Maxima
Данный скрипт, предназначенный для пакета Maxima, проводит процесс ортогонализации Грама ― Шмидта над векторами, заданными в квадратных скобках. Количество векторов и их координат могут быть произвольными.
В данном случае для примера взяты векторы [-2,1,0], [-2,0,1], [-0.5,-1,1].
load(eigen);
x: matrix ([-2,1,0],[-2,0,1],[-0.5,-1,1]);
y: gramschmidt(x);
2.7.3 Реализация для MATLAB
В MATLAB существует встроенная функция [Q,R]=qr(A) которая осуществляет ортогонализацию Грамма-Шмидта столбцов матрицы А. При этом столбцы матрицы предполагаются линейно независимыми.
[Q, R] = qr(A) возвращает матрицы Q с ортонормированными столбцами и обратимую верхнетреугольную матрицу R, так что A=Q*R
3 Литература
- Беклемишев Д. В. Курс аналитической геометрии и линейной алгебры. — М.: Наука.
- Гельфанд И. М. Лекции по линейной алгебре. — 5-е. — М.: Добросвет, МЦНМО, 1998. — 319 с. — ISBN 5-7913-0015-8.
- Ильин В. А., Позняк Э. Г. Линейная алгебра. 6-е изд. — М.: Физматлит, 2010. — 280 с. — ISBN 978-5-9221-0481-4.
- Халмош П. Конечномерные векторные пространства = Finite-Dimensional Vector Spaces. — М.: Физматгиз, 1963. — 263 с.
- Фаддеев Д. К. Лекции по алгебре. — 5-е. — СПб.: Лань, 2007. — 416 с.
4 Ссылки
- https://ru.wikipedia.org/wiki/Процесс_Грама_―_Шмидта Процесс Грама ― Шмидта