Метод Хаусхолдера (отражений) приведения матрицы к двухдиагональной форме
Свойства и структура алгоритма
Общее описание алгоритма
Метод Хаусхолдера (в советской математической литературе чаще называется методом отражений) используется для приведения симметричных вещественных матриц к двухдиагональному виду, или, что то же самое, для разложения [math]\displaystyle{ A=QDU^T }[/math] ([math]\displaystyle{ Q, U }[/math] - ортогональные, [math]\displaystyle{ D }[/math] — нижняя двухдиагональная матрица)[1]. При этом матрицы [math]\displaystyle{ Q, U }[/math] хранятся и используются не в своём явном виде, а в виде произведения матриц отражения[2]. Каждая из матриц отражения может быть определена одним вектором. Это позволяет в классическом исполнении метода отражений хранить результаты разложения на месте матрицы A с использованием двух одномерных дополнительных массивов.
В данной статье рассматривается именно классическое исполнение, в котором не используются приёмы типа сдваивания при вычислениях скалярных произведений.
Математическое описание алгоритма
В методе Хаусхолдера для выполнения разложения матрицы в произведение двухдиагональной и двух ортогональных используются попеременные умножения слева и справа её текущих модификаций на матрицы Хаусхолдера (отражений).
Матрица отражений (Хаусхолдера) - матрица вида [math]\displaystyle{ U=E-2ww^* }[/math], где [math]\displaystyle{ w }[/math] - вектор, удовлетворяющий равенству [math]\displaystyle{ w^{*}w=1 }[/math]. Является одновременно унитарной ([math]\displaystyle{ U^{*}U=E }[/math]) и эрмитовой ([math]\displaystyle{ U^{*}=U }[/math]), поэтому обратна самой себе ([math]\displaystyle{ U^{-1}=U }[/math]).
На [math]\displaystyle{ i }[/math]-м шаге метода с помощью преобразования отражения "убираются" ненулевые поддиагональные элементы в [math]\displaystyle{ i }[/math]-м столбце, а потом (кроме шага с номером [math]\displaystyle{ i=n-1 }[/math]) с помощью преобразования отражения "убираются" ненулевые наддиагональные элементы в [math]\displaystyle{ i }[/math]-м столбце, кроме самого левого из них. Таким образом, после [math]\displaystyle{ n-1 }[/math] шагов преобразований получается хессенбергова [math]\displaystyle{ D }[/math] из разложения матрицы в произведение хессенберговой и двух ортогональных.
На каждом из шагов метода матрицы отражений обычно представляют не в стандартном виде, а в виде [math]\displaystyle{ U=E-\frac{1}{\gamma}vv^* }[/math], где [math]\displaystyle{ v }[/math] находится для левых матриц отражения через координаты текущего [math]\displaystyle{ i }[/math]-го столбца так:
[math]\displaystyle{ s }[/math] - вектор размерности [math]\displaystyle{ n+1-i }[/math], составленный из элементов [math]\displaystyle{ i }[/math]-го столбца, начиная с [math]\displaystyle{ i }[/math]-го.
Если [math]\displaystyle{ (s,s)=0 }[/math], то [math]\displaystyle{ v=e_{i} }[/math], [math]\displaystyle{ \gamma = \frac{1}{2} }[/math].
В остальных случаях по алгоритму вычисляется [math]\displaystyle{ u = \frac{1}{\sqrt{(s,s)}}s }[/math], и далее [math]\displaystyle{ v_{j}=0 }[/math] при [math]\displaystyle{ j\lt i }[/math], [math]\displaystyle{ v_{j}=u_{j-i+1} }[/math] при [math]\displaystyle{ j\gt i }[/math], а [math]\displaystyle{ v_{i}=1 }[/math], если [math]\displaystyle{ u_{1}=0 }[/math] и [math]\displaystyle{ v_{i}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|) }[/math] для остальных значений. При этом [math]\displaystyle{ \gamma = 1+|u_{1}| }[/math].
После вычисления вектора [math]\displaystyle{ v }[/math] подстолбцы справа от ведущего модифицируются по формулам [math]\displaystyle{ x'=x-\frac{(x,v)}{\gamma}v }[/math].
Аналогично для правых матриц отражений [math]\displaystyle{ U=E-\frac{1}{\gamma}vv^* }[/math] [math]\displaystyle{ v }[/math] находится через координаты текущей [math]\displaystyle{ i }[/math]-й строки так:
[math]\displaystyle{ s }[/math] - вектор размерности [math]\displaystyle{ n-i }[/math], составленный из элементов [math]\displaystyle{ i }[/math]-й строки, начиная с [math]\displaystyle{ i+1 }[/math]-го элемента.
Если [math]\displaystyle{ (s,s)=0 }[/math], то [math]\displaystyle{ v=e_{i+1} }[/math], [math]\displaystyle{ \gamma = \frac{1}{2} }[/math].
В остальных случаях по алгоритму вычисляется [math]\displaystyle{ u = \frac{1}{\sqrt{(s,s)}}s }[/math], и далее [math]\displaystyle{ v_{j}=0 }[/math] при [math]\displaystyle{ j\lt i+1 }[/math], [math]\displaystyle{ v_{j}=u_{j-i+2} }[/math] при [math]\displaystyle{ j\gt i+1 }[/math], а [math]\displaystyle{ v_{i+1}=1 }[/math], если [math]\displaystyle{ u_{1}=0 }[/math] и [math]\displaystyle{ v_{i+1}=\frac{u_{1}}{|u_{1}|}(1+|u_{1}|) }[/math] для остальных значений. При этом [math]\displaystyle{ \gamma = 1+|u_{1}| }[/math].