线性方程组与本征值问题(07)：非标准的 QR 算法
\(\newcommand{\R}{\mathbb{R}}\newcommand{\b}{\mathbf}\newcommand{\bi}{\boldsymbol}\newcommand{\mdl}[1]{\left|\!\left|#1\right|\!\right|}\newcommand{\abs}[1]{\left|#1\right|}\)前一篇笔记给出了任意方矩阵的 Hessenberg-Householder 约化方法, 即如何把任意方矩阵变换为上 Hessenberg 矩阵; 这一节中讨论如何对上 Hessenberg 矩阵进行 QR 分解, 它比普通方阵具有什么样的优势, 并引入shifted-QR 算法解决非对角元收敛过慢的问题.
Givens 变换矩阵
请原谅我一直引入新概念, 但这对后面的讨论是有好处的, 所谓 Givens 变换矩阵, 其实就是在第\(j,k\)维度的子空间上做一个类似\(x,y\)平面上的转动. \[ \text{e.g.}\qquad\b{G}(1,2,\theta)=\begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \\ & & \b{I}_{n-2} \end{pmatrix} \] 一般形式比较繁琐, 按下不表. 不难知道, 将 Givens 矩阵左乘上任意方矩阵做变换, 只有第\(j,k\)两行发生改变; 右乘, 则只有\(j,k\)两列改变.
Hessenberg-QR 分解
给定一个上 Hessenberg 矩阵\(\b{H}\)(具体形式如下), 并且做一个 Givens 转动:
\[ \begin{pmatrix} * & * & * & \cdots & * & * & * \\ \boxed{*} & * & * & \cdots & * & * & * \\ 0 & * & * & \cdots & * & * & * \\ 0 & 0 & * & \cdots & * & * & * \\ 0 & 0 & 0 & * & * & * & * \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & * & * \end{pmatrix} \overset{\b{G}^\b{T}(1,2,\theta_1)}{\longrightarrow} \begin{pmatrix} * & * & * & \cdots & * & * & * \\ \boxed{0} & * & * & \cdots & * & * & * \\ 0 & * & * & \cdots & * & * & * \\ 0 & 0 & * & \cdots & * & * & * \\ 0 & 0 & 0 & * & * & * & * \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & * & * \end{pmatrix} \]
注意: Givens 矩阵显然会修改前两行的所有非零值, 我们设定一个转角\(\theta_1\), 使得框选出的元恰好变为零元, 其它元未必变成零元.
接着, 再在第\(2,3\)行做一个 Givens 转动, 同样设定一个转角\(\theta_2\), 恰好把框选出的下对角元变为零元:
\[ \begin{pmatrix} * & * & * & \cdots & * & * & * \\ 0 & * & * & \cdots & * & * & * \\ 0 & \boxed{*} & * & \cdots & * & * & * \\ 0 & 0 & * & \cdots & * & * & * \\ 0 & 0 & 0 & * & * & * & * \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & * & * \end{pmatrix} \overset{\b{G}^\b{T}(2,3,\theta_2)}{\longrightarrow} \begin{pmatrix} * & * & * & \cdots & * & * & * \\ 0 & * & * & \cdots & * & * & * \\ 0 & \boxed{0} & * & \cdots & * & * & * \\ 0 & 0 & * & \cdots & * & * & * \\ 0 & 0 & 0 & * & * & * & * \\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0 & * & * \end{pmatrix} \]
这样, 我们一共经过\(n-1\)个 Givens 转动, 上 Hessenberg 矩阵被变换为上三角矩阵:
\[ \b{U}=\left(\prod_{j=1}^{n-1}\b{G}^\b{T}(n-j,n-j+1,\theta_{n-j})\right)\b{H} \]
然后再对\(\b{U}\)右乘一系列 Givens 矩阵, 得到新的矩阵:
\[ \b{H}^{(1)}=\b{U}\prod_{j=1}^{n-1}\b{G}(j,j+1,\theta_j) \]
不难验证, \(\b{H}^{(1)}\)仍然是上 Hessenberg 矩阵. 注意到 Givens 矩阵是正交的, 因此上面两个公式就是对\(\b{H}\)的 QR 分解. 把一个矩阵\(\b{A}\)先进行 Hessenberg 约化得到\(\b{H}\), 再像这样作 QR 分解, 称之为 Hessenberg-QR 分解.
计算量
确定转角比较简单, 关键是变换过程. 每乘上一个 Givens 矩阵需要改变两行, 即:
- 第\(1,2\)行, 计算\(2n\)个矩阵元
- 第\(2,3\)行, 计算\(2n-2\)个矩阵元
- ...
- 第\(n-1,n\)行, 计算\(4\)个矩阵元
共计算\((n+2)(n-1)\)个矩阵元, 每个矩阵元需要两次乘法, 一次加法, 计算量为\(\lambda(n+2)(n-1)\), \(2<\lambda<3\).
又因为 Givens 矩阵有左乘和右乘两次乘法, 所以计算量约为\(4n^2\sim6n^2\)(按照乘法次数计算). 这比标准的 QR 分解少多了. 要知道, 以得到本征谱或者再加上全体本征矢为目的时, 分解原始矩阵\(\b{A}\)和分解\(\b{H}\)没什么区别, 两者本征值一样, 本征矢相差一个总的旋转矩阵. 因此在迭代次数比较多时, Hessenberg-QR 分解比原始分解方法节约了计算量.
Hessenberg-QR 分解的局限性
这实际上也是(标准)QR分解的局限性. 第5节最后给出的定理提到, 矩阵非对角元的在 QR 分解下的收敛速度为 \[ \abs{t_{ij}^{(k)}}\sim\abs{\frac{\lambda_i}{\lambda_j}}^k, \abs{\lambda_i} \le \abs{\lambda_j} \] 注意到它是可以取等的, 因此, 如果具有一对互为相反数的实本征值, 则它们的交叉项无法收敛到\(0\); 如果模相近, 则交叉项收敛非常慢, 这都导致最终无法得到实 Schur 形式——它要求所有实本征值都以对角元而非对角块的形式出现.
shifted-QR 分解
为了解决本征值模相近或相等时给迭代带来的困难, 一个比较自然的想法就是给矩阵加上对角项的漂移: \[ \begin{align} \b{Q}^{(k)}\b{R}^{(k)}&=\b{T}^{(k-1)}-\mu\b{I} \\ \b{T}^{(k)}&=\b{R}^{(k)}\b{Q}^{(k)}+\mu\b{I} \end{align} \] 这样, 实际参与分解的本征值为\(\{\lambda_1-\mu,\lambda_2-\mu,\cdots,\lambda_n-\mu\}\), 分解完毕后又把本征值还原.
显然, 这时仍须满足本征值的不等式: \(\abs{\lambda_1-\mu}\ge\abs{\lambda_2-\mu}\ge\cdots\ge\abs{\lambda_n-\mu}\), 注意到每次分解可以采用不同的\(\mu\), 考虑采用\(\mu^{(k)}=t_{n,n}^{(k)}\), 则\(\lambda_n-\mu\)的模始终很小, 这会让\(t_{n,n-1}\)快速趋于\(0\).