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

更新: 28 JUL 2016

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

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

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

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

$$ A=uv^T $$

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

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

$$ y = Ax $$

若$y = 0$则说明$x$是矩阵$A$对应于特征值0的特征向量。若$y eq 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 imes n}$ 与幂乘法

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

$$A=XLambda Y^T$$

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

$$ A=sumlimits_{i=1}^{n}lambda_ix_iy_i^T $$

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

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

$$left |left|dfrac{A^k}{lambda_1^k}-widetilde{A} ight| ight|_2 =left |left|sumlimits_{i=2}^{n}left (dfrac{lambda_i}{lambda_1} ight )^kx_iy_i^T ight| ight|_2 leqslant etaleft|dfrac{lambda_2}{lambda_1} ight|^k ightarrow 0, k ightarrow infty$$

其中$eta=(n-1)maxlimits_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 eq 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} ight|$的大小。

注意,幂乘法能够求出矩阵最大的特征值及其对应得特征向量。如果对应$ extbf{HC}=E extbf{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 imes l$维矩阵。

$Q_k^*Q_k=I_l$,$Q_k$是$n imes l$维矩阵(*表示取厄米共轭,因为$Q_k$是复矩阵),$R$是$l imes 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} ightarrow R_k, k ightarrow 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算法。

原文地址:https://www.cnblogs.com/fnight/p/5716937.html