普通矩阵特征值的QR算法.pdf

普通矩阵特征值的 QR 算法普通矩阵特征值的 QR 算法摘摘要要求矩阵的特征值有多种不同的办法，本文主要介绍用 QR 法求矩阵的特征值，QR 法是目前求中等大小矩阵全部特征值的最有效方法之一，使用于求实矩阵或复矩阵的特征值，它和雅可比法类似，也是一种变换迭代法。关键词：关键词：QRQR 分解分解迭代序列迭代序列特征值特征值MatlabMatlab一一 、QRQR 方法的理论：方法的理论：对任意一个非奇异矩阵（可逆矩阵）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，. Ak QkRkAk1 RkQk(QR分解)(迭代定义)这样，可得到一个迭代序列Ak，这就是 QR 方法的基本过程。二、二、QRQR 方法的实际计算步骤方法的实际计算步骤 . . . . * *第第一步一步 A A 上上HessenbergHessenberg阵阵B B 用用HouseholderHouseholder阵阵 作正交相似变换作正交相似变换 1 / 10. . : : : : : : : : * * 普通矩阵特征值的 QR 算法HouseholderHouseholder 变换：变换：如果 v 给出为单位向量而 I 是单位矩阵，则描述上述线性变换的是 豪斯霍尔德矩阵豪斯霍尔德矩阵 （v*表示向量 v 的共轭转置）H=I -2VV* 1 1* *第二步第二步 2 2 B Bk k Q Qk kR Rk kB B 用用GivenGiven变换产生迭代序列变换产生迭代序列B B R R Q Qk kk k k k 1 1 * * A A （对称阵）（对称阵） 三对角阵三对角阵B B 用用HouseholderHouseholder阵阵 作正交相似变换作正交相似变换 * * * * * * n n * *三、化一般矩阵为三、化一般矩阵为 HessenbergHessenberg 阵阵称形如h h1 1n n 1 1h h1 1n n h h1111h h1212 h h h hh hh h22222 2n n 1 12 2n n 2121h h3232h h3333h h3 3n n H H h hnnnn 1 1h hnnnn 的矩阵为上海森堡（Hessenberg）阵。 如果此对角线元(i=2,3,.n) 全 不 为零，则该矩阵为不可约的上 Hessenberg 矩阵。用 Householder 变换将一般矩阵 A 相似变换为Hessenberg 矩阵。首先，选取 Householder 矩阵，使得经相似变换后的矩阵的第一列中有尽可能多的零元素。为此，应取为如下形式 1 10 0 0 0H H1 1 H H1 1 0 0 0 0 其中H H1 1为 n-1 阶 Householder 矩阵。 a a1111于是有H H1 1AHAH1 1 H H1 1a a1 1T Ta a2 2H H1 1 H H1 1A A2222H H1 1 2 / 10普通矩阵特征值的 QR 算法a a1 1 ( (a a2121, ,a a3131, , ,a an n1 1) )T T, ,a a2 2 ( (a a1212, ,a a1313, , a a2222 , ,a a1 1n n) )T T, ,A A2222 a an n2 2a a2 2n n a annnn 只要取H H1 1使得H H1 1a a1 1 ( ( 1 1,0,0,0),0)T T，就会使得变换后的矩阵H H1 1AHAH1 1的第一列出现 n-2 个零元。同理，可构造如下列形式 Householder 矩阵0 0 1 10 00 0 * * * * 0 01 10 0 * * * *0 0 使得H H2 2H H1 1AHAH1 1H H2 2 * * *H H2 2 0 00 0 H H 2 2 0 00 0 * * * * * * * * * * 如此进行 n-2 次，可以构造 n-2 个 Householder 矩阵H H1 1, ,H H2 2, , ,H Hn n 2 2，使得H Hn n 2 2H H2 2H H1 1AHAH1 1H H2 2H Hn n 2 2 H H。其中 H 为上 Hessenberg 矩阵。特别地，当 A为实对称矩阵，则经过上述正交变换后，H 变为三对角阵。用 Householder 方法对矩阵 A 作正交相似变换， 使 A 相似于上 Hessenberg阵，算法如下：k k 1, 1,2,.,2,., n n 2 2（1）计算H Hk k 1 1 I I 1 1 k k 1 1U U( (k k 1)1)( (U U( (k k 1)1) )T Tn n1 12 22 2 k k 1 1 signsign( (a ak k 1 1，k k)( )( ( (a aik ik) ) ) ) , ,i i k k 1 1 k k 1 1 k k 1 1( ( k k 1 1 a ak k 1 1，k k) )U U( (k k 1)1) (0,.,0,(0,.,0, k k 1 1 a ak k 1 1，k k, ,a ak k 2, 2,k k,.,., a anknk) )（2）计算H Hk k 1 1A A A Aj j k k, ,k k 1, 1, ,n nn n1 1（1 1）t tj j ( ( u ul l( (k k 1)1)a alj lj) )（2 2）i i k k 1, 1, ,n na aij ij a aij ij t tj ju ui i( (k k 1)1) k k 1 1l l k k 1 1（3）计算AHAHk k 1 1 A A3 / 10普通矩阵特征值的 QR 算法i i 1,.,1,.,n n1 1(1)(1)t ti i k k 1 1l l k k 1 1 a a u uil iln n( (k k 1)1)l l(2)(2) j j k k 1,.,1,.,n na aij ij a aij ij t ti iu u( (j jk k 1)1)四、上四、上 HessenbergHessenberg 阵的阵的 QRQR 分解分解对上 Hessenberg 阵只需要将其次对角线上的元素约化为零， 用 Given 变换比用 Householder 变换更节省计算量。1 1、 平面旋转阵（平面旋转阵（GivensGivens 变换阵）变换阵）定义n 阶方阵1 11 1coscossinsin1 1R Ri i, , j j1 1sinsincoscos1 1称为平面旋转阵，或称为 Givens 变换阵。i i j j1 1( ( j j i i) )平面旋转阵R Ri,ji,j的性质：T T（1）R Ri iT T, ,j jR Ri i, ,j j I I, , R Ri i, ,1 1j j R Ri i, ,j j平面旋转阵是非对称的正交阵。（2）R Ri iT T, , j j也是平面旋转阵。（3）det(det(R Ri i, , j j)=1)=1平面旋转阵R Ri,ji,j的作用：（1） 将向量x x= =x x1 1, , x x2 2,.,., x xn n的第 j 个分量约化为零。T TR Ri i, , j j左乘向量x x= =x x1 1, , x x2 2,.,., x xn n只改变 X 的第 i 个分量和第 j 个分量。4 / 10T T普通矩阵特征值的 QR 算法若令y y R Ri i, ,j jx x，有y yi i x xi icoscos x xj jsinsiny yj j x xi isinsin x xj jcoscosy y x xk k 1,.,1,.,n n; ;k k i i, , j jk k , ,k k调整，可将y yj j约化为零。令y yj j 0 0，得tantan所以，取C C coscos x xj jx xi ix xi ix x i i，S S sinsin r rx xi i2 2 x x2 2j jx xj jx x x x2 2i i2 2j j x xj jr r于是y yi i CxCxi i SxSxj j r r x xi i2 2 x x2 2j j, ,y yj j 0 0T TR Ri i, , j jx x= = x x1 1,.,., x xi-1i-1, ,r r, , x xi+1i+1,.,., x xj-1j-1,0,0, x xj+1j+1,.,., x xn n （2） 将向量x x= =x x1 1, , x x2 2,.,., x xn n的第i+1i+1个分量到第 n 个分量约化为零。R Ri i, ,i i1 1x x= =x x1 1,.,., x xi-1i-1, ,r r,0,0, x xi+2i+2,.,., x xn n, , r r x xi i2 2 x xi i2 21 1T TT TR Ri i, ,i i2 2R Ri i, ,i i1 1x x= =x x1 1,.,., x xi-1i-1, ,r r,0,0,0,0, x xi+3i+3,.,., x xn n, ,r r x x x xR Ri i, ,n n2 2i i2 2i i1 1T T x x2 2i i2 2T TR Ri i, ,i i2 2R Ri i, ,i i1 1x x= =x x1 1,.,., x xi-1i-1, ,r r,0,.,0,0,.,0, ,r r x x 2 2i i x x2 2n n（3） 用R Ri i, , j j对矩阵 A 作变换得到的结论R Ri i, , j j左乘 A 只改变 A 的第 i，j 行。R Ri iT T, , j j右乘 A 只改变 A 的第 i，j 列。R Ri i, , j jARARi iT T, ,j j只改变 A 的第 i，j 行和第 i，j 列。2 2、用、用 GivensGivens 变换对上变换对上 HessenbergHessenberg 阵作阵作 QRQR 分解分解(1)(1) b b1111 (1)(1)b b2121 对上 Hessenberg 阵B B (1)(1)b b1212(1)(1)b b2222(1)(1)b bnnnn 1 1 b b1 1(1)(1)n n(1)(1) b b2 2n n (1)(1) b bnnnn 5 / 10普通矩阵特征值的 QR 算法通常用 n-1 个 Givens 变换阵可将它化成上三角矩阵，从而得到 B 的 QR 分解式。具体步骤为：(1)(1)设b b2121 0 0（否则进行下一步） coscos 1 1 sinsin 1 1 取旋转矩阵R R(1,(1,2)2) 0 0 sinsin 1 1coscos 1 10 00 00 01 10 0 0 0 1 1 r r1 1 则R R(1,(1,2)2)B B1 1 (2)(2)b b1212(2)(2)b b2222(2)(2)b b3232(2)(2)b b1313(2)(2)b b2323(2)(2)b b3333(2)(2)b bnnnn 1 1 b b1 1(2)(2)n n(2)(2) b b2 2n n (2)(2) B B2 2b b3 3n n (2)(2) b bnnnn (1)(1)(1)(1)b b1111b b2121(1)(1)(1)(1)其中，coscos 1 1 , ,sinsin 1 1 , ,r r1 1 b b1111 b b2121r r1 1r r1 1（2 2）设b b（否则进行下一步） ，再取旋转矩阵3232 0 0 1 1 coscos 2 2 sinsin 2 2R R(3,(3,2)2) sinsin 2 2coscos 2 21 10 0 1 1 r r1 1 则R R(3,(3,2)2)B B2 2 (3)(3)b b1212r r2 2(3)(3)b b1313(3)(3)b b2323(3)(3)b b3333(3)(3)b b4343b b1 1(3)(3)n n 1 1(3)(3)b b2 2n n 1 1(3)(3)b b3 3n n 1 1(3)(3)b b4 4n n 1 1(3)(3)b bnnnn 1 1 b b1 1(3)(3)n n(3)(3) b b2 2n n (3)(3) b b3 3n n B B3 3(3)(3) b b4 4n n (3)(3)b bnnnn (2)(2)(2)(2)b b3232b b2222(2)(2)2 2(2)(2)2 2其中coscos 2 2 , ,sinsin 2 2 , ,r r2 2 ( (b b2222) ) ( (b b3232) )。r r2 2r r2 2假设上述过程已进行了 k-1 步，有6 / 10普通矩阵特征值的 QR 算法 r r1 1 B Bk k R R( (k k 1, 1,k k) )B Bk k 1 1 ( (k k) )设b bk k，取 1 1k k 0 0k k) )b b1 1( (k k 1 1k k) )b b1 1( (k kk k) )b b1 1( (n n 1 1r rk k 1 1( (k k) )b bk k 1 1k k( (k k) )b bkkkk( (k k) )b bk k 1 1k k( (k k) )b bk k 1 1n n 1 1( (k k) )h hknkn 1 1( (k k) )b bk k 1 1n n 1 1( (k k) )b bnnnn 1 1k k) ) b b1 1( (n n ( (k k) ) b bk k 1 1n n( (k k) ) h hknkn ( (k k) ) b bk k 1 1n n ( (k k) ) b bnnnn 1 1 R R( (k k 1, 1,k k) ) 1 1coscos k k sinsin k ksinsin k kcoscos k k1 1 ( (k k) )( (k k) )b bkkkkb bk k( (k k) )2 2( (k k) )2 2) ) ( (b bk k其中coscos k k , ,sinsin k k 1 1k k，r rk k ( (b bkkkk 1 1k k) ) . .r rk kr rk k r r1 1 于是R R( (k k 1, 1,k k) )B Bk k k k 1)1)b b1 1( (k kk k 1)1)b b1 1( (k k 1 1r rk k( (k k 1)1)b bkkkk 1 1( (k k 1)1)b bk k 1 1k k 1 1( (k k 1)1)b bk k 2 2k k 1 1( (k k 1)1)h hk k 1 1n n 1 1( (k k 1)1)h hk k 2 2n n 1 1( (k k 1)1)b bnnnn 1 1k k 1)1) b b1 1( (n n ( (k k 1)1) b bknkn( (k k 1)1) B Bk k 1 1b bk k 1 1n n ( (k k 1)1) h hk k 2 2n n ( (k k 1)1) b bnnnn 因此，最多做 n-1 次旋转变换，即得 r r1 1 R R(1,2)(1,2)B B ( (n n) )b b1212( (n n) )b b1313( (n n) )b b2323r r2 2H H( (n n) ) R R( (n n, ,n n 1)1)R R( (n n 2, 2,n n 1)1)r r3 3n n) ) b b1 1( (n n( (n n) ) b b2 2n n ( (n n) ) R Rb b3 3n n r rn n T TR Rnnnn 1 1R R QRQR其因为R R( (i i, ,i i 1),(1),(i i 2,3,2,3,T TT T中Q Q R R2121R R3232T TT T, ,n n) )均为正交矩阵，故H H R R2121R R3232T TR Rnnnn 1 1仍为正交矩阵。可算出完成这一过程的运算量约为4,比一般矩阵的 QR 分解的运算量O O( (n n3 3) )少一个数量级。7 / 10普通矩阵特征值的 QR 算法可证明仍是上 Hessenberg 阵，于是可按上述步骤一直迭代下去，这样的得到的 QR 方法的运算量比基本 QR 方法大为减少。具体实例具体实例 5 5 3 32 2 现在我们就用上述所说的 QR 方法求解矩阵A A 6 6 4 44 4 的全部特征值。 4 4 4 45 5 第一步：将 A 化成上 Hessenberg 阵，取u u (6,4)(6,4)T T 6 62 2 4 42 2(1,0)(1,0)T T (6(6 52,4)52,4)T TuuuuT T 1 10 0 I I 2 2T T u u u u 0 01 1 0.9160250.9160250.2773500.277350 0.8320500.832050 0.5547000.554700 2 2 0.2773500.2773500.08397470.0839747 0.5547000.5547000.8320500.832050 0 00 0 1 1 于是H H1 1 0 0 0.8320500.832050 0.5547000.554700 0 0 0.5547000.5547000.8320500.832050 5 51.3867501.3867503.3282003.328200 H H H H1 1AHAH1 1 7.2111027.211102 1.2307681.230768 8.1538408.153840 0 0 0.1538460.1538462.2307672.230767 H 即为与 A 相似的上 Hessenberg 阵。将 H 进行 QR 分解，记B B1 1 H H, ,r r1 1 5 52 2 ( ( 7.21102)7.21102)2 2 8.7749648.774964coscos 1 1 5 5 r r1 1 0.56980.0.56980.sinsin 1 1 0.8217810.821781 0.5698030.569803 0.8217810.8217810 0 取R R(2,1)(2,1) 0.8217810.8217810.5698030.5698030 0 0 00 01 1 8 / 10普通矩阵特征值的 QR 算法 8.7749648.7749641.8015961.8015968.5970898.597089 于是R R(2,1)(2,1)B B1 1 0 00.4383100.438310 1.9110301.911030 0 0 0.1538460.1538462.2307672.230767 再取r r2 2 (0.438310)(0.438310)2 2 ( ( 0.153846)0.153846)2 2 0.464526,0.464526,coscos 2 2 0.4383100.438310 r r2 2 0.943564,0.943564,sinsin 2 2 0.1538460.153846 r r2 2 0.3311890.3311890 00 0 1 1 R R(3,(3,2)2) 0 00.9435640.943564 0.3311890.331189 0 00.3311890.3311890.9435640.943564 于是R R(3,2)(3,2) R R(2,1)(2,1)B B1 1 8.7749648.7749641.8015961.8015968.5970898.597089 R R 0 00.4645260.464526 2.5419822.5419821 1 0 00 01.4719531.471953 0.5698030.5698030.7754030.7754030.2721650.272165 Q Q1 1 R RT T(2,1)(2,1)R RT T(3,2)(3,2) 0.8217810.8217810.5376430.5376430.1887120.188712 0 0 0.3311890.3311890.9435640.943564 3.5194823.5194824.9254914.92549110.84011710.840117 第一次迭代得B B2 2 R R1 1Q Q1 1 0.3817390.3817391.0916271.091627 2.3106532.310653 0 0 0.4874950.4874951.3888831.388883 2.9920322.992032 1.00038531.000385312.01339212.013392 重复上述过程， 迭代 11 次得B B1212 0.0074960.0074962.0046952.0046951.9419711.941971 0 0 0.0003250.0003250.9998950.999895 1 1 2.992032,2.992032, 2 2 2.004695,2.004695, 3 3 0.9998950.999895精确值3,2,13,2,1从得出的结果我们可以看出其结果还是很接近精确值的， 其实我们完全可用 Matlab 实现上述计算。用海森伯格 QR 基本算法求矩阵全部特征值 Matlab 得到运算结果为：源程序如下：function l = hessqrtz(A,M)%海森伯格 QR 基本算法求矩阵全部特征值9 / 10普通矩阵特征值的 QR 算法%已知矩阵:A%迭代步数:M%求得的矩阵特征值:lA=5 -3 2;6 -4 4;4 -4 5;M=11;A = hess(A);for(i=1:M)q,r=qr(A);A = r*q;l = diag(A);End而用 QR 基本算法求矩阵全部特征值的运行结果如下：其 Matlab 源程序如下：function l = qrtz(A,M)%QR 基本算法求矩阵全部特征值%已知矩阵:A%迭代步数:M%求得的矩阵特征值:lA=5 -3 2;6 -4 4;4 -4 5;M=11;for(i=1:M)%M 为迭代步数q,r=qr(A);A = r*q;l = diag(A);end由以上两种运算结果的对比可以看出结果基本一致，而且和精确值相当接近，由此可以得出两种方法均可行。参考文献参考文献1 袁慰平、孙志忠、吴洪伟、闻震出.计算方法与实习.东南大学出版社，20052 李海涛、邓樱等MATLAB 程序设计教程M北京：高等教育出版社，200410 / 10