Метод Хаусхолдера (отражений) приведения матрицы к двухдиагональной форме

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

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

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

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

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

В данной статье рассматривается именно классическое исполнение, в котором не используются приёмы типа сдваивания при вычислениях скалярных произведений.

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

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

Матрица отражений (Хаусхолдера) - матрица вида [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]i=n-1[/math]) с помощью преобразования отражения "убираются" ненулевые наддиагональные элементы в [math]i[/math]-м столбце, кроме самого левого из них. Таким образом, после [math]n-1[/math] шагов преобразований получается правая двудиагональная [math]D[/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}|=|v_{i}|[/math].

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

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

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

Если [math](s,s)=0[/math], то [math]v=e_{i+1}[/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+1[/math], [math]v_{j}=u_{j-i+2}[/math] при [math]j\gt i+1[/math], а [math]v_{i+1}=1[/math], если [math]u_{1}=0[/math] и [math]v_{i+1}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|)[/math] для остальных значений. При этом [math]\gamma = 1+|u_{1}|=|v_{i+1}|[/math].

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

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

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

Всё приведение состоит из [math]2n-3[/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]-м, попеременно с последовательным "обнулением" наддиагональных (кроме первой кодиагонали) элементов строк, начиная с 1й строки и заканчивая [math](n-2)[/math]-й

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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