1、QR迭代法求解矩阵A的特征值沈欢00986096北京大学工学院，北京1008712011 年 10 月 23 日摘摘摘 要要要本文用QR迭代法求解矩阵A的特征值:第一步先用豪斯荷尔德变换将矩阵A化为上海森伯格矩阵AH，第二步再对AH进行QR迭代(使用吉文斯变换)，当迭代满足精度要求时输出n个特征值。1问问问题题题描描描述述述对于给定的矩阵A(存在n个不相等的实特征值),采用QR方法求其n个特征值。本文中的矩阵A由文件”gr900900crg.mm”给出。并以.mm格式描述。2算算算法法法描描描述述述2.1得得得到到到实实实对对对称称称矩矩矩阵阵阵A本文由文件”gr900900crg.mm”得到
2、了以.mm格式描述的系数矩阵A。A矩阵是900900的大型稀疏对称矩阵。于是，在matlaB中，使用”A=zeros(900,900)”语句生成900900的零矩阵。再按照.mm文件中的描述，分别对第i行、第j列的元素赋对应的值，就生成了系数矩阵A，并将A存为.mat文件以便之后应用。由于A矩阵是实对称的900阶矩阵，所以A矩阵有900个实特征值，可以使用QR方法求解。下面描述具体算法。2.2QR算算算法法法QR算法的核心思想是：对一个具有n个不相等的特征值的矩阵A，进行如下变换：A=Q1R1;A1=R1Q1A1=Q2R2;A2=R2Q21Ak=Qk+1Rk+1;Ak+1=Rk+1Qk+1其中
3、Qi是正交矩阵，Ri是可逆的上三角矩阵。由于：Ai+1=QTi+1(Qi+1Ri+1)Qi+1=QTi+1AiQi+1所以Ai+1和Ai相似，具有相同的特征根。当i 时，Ai对角线以下元素趋于零。利用这一原理，可以对给定A进行QR迭代，最终得到上三角矩阵(对角元素不一定非零)，输出对角元素可以作为A的特征值。我们又发现，具有如下上海森伯格形式的矩阵AH：AH=a11a12a13.a1n1a1na21a22a23.a2n1a2n0a32a33a34.00.an,n1ann(1)AH在QR迭代过程中可以保持上海森伯格形式不变。所以，QR迭代算法分两步实现，第一步将矩阵A变换为上海森伯格矩阵AH,第
4、二步对此上海森伯格矩阵AH进行QR迭代，直至满足精度为止，输出结果。在第一步中，将矩阵A变换为上海森伯格矩阵AH，为了保持变换后的AH与A具有相同的特征值，这里采用豪斯荷尔德变换(正交变换)，本步在后文中详细说明；在第二步中，对上海森伯格矩阵AH进行QR迭代，并不去直接找出Q矩阵和R矩阵，而是采用吉文斯变换达到相同的效果，本步在后文中也有详细说明；对于QR迭代的精度，在算法采用：当迭代到第k次时有：|Ai+1,i|Ai,i|+|Ai+1,i+1|accuracy对于i=1,2,n-1成立，其中accuracy是所设精度。认为此时的AH次对角元素已经为零，迭代结束，输出该矩阵对角元素作为A的特征
5、值。QR算法主程序(见程序：function eigenvalue,flag1=QR(A,n)的程序流程图如图1所示。对于方框中的两步：1、用豪斯荷尔德变换将A矩阵化为上海森伯格矩阵AH;2、用吉文斯变换实现一次等效的QR迭代；将分别描述。2.3用用用豪豪豪斯斯斯荷荷荷尔尔尔德德德变变变换换换将将将A矩矩矩阵阵阵化化化为为为上上上海海海森森森伯伯伯格格格矩矩矩阵阵阵AH豪斯荷尔德矩阵H=I 2wwT其中w是单位向量。所以H矩阵是正交、对称矩阵。2?A?AH?QR?Ai=QiRiAi+1=RiQi?QR?An?An?A?A?NoYes?QR?A?图 1:主程序(QR迭代法求矩阵特征值)算法用豪斯
6、荷尔德变换将A矩阵化为上海森伯格矩阵AH是实现下述过程：A1=H1AH1A2=H2A1H2An2=Hn2An3Hn2其中H1、H2 满足:H1a11a21a31.an1=a011a0210.0(2)3H2a012a022a032a042.a0n2=a0012a0022a00320.0(3)于是构造豪斯荷尔德矩阵Hi成为此函数的核心。对于Hi，目的是将之前得到的Ai1中的第i列实现如下变换：Hia01ia02ia03ia0iiai+1,iai+2,i.a0ni=a001ia002ia003ia00iia00i+1,i0.0(4)于是从Ai1中取出n i维列向量 x=a0i+1,i,a0i+2,i
7、,a0niT取=k x k2;和n-i维单位向量 e=1,0,0,0T构造=x ek x e k2于是就得到了(n i)(n i)的豪斯荷尔德矩阵H0i,再将其扩充：Hi=Iii00H0i#(5)得到Hi。实现:Ai=HiAi1Hi将此步依次从i=1至i=n 2操作n-2次，得到上海森伯格矩阵AH,且与原矩阵A相似。本函数(function AH=A-to-AHessenberg(A,n)流程图如图二所示。4i=n-2?n-i?wi?wi?(n-i)*(n-i)?H(n-i)*(n-i)?H(n-i)*(n-i)?HiAi=Hi*Ai-1*HiA?Ai-1Ai?AH?A?AH图 2:用豪斯荷尔
8、德变换将矩阵A化为上海森伯格矩阵AH算法2.4用用用吉吉吉文文文斯斯斯变变变换换换实实实现现现一一一次次次等等等效效效的的的QR迭迭迭代代代在上文，已用豪斯荷尔德变换将矩阵A化为上海森伯格矩阵AH.对于一次QR迭代而言，相当于用吉文斯变换实现等效地实现以下过程：G(n 1,n;n1)G(n 2,n 1;n1)G(1,2;1)AH=R其中G(i 1,i;i1)是吉文斯扩充矩阵，是正交矩阵。该过程等效于：AH=QR其中Q=G(1,2;1)TG(n 2,n 1;n1)TG(n 1,n;n1)T得到上三角矩阵R后，再进行：AHnew=RQ=RG(1,2;1)TG(n 2,n 1;n1)TG(n 1,n
9、;n1)T5但是如果真的将Q矩阵找出，并如上述作用，那么计算开销会非常大。所以我们考虑一个等效的过程。对于第一个过程，实际上是依次将AH矩阵的第i行和第i+1行进行变换(总共进行n-1次变换)，目的是使得AHi+1,i=0。于是，在第一个过程中就只进行第i行和第i+1行的变换。取：ci=AHi,iqAH2i,i+AH2i+1,isi=AHi+1,iqAH2i,i+AH2i+1,i于是，第i行和第i+1行的元素变换为：AH0i,j=ci AHi,j+si AHi+1,jAH0i+1,j=ci AHi+1,j si AHi,jj=1,2,n.对i从1至n-1,重复以上过程,则等效实现QR分解，并得
10、到上三角矩阵R。在第二个过程中，R矩阵右乘上Q矩阵，即：AHnew=RQ=RG(1,2;1)TG(n 2,n 1;n1)TG(n 1,n;n1)T等效于依次对R矩阵的第i列和第i+1列进行列变换。变换后的元素为：R0j,i=ci Rj,i+si Rj,i+1R0j,i+1=ci Rj,i+1 si Rj,ij=1,2,n.对i从1至n-1,重复以上过程,则等效地实现AHnew=RQ过程。至此，完成了一次QR迭代。本函数(AHnew=QR-iteration-once(AH,n)的流程图如图三所示。3结结结果果果输输输出出出在functioneigenvalue,flag 1=QR(A,n)函数
11、中，设置accuracy=103,取n=900,A为本文2.1所得的矩阵，运行该函数。结果显示，flag1=9344。也就是，经过9344次迭代，上海森伯格矩阵次AH对角元素接近0，可由对角元素得到A的特征值。所得结果在同一文件夹中的eigenvalue.mat中图四显示了900个特征值从小到大的排列。表一给出了900个特征值中最大的八个和最小的八个特征值。6?G(i,i+1)?AH?i?i+1?AH i+1,i=0?A=QR.?G(1,2)T,G(2,3)T?G(n,n-1)T?R?i?i+1?Anew=RQ?AHi=1:n-1?RNOYESi=1:n-1NOYES?QR?Anew?AH?QR?图 3:用吉文斯变换等效实现一次QR迭代算法表 1:八个最大的和最小的特征值最大特征值11.959111.959111.928711.928711.878411.878311.867311.8668最小特征值0.06150.15320.15320.24400.30500.30500.39420.394270100200300400500600700800900024681012图 4:A的900个特征值从小到大排列8
