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

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}|[/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}|[/math].