【NA】基于QR分解的特征值迭代法

Francis于1961-1962年利用矩阵的QR分解建立了计算矩阵特征值的QR方法，是计算中小型矩阵全部特征值的最有效方法之一。

• 本篇的主线是第一部分介绍QR分解，第二部分介绍从QR分解引出的特征值QR迭代算法，第三部分讨论QR迭代法的收敛性，第四部分引用UTEP-Math 5330中基于Householder变换的QR分解实现，第五部分做总结以及更多讨论。

QR分解.

• 【引理】设 $\vec{x}=(x_1,x_2,\cdots ,x_i,\cdots ,x_j,\cdots ,x_n{)}^T$，其中 $x_i\cdot x_j\ne 0$，则存在一个Givens变换矩阵 $P_{ij}$ 使得$$P_{ij}\vec{x}=(y_1,y_2,\cdots ,y_i,\cdots ,y_j,\cdots ,y_n{)}^T$$其中$$y_i=\sqrt{x_i^2+x_j^2}$$$$y_j=0$$
• 上述Givens矩阵 $P_{ij}$ 中$$P[i,i]=\cos \theta ,P[i,j]=\sin \theta$$$$P[j,i]=-\sin \theta ,P_{[}j,j]=\cos \theta$$$$\cos \theta =\frac{x_i}{\sqrt{x_i^2+x_j^2}},\sin \theta =\frac{x_j}{\sqrt{x_i^2+x_j^2}}$$
• 实际操作时，可以通过对 $\vec{x}$ 的规范化操作来避免数值溢出。

---

• 【定理】对于 $n$ 阶非奇异矩阵 $A$ 而言，可以分解为正交阵 $Q$ 和上三角阵 $R$ 的乘积，即 $A=QR.$ 该分解是唯一的，当矩阵 $R$ 中的元素满足$$\forall i\in [1,n],R[i,i]>0$$

---

• 【定理】对于 $n$ 阶非奇异矩阵 $A$ 而言，存在一系列平面旋转变换矩阵 $P_k|k=1,2,\cdots ,n-1$，使得 $$P_{n-1}\cdots P_2{P}_{1}A=R$$其中 $R$ 是上三角阵且 $r_{ii}>0,i=1,2,\cdots ,n.$
• 【证明】由于 $A$ 非奇异，所以不可能存在全零列，因此对于第一列，可以做出一系列平面旋转变换 $$P_1=\prod P_{1k}$$使得第一列约化为 $(r_{11},0,\cdots ,0{)}^T$，同理对于后面的第 $m$列，可以做一系列平面旋转变换 $P_k$ 使得该列约化为 $(r_{1m},r_{2m},\cdots ,r_{mm},0,\cdots ,0{)}^T$，最终得到上三角阵，并引入对角阵 D = { ± 1 , ± 1 , ⋯   , ± 1 } D=\{\pm1,\pm1,\cdots,\pm1\} D={±1,±1,⋯,±1} 使得主对角元全正。
• 显然如果令 $$Q^T=\prod _{i=1}^{n-1}{P}_i$$那么$$Q^{T}A=R$$因为 $Q$ 是正交阵，所以 $Q^T=Q^{-1}$，所以有$$A=(Q^T{)}^{-1}R=QR$$
• 具体地如何进行QR分解，得到正交阵 $Q$ 和上三角阵 $R$，主流的有三种方案：Schmidt正交化、Givens变换和Householder变换，在第四部分对前两者有简介。

QR迭代算法.

• 设 $A$ 是 $n$ 阶非奇异阵，记 $A_1=A$，对其进行QR分解得到 $$A_1=Q_1{R}_1$$按照如下公式进行迭代：$$A_{k+1}=R_k{Q}_k=Q_k^T{A}_k{Q}_k=Q_{k+1}{R}_{k+1}\text{\hspace{1em}(*)}$$
• 上述利用矩阵QR分解和递推式 $(\ast )$ 构造矩阵序列 { A k } \{A_k\} {Ak​} 的过程即为QR算法。

