线性方程组与本征值问题(06)：Householder 变换和 Hessenberg 约化

Householder 变换

\(\newcommand{\R}{\mathbb{R}}\newcommand{\b}{\mathbf}\newcommand{\bi}{\boldsymbol}\newcommand{\mdl}[1]{\left|\!\left|#1\right|\!\right|}\)考虑一个矢量\(\bi{v}\in\R\), 构造一个反射矩阵 \[ \b{P}_{\boldsymbol{v}}=\b{I} - \frac{2\bi{v}\bi{v}^\b{T}}{\mdl{\bi{v}}^2} \] 之所以称作反射矩阵, 是因为它和右侧任何矢量的积为 \[ \b{P}_\bi{v}\bi{x}=\bi{x}-\frac{2\bi{v}^\b{T}\cdot\bi{x}}{\bi{v}^\b{T}\cdot\bi{v}}\bi{v} \] 它在\(\bi{v}\)的方向上减去了\(\bi{x}\)分量的两倍, 相当于将其反向; 而正交于\(\bi{v}\)的任何方向上不发生变化. 这个变换称作 Householder 变换. 应当指出的是

1. Householder 变换不改变矢量的模, 因此 Householder 矩阵是一个正交矩阵; 同时它还是对称矩阵.
2. 用矩阵计算变换的时空复杂度通常为\(O(n^2)\), 但 Householder 变换实际上只需要存储\(\bi{v}\), 计算也简化为矢量内积, 复杂度为\(O(n)\).

现在考虑一种特殊情况: \(\bi{x}\)和\(\bi{v}\)相当接近, 以至于只有一个方向上的区别 \[ \bi{v}=\bi{x}\pm\mdl{\bi{x}}\bi{e}_m\qquad\Rightarrow\qquad\b{P}_\bi{v}\bi{x}=\bi{x}-\bi{v}=\mp\mdl{\bi{x}}\bi{e}_m \]

这样的Householder变换将原矢量转动到这个方向上.

同理, 构造矩阵 \[ \b{P}_{(k)}=\begin{pmatrix} \b{I}_{k} & \b{O} \\ \b{O} & \b{P}_\bi{v} \end{pmatrix}\qquad\text{其中, }\bi{v}= \bi{x}_{(n-k)} + \mdl{\bi{x}_{(n-k)}}\bi{e}_{k+1} \] 所谓的\(\bi{x}_{(n-k)}\)指\(\bi{x}\)的后\(n-k\)维构成的矢量. 这样, 有 \[ \b{P}_{(k)}\bi{x}= \begin{pmatrix} \b{I}_{k} & \b{O} \\ \b{O} & \b{P}_\bi{v} \end{pmatrix} \cdot\begin{pmatrix} \bi{x}^{(k)} \\ \bi{x}_{(n-k)} \end{pmatrix} =\begin{pmatrix} \bi{x}^{(k)} \\ \mdl{\bi{x}_{(n-k)}} \\ 0 \\ \vdots \end{pmatrix} \] 其实它还是 Householder 变换, 把 Householder 矢量的维数补全即可 \[ \bi{u}=\begin{pmatrix}\b{0}\\\bi{v}\end{pmatrix}\Rightarrow\b{P}_{(k)}=\b{P}_\bi{u} \] 这种变换保持矢量的前\(k\)维不变, 后\(n-k\)维整体转动到第\(k+1\)个基矢量上.

我们更感兴趣的是以它为转动矩阵的相似变换: \[ \b{P}_{(k)}^\b{T}\b{A}\b{P}_{(k)}=\begin{pmatrix} \b{A}_{11} & \b{A}_{12}\b{P}_\bi{v} \\ \b{P}_\bi{v}\b{A}_{21} & \b{P}_\bi{v}\b{A}_{22}\b{P}_\bi{v} \end{pmatrix} \] 假设对于任意\(k\), \(\b{P}_\bi{v}\)的选取都以\(\b{A}_{21}\)的最后一个列矢量为\(\bi{x}\); 进一步假设\(\b{A}_{21}\)的其它列均为\(0\), 则 \[ \b{P}_\bi{v}\b{A}_{21}=\begin{pmatrix} \b{O} & \begin{matrix} \mdl{\bi{x}} \\ 0 \\ \vdots \\ 0 \\ 0 \end{matrix} \end{pmatrix} \]

Hessenberg 约化

按照这种构造方式, 可以将任意矩阵逐步约化.

1. 对原始矩阵用\(\b{P}_{(1)}\)转动, 则第\(1\)列从第\(3\)行开始变为\(0\);
2. 对更新后的矩阵用\(\b{P}_{(2)}\)转动, 则第\(2\)列从第\(4\)行开始变为\(0\);
3. ...
4. 对更新后的矩阵, 用\(\b{P}_{(n-2)}\)转动, 则第\(n-2\)列只让第\(n\)行变为\(0\).

下面给出\(n=4\)矩阵的约化过程, 每次把\(\b{A}_{21}\)框出: \[ \begin{pmatrix} * & \begin{matrix}* & * & *\end{matrix} \\ \boxed{\begin{matrix} * \\ * \\ * \end{matrix}} & \begin{matrix}* & * & * \\ * & * & * \\ * & * & *\end{matrix} \end{pmatrix}\overset{\b{P}_{(1)}旋转}{\longrightarrow} \begin{pmatrix} \begin{matrix}* & * \\ * & *\end{matrix} & \begin{matrix}* & * \\ * & *\end{matrix} \\ \boxed{\begin{matrix}0 & * \\ 0 & *\end{matrix}} & \begin{matrix}* & * \\ * & *\end{matrix} \end{pmatrix}\overset{\b{P}_{(2)}旋转}{\longrightarrow} \begin{pmatrix} * & * & * & * \\ * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \end{pmatrix} \] 按照这种方式, 得到矩阵的特点为: 下对角线以下的元都必定为\(0\), 其它元不定. 这种矩阵叫做 上 Hessenberg 矩阵, 约化方案称作 Hessenberg-Householder 约化.

计算量

1. 每轮约化需要得到该列对应的矢量\(\bi{v}\), 它需要计算列矢量的模, \(n\)次乘法
2. 这时还可以仅仅存储矢量, 因为\(\b{P}_{(k)}=\b{I}-\displaystyle\frac{2\bi{u}\bi{u}^\b{T}}{\mdl{\bi{u}}^2}\), \(\bi{u}\)是包含了前\(k\)维的零元的完整矢量. 从而给每列计算变换结果时, 是\(n-k\)次乘法; 它从第\(k\)列计算到第\(n\)列, 因此每轮的计算量约为\((n-k+1)(n-k)\).

可见, 总计算量是\(O(n^3)\). 具体地, 求和给出 \[ \sum_{k=1}^{n-2}(n-k+1)(n-k)=\frac{(n-2)(n^2+2n+3)}{3}\qquad(n\ge3) \] 大约算了\(n^3/3\)次, 复杂度与一轮 QR 分解相当. 看起来是做无用功, 实则约化为上 Hessenberg 形式对迭代非常有利, 具体怎样有利于计算是下一节将要讨论的内容.