[矩阵计算]求解矩阵特征值的QR方法

更新： 28 JUL 2016

求矩阵特征值问题的重要性自不待言，在此仅举一例，即在量子力学中解$\textbf{HC}=E\textbf{C}$特征方程求本征值（能量）和本征矢。本人目的就在于此，求矩阵特征值是一个总体思路。然而实际上哈密顿量$\textbf{H}$通常是一个巨大的稀疏矩阵，采用最经典的QR方法实际上不可行。此处介绍QR方法作为后文介绍Lanczos方法的知识铺垫。

参考教材：徐树方：矩阵计算的理论与方法

 

一、秩1矩阵的特征值问题

秩1矩阵$A$能够表示为

$$ A=uv^T $$

其中$u$，$v$是$\mathbb{C}^n$中的非零向量。显然$A$是$n\times n$的矩阵，但是只有一个非零的特征值$v^Tu$。

现取$\mathbb{C}^n$中任意一非零向量$x$，做

$$ y = Ax $$

若$y = 0$则说明$x$是矩阵$A$对应于特征值0的特征向量。若$y \neq 0$则$y$是对应于特征值$v^Tu$的特征向量，因为

$$ Ay=AAx=uv^Tuv^Tx=u(v^Tu)v^Tx=(v^Tu)uv^Tx=(v^Tu)y. $$

综上，求秩1矩阵的特征向量可以在线性空间中任取一组正交基，通过如上的乘法即可得到全部的特征值（n-1个0）以及特征函数（n-1个是原来的正交基，有一个变换成了y）。

 

二、一般矩阵$A \in \mathbb{C}^{n\times n}$ 与幂乘法

非亏矩阵$A$可以表示为

$$A=X\Lambda Y^T$$

其中$XY^T = I, \Lambda = diag(\lambda_1, \lambda_2, \cdots, \lambda_n)$。用$x_i$，$y_i$表示$X$，$Y$的第$i$列，则$A$可以写成

$$ A=\sum\limits_{i=1}^{n}\lambda_ix_iy_i^T $$

即$n$个秩1矩阵之和。

不妨设$|\lambda_1|>|\lambda_2|\geqslant \cdots\geqslant |\lambda_n|$，并记$\widetilde{A} =x_1y_1^T$是第1个秩1矩阵，则有

$$\left |\left|\dfrac{A^k}{\lambda_1^k}-\widetilde{A}\right|\right|_2 =\left |\left|\sum\limits_{i=2}^{n}\left (\dfrac{\lambda_i}{\lambda_1}\right )^kx_iy_i^T\right|\right|_2 \leqslant \eta\left|\dfrac{\lambda_2}{\lambda_1}\right|^k \rightarrow 0,\ k\rightarrow \infty$$

其中$\eta=(n-1)\max\limits_i||x_iy_i^T||_2.$

以上结果表明，当k充分大的时候，$\dfrac{A^k}{\lambda_1^k}$与秩1矩阵$\widetilde{A}$非常接近，由上文可知，若取$\mathbb{C}^n$中某一非零向量$x$满足$\dfrac{A^k}{\lambda_1^k}x \neq 0$，则得到矩阵$A$的一个很好的近似特征向量$u^{(k)}=\dfrac{A^k}{\lambda_1^k}x$

如果作为一种求近似特征向量的方法，以上思路至少有两个问题，但是都可以解决：

问题1. $\lambda_1$并不知道；解决：$\lambda_1$实际上只对向量的模长起作用，可以用归一化替代；

问题2. $A^k$计算量太大；解决：迭代计算$A^kx$。

至此，可以设计迭代求最大特征值对应特征向量的迭代格式：

$ y_k=Au_{k-1}, $

$ \mu_k=\zeta_j^{(k)}, $  $\zeta_j^{(k)}$是$y_k$的模最大分量 //

$ u_k=y_k/\mu_k,\qquad k=1,2,\cdots $

其中$u_0 \in \mathbb{C}^n$是任意给定的初始向量，$||u_0||_\infty =1$

只要$ u_0 $选择较为合理，则迭代产生的$\{u_k\}$收敛于对应于最大特征值$\lambda_1$的特征向量$v_1$，$\{\mu_k\}$收敛于$\lambda_1$，收敛速度依赖于$\left|\dfrac{\lambda_2}{\lambda_1}\right|$的大小。

注意，幂乘法能够求出矩阵最大的特征值及其对应得特征向量。如果对应$\textbf{HC}=E\textbf{C}$特征方程求最小特征值，这样的思路是可以借鉴的。这一点在后面的文章中会涉及到。

 

三、幂乘法的几何意义与正交迭代法

上面的幂乘法相当于是将一个由任意非零向量$u_0$张成的一维子空间$S=span\{u_0\}$迭代产生一维子空间序列

$$S, AS,\cdots, A^kS,\cdots$$

随着$k$增大$A^kS$将会逼近$A$的特征子空间$span\{v_1\}$。

将这一思路推广，如果对$l$维的任意子空间$S$做迭代形成序列，最终能够收敛到$l$维的特征子空间。（证明略）

迭代格式如下：

$ Z_k=AQ_{k-1}, $

$ Q_kR_k=Z_k, \qquad k=1,2,\cdots$ 

其中$Q_0$是任意的l维子空间$S$的一组正交基，是$n \times l$维矩阵。

$Q_k^*Q_k=I_l$，$Q_k$是$n \times l$维矩阵（*表示取厄米共轭，因为$Q_k$是复矩阵），$R$是$l\times l$维的上三角矩阵。最终迭代将得到$l$维特征子空间的一组正交基。这种迭代法即正交迭代法。

 

四、QR迭代法

由上节自然想到，取$l = n$ 为矩阵$A$的维度时，将会得到全部的特征向量和特征值。定义

$$ T_k=Q_k^*AQ_k $$

则$T_k=Q_k^*Z_{k+1}=Q_k^*Q_{k+1}R_{k+1} \rightarrow R_k,\ k\rightarrow \infty$，得到的是矩阵$A$的Schur分解。所以改进迭代方案：

$ T_0=A, $

$ Q_kR_k=T_{k-1},$ （QR分解）

$ T_k=R_kQ_k, \qquad k=1,2,\cdots$ 

此即复空间的QR迭代法。

 

五、实空间QR迭代法

实空间得到Schur分解虽然同样可行，但是实际上效率较低。通常采用上Hessenberg化，并且上Hessenberg矩阵的QR分解可以用Givens变换实现（在此不展开）。

具体的实现算法有双重步位移QR算法。