---

• $A_{k+1}=Q_k^T{A}_k{Q}_k$
• $A_{k+1}=(Q_1{Q}_2\cdots Q_k{)}^T{A}_1(Q_1{Q}_2\cdots Q_k)$
• 若记 $$\widetilde{Q_k}=\prod _{i=1}^k{Q}_i,\widetilde{R_k}=\prod _{i=k}^1{R}_i$$那么矩阵 $A^k$ 的QR分解式为 $$A^k=\widetilde{Q_k}\widetilde{R_k}$$
• 【证明】数学归纳法
  ①当 $k=1$ 时显然有 $$A=Q_1{R}_1=\widetilde{Q_1}\widetilde{R_1}$$
  ②设 $A^{k-1}$ 分解式为 $$A^{k-1}=\widetilde{Q_{k-1}}\widetilde{R_{k-1}}$$
  ③那么 $$\begin{array}{rl}\widetilde{Q_k}\widetilde{R_k}& =Q_1{Q}_2\cdots Q_{k-1}(Q_k{R}_k)R_{k-1}\cdots R_2{R}_1\\ & =Q_1{Q}_2\cdots Q_{k-1}{A}_k{R}_{k-1}\cdots R_2{R}_1\\ & =\widetilde{Q_{k-1}}{A}_k\widetilde{R_{k-1}}\\ & =\widetilde{Q_{k-1}}(\widetilde{Q_{k-1}}{)}^{T}A\widetilde{Q_{k-1}}\widetilde{R_{k-1}}\\ & =A{A}^{k-1}\\ & =A^k\end{array}$$

收敛性.

