普通矩阵特征值的QR算法

普通矩阵特征值的QR算法摘要求矩阵的特征值有多种不同的办法，本文主要介绍用QR法求矩阵的特征值，QR法是目前求中等大小矩阵全部特征值的最有效方法之一，使用于求实矩阵或复矩阵的特征值，它和雅可比法类似，也是一种变换迭代法。关键词：QR分解迭代序列特征值Matlab一、QR方法的理论：对任意一个非奇异矩阵（可逆矩阵）A，可以把它分解成一个正交阵Q和一个上三角阵的乘积，称为对矩阵A的QR分解，即A=QR。如果规定R的对角元取正实数，这种分解是唯一的。若A是奇异的，则A有零特征值。任取一个不等于A的特征值的实数μ，则A-μI是非奇异的。只要求出A-μI的特征值和特征向量就容易求出矩阵A的特征值和特征向量，所以假设A是非奇异的，不是一般性。设A=A1，对A1作QR分解，得A1=Q1R1，，交换该乘积的次序，得A2=R1Q1=，由于Q1正交矩阵，A1到A2的变换为正交相似变换，于是A1和A2就有相同的特征值。一般的令A1=A，对k=1，2，3，…..)()(1迭代定义分解kkkkkkQRAQRRQA这样，可得到一个迭代序列{Ak}，这就是QR方法的基本过程。二、QR方法的实际计算步骤HouseholderAHessenbergB用阵作正交相似变换上第阵一步............*::::*Householder变换：如果v给出为单位向量而I是单位矩阵，则描述上述线性变换的是豪斯霍尔德矩阵（v*表示向量v的共轭转置）H=I-2VV*1kkkGivenkkkBQRBBRQ用变换产生迭代序列第二步12***nHouseholderAB用阵作正交相似变换（对称阵）三对角阵***三、化一般矩阵为Hessenberg阵称形如11121112122212323331nnnnnnnnnhhhhhhhhhhhHhh的矩阵为上海森堡（Hessenberg）阵。如果此对角线元(i=2,3,….n)全不为零，则该矩阵为不可约的上Hessenberg矩阵。用Householder变换将一般矩阵A相似变换为Hessenberg矩阵。首先，选取Householder矩阵，使得经相似变换后的矩阵的第一列中有尽可能多的零元素。为此，应取为如下形式1110000HH其中1H为n-1阶Householder矩阵。于是有111211111221TaaHHAHHaHAH121311212131(,,,),(,,,),TTnnaaaaaaaa222222nnnnaaAaa只要取1H使得111(,0,,0)THa，就会使得变换后的矩阵11HAH的第一列出现n-2个零元。同理，可构造如下列形式Householder矩阵22100001000000HH使得2112**************HHAHH如此进行n-2次，可以构造n-2个Householder矩阵12,,,HH2nH，使得221122nnHHHAHHHH。其中H为上Hessenberg矩阵。特别地，当A为实对称矩阵，则经过上述正交变换后，H变为三对角阵。用Householder方法对矩阵A作正交相似变换，使A相似于上Hessenberg阵，算法如下：1,2,...,2kn（1）计算(1)(1)111221111111(1)112,1()()(()),()(0,...,0,,,...,)kkTkknkkkikikkkkkkkkkkkknkHIUUsignaaaUaaa，，，（2）计算1kHAA(1)11(1),1,,11()21,,nkjlljlkkkijijjijkkntuaiknaatu（）（）（3）计算1kAHA(1)11(1)1,...,1(1)(2)1,...,nkiilllkkkijijijintaujknaatu四、上Hessenberg阵的QR分解对上Hessenberg阵只需要将其次对角线上的元素约化为零，用Given变换比用Householder变换更节省计算量。1、平面旋转阵（Givens变换阵）定义n阶方阵,11cossin11sincos11()ijiRjji称为平面旋转阵，或称为Givens变换阵。平面旋转阵Ri,j的性质：（1）1,,,,,TTijijijijRRIRR平面旋转阵是非对称的正交阵。（2）,TijR也是平面旋转阵。（3）,ijRdet()=1平面旋转阵Ri,j的作用：（1）将向量,,...,xxxxT12n=的第j个分量约化为零。,ijR左乘向量,,...,xxxxT12n=只改变X的第i个分量和第j个分量。若令,ijyRx，有cossinsincos1,...,;,iijjijkkyxxyxxyxknkij,调整，可将jy约化为零。令0jy，得tanjixx所以，取22cosiiijxxCrxx，22sinjjijxxSrxx于是22,0iijijjyCxSxrxxy,,...,,,,...,,0,,...,ijRxxxrxxxxT1i-1i+1j-1j+1n=（2）将向量,,...,xxxxT12n=的第i+1个分量到第n个分量约化为零。22,11,...,,,0,,...,,iiiiRxxxrxxrxxT1i-1i+2n=,2,122212,...,,,0,0,,...,,iiiiiiiRRxxxrxxrxxxT1i-1i+3n=,,2,122,...,,,0,...,0,iniiiiinRRRxxxrrxxT1i-1=（3）用,ijR对矩阵A作变换得到的结论,ijR左乘A只改变A的第i，j行。,ijRT右乘A只改变A的第i，j列。,,ijijRART只改变A的第i，j行和第i，j列。2、用Givens变换对上Hessenberg阵作QR分解对上Hessenberg阵(1)(1)(1)11121(1)(1)(1)21222(1)(1)1nnnnnnbbbbbbBbb通常用n-1个Givens变换阵可将它化成上三角矩阵，从而得到B的QR分解式。具体步骤为：设(1)210b（否则进行下一步）取旋转矩阵1111cossin00sincos00(1,2)0011R则(2)(2)(2)112131(2)(2)(2)22232(2)(2)(2)1232333(2)(2)1(1,2)nnnnnnnrbbbbbbRBBbbbbb其中，(1)(1)(1)(1)1121111112111cos,sin,bbrbbrr设2320（）b（否则进行下一步），再取旋转矩阵222210cossinsincos(3,2)11－R则(3)(3)(3)(3)11213111(3)(3)(3)223212(3)(3)(3)3331323(3)(3)(3)43414(3)(3)1(3,2)nnnnnnnnnnnnrbbbbrbbbbbbRBBbbbbb其中(2)(2)(2)2(2)23222222223222cos,sin,()()bbrbbrr。假设上述过程已进行了k-1步，有1(1,)kkBRkkB()()()()1111111()()()11111()()()1()()()1111()()1kkkkkknnkkkkkkknknkkkkkknknkkkkkknknkknnnnrbbbbrbbbbhhbbbbb设()10kkkb，取11(1,)cossinsincos1kkkkRkk其中()()1cos,sinkkkkkkkkkkbbrr，()2()21()().kkkkkkkrbb于是(1)(1)(1)11111(1)(1)1(1)(1)(1)111111(1)(1)(1)21212(1)(1)1(1,)kkkkknkkkkkknkkkkkkkknknkkkkkknknkknnnnrbbbrbbRkkBBbhbbhhbb因此，最多做n-1次旋转变换，即得()(,1)(2,1)(1,2)nHRnnRnnRB()()()112131()()2232()33nnnnnnnnnnrbbbrbbRrbr因为(,1),(2,3,,)Riiin均为正交矩阵，故21321TTTnnHRRRRQR其中21321TTTnnQRRR仍为正交矩阵。可算出完成这一过程的运算量约为4,比一般矩阵的QR分解的运算量3()On少一个数量级。可证明仍是上Hessenberg阵，于是可按上述步骤一直迭代下去，这样的得到的QR方法的运算量比基本QR方法大为减少。具体实例现在我们就用上述所说的QR方法求解矩阵532644445A的全部特征值。第一步：将A化成上Hessenberg阵，取22(6,4)64(1,0)(652,4)TTTu10201TTuuIuu0.9160250.2773500.8320500.55470020.2773500.08397470.5547000.832050于是110000.8320500.55470000.5547000.832050H1151.3867503.3282007.2111021.2307688.15384000.1538462.230767HHAHH即为与A相似的上Hessenberg阵。将H进行QR分解，记2211,5(7.21102)8.774964BHr111cos50.56980.sin0.821781r取0.5698030.8217810(2,1)0.8217810.5698030001R于是18.7749641.8015968.597089(2,1)00.4383101.91103000.1538462.230767RB再取于是1(3,2)(2,1)RRB18.7749641.8015968.59708900.4645262.541982001.471953R10.5698030.7754030.272165(2,1)(3,2)0.8217810.5376430.18871200.3311890.943564TTQRR第一次迭代得2113.5194824.92549110.8401170.3817391.0916272.31065300.4874951.388883BRQ重复上述过程，迭代11次得122.9920321.000385312.0133920.0074962.0046951