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

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

Граф алгоритма без отображения входных и выходных данных. n=4. Зелёным выделены операции типа a+bb, тёмно-синим и сиреневым - типа a+bc, красным - вычисление [math]\gamma , v_{i}[/math], жёлтым - завершение вычисления [math]\frac{(x,v)}{\gamma}[/math]

Как уже сказано в описании ядра, основная часть - вычисления на каждом шагу скалярных произведений [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.5 Схема реализации последовательного алгоритма

Последовательность выполнения алгоритма обычно записывается как последовательное "обнуление" поддиагональных элементов столбцов, начиная с 1-го столбца и заканчивая предпоследним [math](n-1)[/math]-м.

При этом в каждом "обнуляемом" [math]i[/math]-м столбце "обнуляются" сразу все его поддиагональные элементы одновременно, с [math](i+1)[/math]-го до [math]n[/math]-го.

Каждое "обнуление" [math]i[/math]-го столбца состоит из двух шагов: а) вычисление параметров матрицы отражения [math]U_{i}[/math] такой, чтобы при умножении на неё слева "обнулились" все поддиагональные его элементы; б) умножение слева матрицы отражения [math]U_{i}[/math] на текущую версию матрицы.

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

В последовательной версии основная сложность алгоритма определяется прежде всего вычислениями скалярных произведений векторов, а также модификаций векторов вида [math]x'=x-\alpha v[/math], причем над векторами убывающей по ходу алгоритма размерности. Они, если не учитывать возможную разреженность, составляют (в главном члене) по [math]2n^3/3[/math] операций действительного умножения и сложения/вычитания.

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

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

На рисунке приведён граф алгоритма метода Хаусхолдера в наиболее его быстром (с параллельной точки зрения) варианте, использующем то, что с точностью до множителя ведущий вектор матрицы отражения отличается отличается от подстолбца, где выполняется очередное исключение, только одним элементом.

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

Для понимания ресурса параллелизма в разложении матрицы порядка [math]n[/math] методом Хаусхолдера нужно рассмотреть критический путь графа.

Как видно из описания разных вершин, вычисления при "обнулении" [math]i[/math]-го столбца параметров отражения и скалярных произведений состоят из основной части - ветви длиной по [math]n-i[/math] умножений и сложений - и коррекции вычислений, которые составляют [math]O(1)[/math] операций.

Поэтому по грубой (без членов низших порядков) оценке критический путь метода Хаусхолдера будет идти через [math]\frac{n^2}{2}[/math] умножений и [math]\frac{n^2}{2}[/math] сложений/вычитаний.

Поэтому в параллельном варианте, как и в последовательном, основную долю требуемого для выполнения алгоритма времени будут определять операции вида [math]a+bc[/math].

При классификации по высоте ЯПФ, таким образом, метод Хаусхолдера относится к алгоритмам с квадратичной сложностью. При классификации по ширине ЯПФ его сложность будет также квадратичной (без расширения ряда ярусов, связанных с векторными операциями сложения, пришлось бы увеличить вдвое длину критического пути).

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

Входные данные: плотная квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]).

Объём входных данных: [math]n^2[/math].

Выходные данные: правая треугольная матрица [math]R[/math] (ненулевые элементы [math]r_{ij}[/math] в последовательном варианте хранятся в элементах исходной матрицы [math]a_{ij}[/math]), унитарная (ортогональная) матрица Q - как произведение матриц Хаусхолдера (отражения) (их вектора нормалей к плоскостям отражения в последовательном варианте хранятся в поддиагональных элементах исходной матрицы [math]a_{ij}[/math] и в одном дополнительном столбце размерности n).

Объём выходных данных: [math]n^2+n[/math].

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

Соотношение последовательной и параллельной сложности, как хорошо видно, является линейным, что даёт определённый стимул для распараллеливания. Однако у наискорейшей ЯПФ ширина квадратична, что указывает на дисбаланс между загруженностями устройств при попытке её реально запрограммировать. Поэтому более практично даже при хорошей (быстрой) вычислительной сети оставить количество устройств (например, узлов кластера) линейным по размеру матрицы, что удвоит критический путь реализуемой ЯПФ.

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

Алгоритм в рамках выбранной версии полностью детерминирован.

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

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

2.1 Особенности реализации последовательного алгоритма

В варианте с кратчайшим критическим путём графа алгоритма (с использованием зависимости между обнуляемым вектором и направляющим вектором отражения) метод Хаусхолдера (отражений) QR-разложения квадратной вещественной матрицы на Фортране 77 можно записать так:

	DO  I = 1, N-1
            DO K = I, N
              SX(K)=A(N,I)*A(N,K)
            END DO
            DO J = N-1, I, -1
               DO K = I, N
                  SX(K)=SX(K)+A(J,I)*A(J,K)
               END DO
            END DO
            IF (SX(I).NE.0.) THEN
               R = SQRT (SX(I))
               IF (A(I,I).NE.0.) THEN
                   
               ELSE

               END IF               
            ELSE
               S1(I)=1.
               SX(I)=.5
               DO K = I+1, N
                  A(K,I)=-A(K,I)
               END DO
            END IF
		
	END DO

В этом варианте R расположена в верхнем правом треугольнике массива A,