• 对于 $n$ 阶非奇异阵 $A$，若 $A$ 的特征值满足 $①|{\lambda }_1|>|{\lambda }_2|>\cdots >|{\lambda }_n|$；$②A$ 有标准型 $A=XD{X}^{-1}$，其中 $D=diag{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$，且 $X^{-1}$ 有三角分解 $X^{-1}=LU$，则由QR算法产生的矩阵序列 { A k } \{A_k\} {Ak​} 有下列极限：$$\underset{k\to \infty }{\lim }{A}_k=R$$ 其中 $R$ 为上三角阵，并且 $r_{ii}={\lambda }_i.$
• 【证明】对于非奇异矩阵 $A$，存在可逆矩阵 $X$ 使得 $X^{-1}AX=D$，即 $A=XD{X}^{-1}$，因此有$$A^k=X{D}^k{X}^{-1}\text{\hspace{1em}(1)}$$条件中认为 $X^{-1}$ 存在LU分解 $X^{-1}=LU$，因此 $(1)$ 式可写为 $$A^k=X{D}^{k}LU=X(D^{k}L{D}^{-k})D^{k}U\text{\hspace{1em}(2)}$$将 $(2)$ 中 $D^{k}L{D}^{-k}$ 做加性分解，得到$$D^{k}L{D}^{-k}=E+{\Lambda }_k$$其中 ${\Lambda }_k$ 是一个下三角阵，且其主对角元 $\Lambda [i,i]=0$，非零元素表达式为$$\Lambda [i,j]=L[i,j]\cdot (\frac{{\lambda }_i}{{\lambda }_j}{)}^k,i>j$$由于特征值按绝对值降序排列，因此有下面的极限存在$$\underset{k\to \infty }{\lim }\Lambda =0$$这里的 $0$ 代表零矩阵。由于矩阵 $X$ 非奇异，对其进行QR分解，有$$X=QR$$至此 $(2)$ 式可以写为$$A^k=QR(E+{\Lambda }_k)D^{k}U=Q(E+R{\Lambda }_k{R}^{-1})R{D}^{k}U\text{\hspace{1em}(3)}$$当 $k\to \infty$ 时，矩阵 $R(E+{\Lambda }_k)$ 可逆，从而 $R(E+\Lambda )R^{-1}=E+R{\Lambda }_k{R}^{-1}$ 可逆，对其进行QR分解，有$$E+R{\Lambda }_k{R}^{-1}=Q_k^{\prime }{R}_k^{\prime }$$从而 $(3)$ 式写为$$A^k=(Q{Q}_k^{\prime })(R_k^{\prime }R{D}^{k}U)\text{\hspace{1em}(4)}$$显然 $(4)$ 式已经是 $A$ 的QR分解，为确保唯一性，引入对角阵 D ∗ = d i a g { ± 1 , ± 1 , ⋯   , ± 1 } D^*=diag\{\pm1,\pm1,\cdots,\pm1\} D∗=diag{±1,±1,⋯,±1}将 $(4)$ 写为$$A^k=(Q{Q}_k^{\prime }{D}^{\ast })(D^{\ast }{R}_k^{\prime }R{D}^{k}U)\text{\hspace{1em}(5)}$$可得$$\widetilde{Q_k}=Q{Q}_k^{\prime }{D}^{\ast }$$$$\widetilde{R_k}=D^{\ast }{R}_k^{\prime }R{D}^{k}U$$因此 $$A_{k+1}=(Q_1{Q}_2\cdots Q_k{)}^T{A}_1(Q_1{Q}_2\cdots Q_k)=(\widetilde{Q_k}{)}^{T}A\widetilde{Q_k}=D^{\ast }(Q_k^{\prime }{)}^T{Q}^{T}AQ{Q}_k^{\prime }{D}^{\ast }\text{\hspace{1em}(6)}$$因为 $A=XD{X}^{-1},X=QR$，所以 $Q^{T}AQ=RD{R}^{-1}$，因此 $(6)$ 式可以写为$$A_{k+1}=D^{\ast }(Q_k^{\prime }{)}^T(RD{R}^{-1})Q_k^{\prime }{D}^{\ast }\text{\hspace{1em}(7)}$$记 $$R_0=RD{R}^{-1},g_k=Q_k^{\prime }{D}^{\ast }$$$$R_0[i,i]={\lambda }_i$$所以 $(7)$ 式可以写为$$A_{k+1}=g_k^T{R}_0{g}_k\text{\hspace{1em}(8)}$$由于$$(E+{\Lambda }_k\to E)\Rightarrow (R(E+\Lambda )R^{-1}=E+R{\Lambda }_k{R}^{-1}\to E)$$所以 $Q_k^{\prime }\to E$，所以 g k → D ∗ = d i a g { ± 1 , ± 1 , ⋯   , ± 1 } g_k\rightarrow D^*=diag\{\pm1,\pm1,\cdots,\pm1\} gk​→D∗=diag{±1,±1,⋯,±1}所以 $k\to \infty$ 时，$$\lim A_{k+1}[i,i]={\lambda }_i,i=1,2,\cdots ,n$$其收敛速度由 $||{\Lambda }_k|{|}_{\infty }$ 决定，有上限$$||{\Lambda }_k|{|}_{\infty }\le C\cdot \underset{j\in [1,n-1]}{\max }|\frac{{\lambda }_{j+1}}{{\lambda }_j}{|}^k$$
• 上述收敛速度定义式 $r_n=|{\lambda }_n/{\lambda }_{n-1}|$ 与幂法中的收敛速度定义式相似度极高，当 $r_n\approx 1$ 时收敛极慢，可以效仿幂法中采取原点位移策略，即对于矩阵 $A-s\cdot E$ 使用QR分解，那么收敛速度即为 $$r_n^{\ast }=\left|\frac{{\lambda }_n-s}{{\lambda }_{n-1}-s}\right|$$上述算法称为带原点位移的QR迭代算法。

【Reference】UTEP-Math 5330.

• 逻辑极度清晰，可惜我弄丢了上述Lecture Note的网址，关于Householder变换可以参考《Householder变换》。
• 另外注意到Note中讨论了 $a_{11}=\sqrt{a_{11}^2+a_{21}^2+\cdots +a_{n1}^2}=0$ 的情况，不难发现上述等式成立意味着 $a_{i1}=0,i=1,2,\cdots ,n$，那么 $A$ 是奇异矩阵，无法进行QR分解。
• 第一个重点是 Algorithm3.1 的伪代码，一般情况下我们假定 $A$ 是 $n$ 阶非奇异阵，即 $m=n$，伪代码中 $\beta$ 是当次迭代列向量的模值，$\gamma$ 是向量 $v$ 的分母，后续计算出向量 $v$ 以及对应的豪氏矩阵 $H_v=E-2v{v}^T$与矩阵 $A$ 进行左乘 $A_1=H_v{A}_0$，变量 $Q$ 初值是 $E_m$ 单位阵，迭代过程中对 $H_v$ 左乘，所以最终停止时 $$Q=\prod _{i=1}^{n-1}{H}_i$$即为QR分解中 $A=QR$ 的 【$Q$】，而此时保存 $A$ 的单元被 【$R$】所覆盖。
• 11-13行代码作用是在矩阵本身已是上三角化时直接进行下一次迭代。这篇Note的一个特点是，开篇已经假定 $A$ 是非奇异阵，但后续的分析和伪码又是针对 $A\in R^{m\times n}$ 考虑的，个人猜测是出于算法健壮性。
• Note的最后一部分给出了基于Householder变化进行QR分解复杂度分析：
  
  
• 上半部分关于 $A$ 阵迭代无需显式计算豪氏矩阵以及 $Q$ 阵迭代的浮点运算 flops 挺靠谱的，后半部分关于整个算法的复杂度分析持保留态度，尤其是最后给出的 $\frac{2}{3}{r}^3$，这里不妨认为 $$r=\underset{m=n}{\min }(m,n)=n$$那么复杂度为 $\frac{2}{3}{r}^3$，但这个已经出现 $4p^2$ 项，求和后应为 $\frac{4}{3}{r}^3$，个人倾向于这个答案，大致的 $O(n^3)$ 是无疑的。
• 给出结果为 $\frac{4}{3}{n}^3$ 的文献如下，来自UTexas(这俩名字应该是一个学校吧)：
  
  
  

Summary.

• 需要明确的是，第四部分内容是对于QR分解以及基于 Householder变换实现QR分解的，并不是用于求解矩阵全部特征值的QR迭代法。
• 【存疑】QR迭代法思路之一就是在迭代公式中使用Householder变换实现每一次迭代的QR分解，使用Householder变换进行QR分解的时间复杂度为 $\frac{2}{3}{n}^3$，大约是第一部分构造性证明中使用Givens变换的 $\frac{2}{3}$ 倍。
• 另一种QR迭代法的思路是首先使用Householder变换，将一般实矩阵(实对称矩阵)约化为上海森伯格矩阵(对称三对角矩阵)，而后再进行QR迭代，这样做的优势在于可以使用对于稀疏矩阵较为高效的Givens变换阵(常数时间即可计算出)实现QR分解，而无需使用需要 $n^2$ 才能计算出的豪斯荷尔德矩阵。

---

• 【亟待解决】直接使用Householder变换进行的QR分解时间复杂度，有文献表明是 $\frac{4}{3}{n}^3$，也有文献认为是 $\frac{2}{3}{n}^3$，个人倾向于前一种说法，理由在上面。
• 另外的两种QR分解方法，格拉姆-施密特正交化GS(以及修正格拉姆-施密特正交化MGS)和吉文斯变换，其复杂度也为 $O(n^3)$，有说法认为前者的常数因子为 $1$，后者为 $\frac{4}{3}.$
• 使用豪斯荷尔德变换将实矩阵约化为上Hessenberg矩阵需要的乘法次数：$\frac{5}{3}{n}^3.$
• 那么下图中的说法就不能成立：