Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы, вещественный точечный вариант

QR-разложение методом Хаусхолдера (отражений)
Последовательный алгоритм
Последовательная сложность	[math]\frac{4 n^3}{3}[/math]
Объём входных данных	[math]n^2[/math]
Объём выходных данных	[math]n(n + 1)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(n^2)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/math]

Основные авторы описания: А.В.Фролов

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

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

Метод Хаусхолдера (в советской математической литературе называется также методом отражений) используется для разложения матриц в виде [math]A=QR[/math] ([math]Q[/math] - унитарная, [math]R[/math] — правая треугольная матрица). При этом матрица [math]Q[/math] хранится и используется не в своём явном виде, а в виде произведения матриц отражения. Каждая из матриц отражения может быть определена одним вектором. Это позволяет в классическом исполнении метода отражений хранить результаты разложения на месте матрицы A с использованием минимального одномерного дополнительного массива.

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

В методе Хаусхолдера для выполнения [math]QR[/math]-разложения матрицы используются умножения её слева на матрицы Хаусхолдера (отражений).

Матрица отражений (Хаусхолдера) - матрица вида [math]U=E-2ww^*[/math], где [math]w[/math] - вектор, удовлетворяющий равенству [math]w^{*}w=1[/math]. Является одновременно унитарной ([math]U^{*}U=E[/math]) и эрмитовой ([math]U^{*}=U[/math]), поэтому обратна самой себе ([math]U^{-1}=U[/math]).

На [math]i[/math]-м шаге метода с помощью преобразования отражения "убираются" ненулевые поддиагональные элементы в [math]i[/math]-м столбце. Таким образом, после [math]n-1[/math] шагов преобразований получается матрица [math]R[/math] из [math]QR[/math]-разложения.

На каждом из шагов метода матрицу отражений обычно представляют не в стандартном виде, а в виде [math]A=E-\frac{1}{\gamma}vv^*[/math], где [math]v[/math] находится через координаты текущего [math]i[/math]-го столбца так:

[math]s[/math] - вектор размерности [math]n+1-i[/math], составленный из элементов [math]i[/math]-го столбца, начиная с [math]i[/math]-го.

Если [math](s,s)=0[/math], то [math]v=e_{i}[/math], [math]\gamma = \frac{1}{2}[/math].

В остальных случаях вычисляем [math]u = \frac{1}{\sqrt{(s,s)}}s[/math], и далее [math]v_{j}=0[/math] при [math]j\lt i[/math], [math]v_{j}=u_{j-i+1}[/math] при [math]j\gt i[/math], а [math]v_{i}=1[/math], если [math]u_{1}=0[/math] и [math]v_{i}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|)[/math] для остальных значений. При этом [math]\gamma = 1+|u_{1}|[/math].

После вычисления вектора [math]v[/math] подстолбцы справа от ведущего модифицируются по формулам [math]x'=x-\frac{(x,v)}{\gamma}v[/math].

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

Основную часть алгоритма составляют вычисления на каждом шагу скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстолбцов [math]x[/math]справа от текущего, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math]. Это используется при программировании метода во многих библиотеках для его конструирования из стандартных подпрограмм (например, из BLAS).

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

Как уже сказано в описании ядра, основная часть - вычисления на каждом шагу скалярных произведений [math](s,s)[/math] и [math](x,v)[/math] для всех подстолбцов [math]x[/math]справа от текущего, а также векторные операции [math]x'=x-\frac{(x,v)}{\gamma}v[/math]. При этом, однако, строгая последовательность выполнения этих трёх подшагов не обязательна, в силу связи получаемых векторов [math]s[/math] и [math]v[/math] можно одновременно с [math](s,s)[/math] вычислять и произведения [math](x,s)[/math] с последующим выражением через них [math](x,v)[/math]. Это позволяет почти вдвое уменьшать критический путь графа алгоритма.

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

Граф алгоритма без отображения входных и выходных данных. n=4.