(完整word版)第9章矩阵特征值的数值解法

1、1第 9 章 矩阵特征值的数值解法9.1 引言矩阵特征值问题有广泛的应用背景.例如动力系统和结构系统中的振动问题、电力系统的静态稳定分析上、工程设计中的某些临界值的确定等，都归结为矩阵特征值问题.数学中诸如方阵的对角化及解微分方程组等问题，都要用到特征值的理论.本章介绍n阶实矩阵A Rn n的特征值与特征向量的数值解法 .定义 9.1.1 已知n阶实矩阵A二(aij)Rn n，如果存在常数和非零向量x,使AxVx或(A I )x=0(9.1.1)那么称为A的特征值(eigenvalue)，x为A的相应于的特征向量(eigenvector).多项式称为特征多项式(characteristic p

2、olynomial)，det(A - I )=0称为特征方程( (characteristic equation).注式(9.1.3)是以，为未知量的一元n次代数方程, 多项式.显然，A的特征值就是特征方程(9.1.3)的根.特征方程(9.1.3)在复数范围内恒有解, 其个数为方程的次数(重根按重数计算)，因此n阶矩阵A在复数范围内有n个特征值.除特殊情况(如n =2,3或A为上(下)三角矩阵)外，一般不通过直接求解特征方程(9.1.3)来求A的特征值，原因是这样的算法往往不稳定.在计算上常用的方法是幕法与反幕法和相似变换方法.本章只介绍求矩阵特征值与特征向量的这两种基本方法.为此将一些特征值

3、和特征向量的性质列在此处定理 9.1.2 设n阶方阵A =(aj)n n的特征值为1,2,川，n，那么(1)12Vn二知a?271 - ann；(2)i 2丨In= detA.定理 9.1.3 如果是方阵A的特征值，那么(1)k是Ak的特征值，其中k是正整数；(2)当A是非奇异阵时，丄是A 的特征值.(3)Pn( )是Pn( A)的特征值，其中Pn(x)是多项式Pn(x)二a盼a?x2丨1|anXn.定义 9.1.4 设A, B都是n阶方阵.若有n阶非奇异阵P，使得P二AP = B，则称矩 阵A与Bai2Pn( )二det(A- I)二a2ia22 -IIIIIIainIa2n(9.1.2)a

4、n1an2III(9.1.3)pn( )二det( A釘)是的n次2相似(similar)，P二AP称为对A进行相似变换(similarity transformation) ，P称 为相似变换矩阵( (similarity transformation matrix).定理 9.1.5 若矩阵A与B相似，则A与B的特征值相同.定理 9.1.6 如果A是n阶正交矩阵，那么AJ= AT，且det A= 1或-1;(2)若y = Ax，贝 Vy2=x2,即xT= yTy.定理 9.1.7 设A是任意n阶实对称矩阵，则(1)A的特征值都是实数；(2)A有n个线性无关的特征向量定理 9.1.8 设A是

5、任意n阶实对称矩阵，则必存在n阶正交矩阵P，使得P AP =PTAP=上，其中：二diag(n)是以A的n个特征值仆匕,川，n为对角元素的对角矩阵定理 9.1.9 (圆盘定理)矩阵A二(aj)n n的任意一个特征值至少位于复平面上的几个圆 盘nDi=z|z-aii,i=1,2，山,n,j=, j HJ中的一个圆盘上。9.2 幂法与反幂法9.2.1幕法及其加速9.2.1.1 幕法幕法是计算矩阵按模最大特征值(largest eigenvalue in magnitude)及相应特征向量的迭代法该方法稍加修改，也可用来确定其他特征值幕法的一个很有用的特性是：它不仅可以求特征值，而且可以求相应的特征

6、向量 实际上，幕法经常用来求通过其他方法确定的特征 值的特征向量下面探讨幕法的具体过程设矩阵A Rnn的n个特征值满足人|丸2即乜|兰川公卩 n|启0，(921)且有相应的n个线性无关的特征向量x2l(,xn，则x2l(,xn构成n维向量空间Rnfn1的一组基,因此Rn= ?zz = H%儿，ER, i=1,2|,nj.3在Rn中选取某个满足- 0的非零向量nZo - 7 JXi.i4用矩阵A同时左乘上式两边，得nnAZo:-i AXi八：iXi.i=1i=1再用矩阵A左乘上式两边，得n2、2A Zo二 :i iXi.i丄这样继续下去，一般地有nAkZo - 7、，k =d,2,|）|.i二k

7、记Zk= AZk=A Zo, k =1,2,川，则由式(922)得n kZk1 Xi.式（9.2.4）表明随着k的增大，序列:丰越来越接近A的对应于特征值1的特征向量Xi的：1倍，由此可确定对应于1的特征向量X1.当k充分大时，可得X1的近似值.上述收敛速度取决于比值.事实上，由式（9.2.3）知,再由式（9.2.1）得(9.2.2)nn=E昭,=雪1X,吃i人ii=2l入k=1,2,川.(9.2.3)由假设（921），结合式（923），得(9.2.4)于是对充分大的k有(9.2.5)Zkik(9.2.6)3(9.2.7)4结合式（9.2.6）和式（9.2.7）知，序列 Q/扎：收敛速度取决于

8、比值|妇/人5F 面计算i.由式（923）知k1k1Zk 1 AZk= AZQ=1当 k 充分大时,zk+&霜妆X.结合式（9.2.5），得这表明两个相邻向量大体上只差一个常数倍，这个倍数就是A的按模最大特征值.记/和(n)；血若Zk, Zk, Zk,则有即两个相邻的迭代向量所有对应分量的比值收敛到1.定义 9.2.1 上述由已知非零向量ZQ及矩阵A的乘幕Ak构造向量序列zk来计算A的 按模最大特征值1及相应特征向量的方法称为幕法（power method），其收敛速度由比值丫 = |為/人|来确定，戈越小，收敛越快.注 由幕法的迭代过程（9.2.3）容易看出，如果 卩 q：1 （或1

9、），那么迭代向量zk的 各个非零的分量将随着k-:而趋于无穷（或趋于零），这样在计算机上实现时就可能上溢（或下溢）.为了克服这个缺点，需将每步迭代向量进行规范化：记yk二Azk厂yk1), yk2),|, ykn)若存在yk的某个分量ykjQ），满足y：jQ）二徑務ykj），则记max（ yk） =ykjQ）.将yk规范化yk；max( yj,这样就把zk的分量全部控制在-1,1中.例如，设yk=（-2, 3, Q, -5, 1）T，因为yk的所有分量中，绝对值最大的的是-5，所以max（yQ _ -5，故Zk二y max（yQ =（Q.4, -Q.6, Q, 1, -Q.2）T.综上所述，得

10、到下列算法：算法 9.2.1 （幕法）设A是n阶实矩阵，取初始向量ZQ，Rn，通常取ZQ=（1,1,| l|,1）T， 其迭代过程是：对k =1,2,川，有yk二Amk=max( yQ,Xii =21,2,川，n，(9.2.8)(9.2.9)(j)6定理 9.2.1 对式（9.2.9）中的序列 和有7其收敛速度由 Y=再/人确定. 证明由迭代过程(9.2.9)知Zk二y/mk二AZk j/mAykjmkmk二A=|( = AkZo | min其中Z0冷Xi.若1 =0，则由(9.2.3)知：AkZ0二1kid(9.2.11)得nna-x1 7aii =2于是4Hkmax x1(9.2.10)n

11、j 、a1x書ai-1k1Ximaxmk二max( yQ二ii =2pmm/nk1 Xi“丿Jaii =2na1x1(9.2.15)Xi(9.2.16)4. mg-nk1ct1%+ 迟aixi,代入式-i-宀丿J(nZn、/k、max口兀+瓦iXi =2严J(9.2.12)Ximax x1(9213)AykjA2Z2kA Zoyk= Azk 1kmax( y) max(AZk/) max(A %)二Ak矶max( Akz).(9.2.14)xi)(9.2.11)Xir 為亠二i=28由式(926)和式(927)知：上述收敛速度由 Y =卩吃/州确定.证毕!例 9.2.1 用幕法求方阵12 3A

12、= 213335_按模最大特征值及相应的特征向量，要求mk-mkc10.解选取初始向量zo=(1,1,1)T，按式(929)迭代，结果见表 9.2.1.表 9.2.1kZkmkg -叫1(0.545455, 0.545455, 1)T112(0.560441, 0.560441, 1)T8.27272.72733(0.559787, 0.559787, 1)T8.36279汇1。24(0.559818, 0.559818, 1)T8.358740乂5(0.559817, 0.559817, 1)T8.35890.2X10，因此，所求按模最大特征值8.3589，相应特征向量X - (0.5598

13、17, 0.559817, 1)T.事实上，A的按模最大特征值.198.358898(，相应特征向量为=(0.55981649川，0.55981649 I, 1)T,故所得结果具有较高的精度.9.2.1.2 幕法的加速从上面的讨论可知，由幕法求按模最大特征值，可归结为求数列mk的极限值，其收敛速度由Y=卩2/卅确定.当丫=彼/扎1接近1时，收敛速度相当缓慢.为了提高收敛速度，可以利用外推法进行加速因为序列mk的收敛速度由丫=卩2/州确定，所以若mk收敛，当k充分大时，则有mk-打=0，或改写为-1、kmu*,其中c是与k无关的常数由此可得(9.2.17)这表明幕法是线性收敛的由式(9217)知

14、9g 1九1mk九1mk 2 1mk 1 - 1由上式解出1，并记为mk 2，即机22mk 2mk - mk 1(mk1k)2,(9.2.18)mk 2-2mk 1mkmk 2- 2叫1mk这就是计算按模最大特征值的加速公式将上面的分析归结为如下算法：算法 9.2.2 (幕法的加速) )设A是n阶实矩阵，给定非零初始向量ZoRn，通常取Z =(1,1,川,1T.对k=1,2，用迭代式yk= Azkj, mk=max( yk), Zk= yk/mk求出mi, m2及Z|, Z2.再对k二3,4, |，迭代过程为Yk=Azki,mk=max( y)(mik-m2)2(9219)削 _mk _2 -

15、,耳二yjk.当mk-mk；：，(；是预先给定的精度)时，迭代结束，并计算Zk；否则继续迭代，直至满足迭代停止条件mk- mk有关加速收敛技术，读者请参考文献11.8.2.2反幕法及其加速反幕法是计算矩阵按模最小特征值及相应特征向量的迭代法，其基本思想是：设矩阵A Rnn非奇异，用其逆矩阵A代替A,矩阵A的按模最小特征值就是矩阵A 的按模最大特征值这样用A4代替A做幕法，即可求出A的按模最大特征值，也就是矩阵A的 按模最小特征值；这种方法称为反幕法(inverse power method).、d1因为矩阵A非奇异，所以由AXj = Xj可知：A為xi.这说明：如果A的特征值满足人屮2忙川屮n

16、-1,|人0，(9.2.20)10那么A的特征值满足且对A的应于特征值 i的特征向量xi也是A的对应于特征值 1 i的特征向量.由上述分析知：对AJ应用幕法求按模最大的特征值1 n及相应的特征向量xn，就是求A的按模最小的特征值，n及相应的特征向量Xn.算法 9.2.3 (反幕法) )任取初始非零向量Rn，通常取Zo=(1,1,|,1)T.为了避免求A，对k =1,21，将迭代过程(929)改写为：Ayk = Zj“ mk=max( yj(9.2.22)k=ymk.仿定理 9.2.1 的证明，可得：定理 9.2.2 对式(9.2.22)中的序列Uyk 厂v.9.2.3原点平移法为了提高收敛速度

17、，下面介绍加速收敛的原点平移法设矩阵B = A - pl，其中p是一个待定的常数，A与B除主对角线上的元素外，其他兀素相同.设A的特征值为1, 2I, n，则B的特征值为r-P，2Pl(,n-P， 且A与B的特征向量相同.9.2.3.1 原点平移法的幕法设1是A的按模最大的特征值，选择p，使I2(9221)11-p -p X-p , i =3,4,|,n,及12对B应用反幕法，得：算法 9.2.5 对k=1,2,Hl，有,一2-P入一P-2 (9.2.24)对B应用幕法，得算法 9.2.4 对k =1,2川|，有yk=(A- pi)zkJ, mk=max( yjZk= yk；m,(9.2.25

18、)lim m% - p,lim Zk-,k_：:- max x1其收敛速度由(2- P)：(1-P)确定由式(9.2.24)知：在计算B的按模最大特征值，-p 的过程(9.2.25)中,加速；算法 9.2.4 又称为原点平移下的幕法 (power method with shift).923.2 原点平移下的反幕法设n是A的按模最小的特征值，选择p，使(9.2.26)收敛速度得到入 -P|V人4-P|半厂P，iT，2,|H, n2.及(9.2.27)人- Pi.人n J一P%(9.2.28)若矩阵B = A - pI可逆，贝V B4的特征值为1- P且有1,川，-,几2P州-P1,21,2川,

19、n2P盒i一p(9.2.29)13(A - pl)yk二zkJ,5k=max(yk),吕=yjmk,其收敛速度由(- p). ( nd - P)确定.的过程(9.2.30)中，收敛速度得n_ P到加速.算法 9.2.5 又称为 原点平移下的反幕法( (inverse power method with shift).定义 9.2.8 原点平移下的幕法与原点平移下的反幕法统称为原点平移法.注有的资料上的原点平移法专指原点平移下的反幕法；而有的资料上的反幕法指的就是原点平移下的反幕法.原点平移下的反幕法(算法 9.2.7)的主要应用是：已知矩阵的近似特征值后，求矩阵的特征向量.其主要思想是：若已知

20、A的特征值 韦的近似特征值为 叱，则A-*ml的特征值就是- .m(i =1,2j|,m)，其中i(i =1,2| ,m)是A的特征值.而按模最小的特征 值是- $m，相应的特征向量与A的特征向量相同.利用公式(9.2.30)可求出曲，皐,并且有只要近似值 叱适当好，收敛速度就很快，往往迭代几次就能得到满意的结果.例 9.2.2 已知特征值的近似值：Y 二-0.3589，用原点平移下的反幕法求方阵12 3A= 213335一得对应特征值的特征向量.解取p = -0.3589，对矩阵(9230)kmmkkimzk二Ximax x-i (9.2.31)由式(9.2.28)知：在计算B 的按模最大特

21、征值nfmHki休)确定.于是当k充分大时，可取其收敛速度由1141.358923A pl = A+ 0.35891 =21.35893335.3589迭代公式(9.2.30)中的yk是通过解方程组(A - pI ) yk二Zkj100330.6667100-0.64110.4530-11 -00P( A -1.26791) y Pz，即LUyPz.令Lvk二Pzk4, Uyk二Vkj选取Z，使Uy LPz=(1,1,1)T，得y1-(290929.45, 290927.56, -325732.90)T，mb = max( yj = -325732.90，乙=比伽=(-0.8932, -0.8

22、931, 1)T.1由Uy?= LPR得y2=(-845418.49, -845418.49, 946558.42)T，m2= max( y2) = 946558.42，Z2二y2.mi2=(-0.8932, -0.8932, 1)T.由于z,与Z2的对应分量几乎相等，故A的特征值为1 10.35890.35889894354117,m2946558.42相应的特征向量为x =z2=(-0.8932, -0.8932, 1)T.而矩阵A的一个特征值为=4 - 19二0.358898943540674 |，相应的特征向量为(-0.89315,求得的.为了节省工作量，可先将A _ pI进行 LU

23、分解.在 LU 分解中尽量避免较小的Urr当除数，通常可以先对矩阵A-pI的行进行调换后再分解.为此，我们可用0 0 1P = 010乘A- pl后再进行 LU 分解，即10 00P (A pl )= 0L111.358921.358933*_33=25.35891.358935.35891.35893=LU5.3589-0.5726-3.07 10上15-0.89315, 1)T，由此可见得到的结果具有较高的精度.9.3 QR 算法上一节我们介绍了求矩阵特征值的幕法和反幕法.幕法主要用来求矩阵的按模最大特征值，而反幕法主要用于求特征向量本节将介绍幕法的推广和变形一一 QR 算法，它是求一般中

24、小型矩阵全部特征值的最有效的方法之一，其基本思想就是利用矩阵的QR 分解矩阵A：= Rn n的QR分解就是：用 Householder 变换将矩阵A分解成正交矩阵Q与上三角矩 阵R的乘积，即A二QR.下面首先介绍 Householder 变换.9.3.1豪斯霍尔德变换定义 9.3.1 设B = (bj)nn是n阶方阵，若当i1时，bj=0，则称矩阵B为上Hessenberg 矩阵(Hessenberg matrix)，又称为准上三角矩阵，它的一般形式为川bm!川b2nIIIb3n.(9.3.1)：:时，Af 收敛到一个以 A 的特征值i(i=1,2l(, n)为主对角线元素的上 三 角矩阵.首

25、先介绍用 QR 分解( (QR decomposition 或 QR factorization);即 用Householder 变换将矩阵A分解成正交矩阵Q与上三角矩阵R的乘积，即A = QR.9.3.2.1 QR 分解H2= I- JU2U2011由此得正交矩阵10 00 1I2101 00H2 -00-0.99010.1415J000.14150.9899-1-300-32.3333-0.47140=Q2A2Q20-0.47141.5733-1.5000-00-1.50000.50001000 -0.66670.23570.7070使得QTAQ.02.3333-0.4714-0.4714

26、1.1667-1.5000-1.50000.50000-0.66670.2357-0.70710-0.3333-0.94290.000122算法 9.3.2(QR 分解) )23第 1 步 记A的第 1 列为x(a1(l),a211)|,ani1)T,A = A=佝(叫n利用式(934):fn2=sg询(1) Q (a )Z 丿3 = %其中e =(1,0川,0)TERn,二严1(匚1事),H1= Ir1Ju1u1r,构造出的H1是n阶对称正交矩阵，使得H1 &“丁待，从而有于是有二Hn4代厂Hn4Hn代二HnnJ!( (H1A二H.斗HH 4令R =An,Q二H1H2山Hn4，其中Q

27、是对称正交矩阵，则R = QA.因为Q是对称正交1汀10A2= H1A-)= *-0af?III盘a22)FIIIa22nIFan2IIIr(2)ann第 2 步记X2=(a2?, a32),|,an?)T，同理可构造出n -1阶对称正交矩阵H2，使得H2x2=2e2，其中二2- sgn(a22) |二. 二n皿2，ai2i=2丿e2=(1,0川|,0)TR巴1)若记H2 =山,它仍是对称正交矩阵，于是有IH2丿-CT1a；? anIH06a23)IHa23)=H2A2=0+0(3)a33IHa3n+0r0(3)an3IHf(3)ann如此继续下去，直到完成第n -1步后，得到上三角矩阵24矩

28、阵，所以得注 1 若A非奇异，则上三角矩阵R也非奇异，从而R的主对角线元素不为零例 9.3.3 求矩阵1 0 -1A= 214-2 30的 QR 分解A =QR，并使矩阵R的主对角线上的元素都是正数。解对A运用算法 9.3.2.第 1 步记A二A,洛=(1,2, -2)T，则 sig n( 1)一1222(-2)2= 3,P1=r1(耳十玄们)=3(3十1) =12,(3, 0, 0)T=(4, 2,-2)T,e(1,0, 0)T,-3A2=H1A1=0.0第 2 步 记x(5 3, 73)T，二2二sign(5 3K (5 3)2(7 3)2二2.86744,鶴二6(6 a2?) =2.86

29、744(2.86744 - 53) =13.0013,业二X2二2& =(5.3, 7 3)T(2.86744, 0)(4.53411, 2.33333)T,-1!|4.53411(4.53411,2.33333)13.00132.33333/-0.58124 -0.81373-0.813730.58124-注 2 若要求R的主对角线元素取正数，则A的QR分解是唯一的5 二 二恂=(1, 2,-2)T1TH I u1u1100101112412-2 J(4, 2, -2)=1-1-2_2-221211,2-7 310.3.2 31 010 1一2二I -次Uj253 1.33333A3

30、= H2A2= 0-2.867440 0-10 0为了使R的主对角线上的元素都是正数，取H3=0-10，显然H3是正交阵,0 0 -1且3-1.33333 2.33333A4= H3A3=02.867442.47995.卫02.32494一令3-1.333332.33333R =傀=02.867442.47995，002.32494 _0.333330.15499Q=H1H2H3=0.666670.65874-0.666670.736239.322 QR 算法算法 9.3.3（QR 算法）第 1 步 令A二A，利用算法 9.3.2 将A1进行QR分解，得A二QR，其中Q“为正交矩阵，&

31、为上三角矩阵；然后将Q1与R1逆序相乘，得A?二R1Q1.1 1因为艮二Q1 A，所以有A：二RQ二Q1AA，即A2与A相似.第 2 步 以A?代替A1，再作QR分解，得A?二Q2R2，其中Q2为正交矩阵，R2为.-100 *10!00-0.58124-0.81373H2一0-0.813730.58124H2-2.33333-2.47995-2.32494-0.929980.34874-0.11625了解了QR分解后，下面介绍QR算法(QR algorithm).26上三角矩阵；再将Q2与R2逆序相乘，并记27A3= R2Q2.1 1因为R2=Q2_A2，所以A3= R2Q2= Q2A2Q2，

32、即A3与A2相似.依此类推，可得QR算法公式：对k =1,21,(9.3.10)因为Rk 1- Qk丄Ak 1,所以Ak= RkjQk J- Qk 1AkjQk,即Ak与Ak J相似.故序列：Ak !这里，我们不加证明地给出 QR 算法收敛的充分条件:定理 9.3.9 (QR 算法的收敛性) )设A二(aj- Rnn,是由 QR 算法产生的矩阵序 列，其中代=(a(k).若(1) A 二A的特征值打(i =1,2,|(, n)满足入卜|码|卄九 n| 0 ；A =P二DP,其中D二diag1,，2ln)，且P有三角分解P = LU(L是单位F 三角矩阵，U 是上三角矩阵)，则即代攵敛到一个以A

33、的特征值打(i h,2l(, n)为主对角线元素的上三角矩阵推论 若矩阵A Rn n是对称矩阵，且满足定理9.3.9 中的条件，则由 QR 算法产生的矩阵序列 冲收敛到对角矩阵D = diag( 1,M,.9.3.3带原点位移的QR算法QR 算法收敛速度是线性的，需要提高收敛速度.经分析知：定理 9.3.9 中lim an?kSC1，当n越小，收敛速度越快这启发我们把 9.2 节介绍的原点平移法加速幕法和反幕法的思想运用到QR 算法的收敛加速，这就是下面将要介绍的带原点位移的 QR 算法(QR algorithm with shift).因为求上 Hessenberg 矩阵的特征值比一般的方阵

34、简单，所以在下面的算法中首先将所考察的矩阵化成上 用位移加速方法进行加速lim a(k)k_ -0, i j,i, i = j,(9311)的速度依赖于比值Yn=丸nHessenberg 矩阵，然后再2829其中；是预先给定的精度.算法 9.3.4(带原点位移的 QR 算法) )第 1 步 将矩阵A化成上 Hessenberg 矩阵A1.第 2 步 对k =1,2川 I，给定原点位移数列d，进行迭代加速A f I gRk(QR分解),.Ak 1 =RkQ -Qk 1 & 1Sk1，(9312)并称由Ak到Ak“的变换为带原点位移的 QR 变换( (QR transformation w

35、ith shift)，且记=(a(k)Rn.其中g 的选取方法主要有(k)(a)取Sk= ann；(b)取Sk为Ak的右下角主子矩阵a(k)n d, nVa(k)n,nd(k)的 2 个特征值中最靠近ann的那一个注 1 由算法 9.3.4 产生的Ak 1与Ak相似，即Ak 1二QkAkQk二 QjARQk;且它们都 是上HeSSenberg 矩阵.注 2 计算实践表明，取方法(b)中SR的算法比取方法(a)中SR的算法收敛速度更快.特别 是对称矩阵，带方法(b)中SR位移的 QR 算法是无条件收敛的，且收敛阶为3.注 3 在迭代过程中,当ankn：-o时，由定理 9.3.4 知：anR)的近

36、似特征值.而lim a? = n的速度依赖于比值nk)二k_ i:收敛速度越快；因此平移值取sk=ank)显然是一注 4 在迭代过程中，当ani/肚0时，取a(k)an,na(k)：-n，因此可取(n-Sk)；ann)作为Ank)越小，的 2 个特征值作为A的特征值-n4, n的近似值.注 5 判别an仁和諜心约等于 0 的标准可取为,antz (k)min、an 4,nr,antz30注 6 当求得矩阵A的特征值n 或n4, n 后，可以将AR作降阶处理：对 入左上角的31n _1阶或n - 2阶主子矩阵继续运用带原点位移的QR 算法，以求A的其他特征值，这样可以大大节省运算量9.4 Jac

37、obi 方法上一节我们介绍了 Householder 变换将矩阵化为上 Hessenberg 矩阵的方法.如果矩阵是 实对称矩阵，用平面旋转变换将矩阵化为上Hessenberg 矩阵比用 Householder 变换更好.本节介绍的 Jacobi 方法就是这种方法.它主要用于求实对称矩阵的全部特征值和特征向量.Jacobi 方法(Jacobi method)的基本思想是：对实对称矩阵A = (a人n，定存在正交矩阵Q，使1TQ AQ二Q AQ = D,(941)其中D二diag( id川，n)，j (j =1,2|,n)就是矩阵A的特征值，而正交矩阵Q的 第j列就是对应于j的特征向量.由此可见

38、，Jacobi 方法的实质和关键就是找一个正交矩阵Q，将A化为对角矩阵.9.4.1 Jacobi方法定义 9.4.1 设A= (aj)n n是n阶实对称矩阵,(i )G(i, j门)二为旋转矩阵( (rotation matrix)，或 Givens 矩阵(Givens matrix)，简记为Gij.对A进行的变换(943)称为 Givens 旋转变换(Givens rotation).注 1 Give ns 矩阵是在n阶单位矩阵I的第i行第i列、第i行第j列、第j行第i列、n阶矩阵cos：III1HIIII(i)(9.4.2)IIIHI1IIIcos：(j)GijAGT32第j行第j列的交叉

39、的位置上分别换上rii二COST、rij二sin、5 二-sin、5 二cos而33成的注 2 Give ns 矩阵是正交矩阵，变换(943)是正交相似变换注 3 Jacobi 方法就是通过一系列Give ns 旋转变换，把A化为对角矩阵，从而求得特征值及相应的特征向量的方法.因此，Jacobi 方法也称为平面旋转法(plane rotationmethod).下面介绍具体介绍将n阶实对称矩阵化为对角矩阵的Jacobi 方法.设A= (aj)nn是n阶实对称矩阵，记(1)TA1=(aij)n n= GijAGij.(944)因为TTTTAipGjAGj = GjAGj =A1,所以Ai仍是对称

40、矩阵.通过直接计算可得基文斯(J. W. Givens 1910 年 12 月 14 日-993 年 3 月 5 日)是美国数学家，计算机领域的先驱之一a(1)=ancos2日 +ajjsin2日+ 2丙cos日sin, a(1)=aiisin2v a” cos2二-2aijcos：sin二a(1)=a(1)=a“ cos日+ajisin日，i, j,a(1)二a二-aHsin二-aicos,l G, j,alm =aml =aml,m,_ i,j,a =a(i1-2(ajj-aii)sin 28+aq(cos28-sin28).不难看出，A经过Gj作用后，与A相比，只有A1的第i行、第i列、

41、第j行、第j列的 元素发生了变化，而其他元素与A的相同.特别地，若aj =0，由(9.4.5)式中最后的式子知：若取二满足关系式可使=a=0，也就是说，用Gij对A进行变换，可将A的 2 个非主对角线元素aij和aji化为零.Jacobi 方法的一般过程是：记A0二A，选择A的一对最大的非零非主对角线元素aj和aji，使用 Give ns 矩阵Gj对A作正交相似变换得A1，可使A1的这对非零非主对角线元 素a(1)a(1)二0.再选择A1的一对最大的非零非主对角线元素作上述正交相似变换得A2，可使A2的这对非零非主对角线元素化为零如此不断地做下去，可产生一个矩阵序列A = A0, A1JII,

42、 Ak,UI.注 虽然A至多只有n(n_ 1)2寸非零非主对角线元素，但是不能期望通过n(n -1) 2次变换使A对角化.因为每次变换能使一对非零非主对角线元素化为零，例如，(9.4.5)cotR引21 -tan22ta nTt _ Tt,44(9.4.6)34aj和aji化为零但在下一次变换时，它们又可能由零变为非零不过可以证明，如此产生的矩阵序列A = A0, A1J|, Ak,|将趋向于对角矩阵，即Jacobi 方法是收敛的而这个对角矩阵的主对角线元素就是矩阵A的特征值用 Jacobi 方法求矩阵A的特征值的 步骤为：(1)记A二A，在矩阵A中找出按模最大的非主对角线元素aij，取相应的

43、 Give ns 矩阵Gj，记为Gi= Gj；2 2(2)由条件佝)-aii)sin 2门-2aj(cos二- sin R = 0,定出sin,cos.为避免使 用三角函数，令daii _ajjd= 二，2aj,匕川)，n) =limQ：AQk二QTAQ.k_sc因为Q是正交矩阵，所以QT= Q，若记Q = (Q(YQ(川)Q，其中Q(j)(j -1.2JH, n)是Q的第j列，则有A Q Q, D即A(Q(Q,)川Q(n 2=,Q Q) f|(QnD,，得，)(AQAQ，川,AQ(n)=GQ崩Q(2),川，九nQ(n)，亦即AQ一jQ(j), j =1,2,川，n.这说明正交矩阵Q的第j列就

44、是对应于 j 的特征向量,例 9.4.1 利用 Jacobi 方法求矩阵-_1-201A = -2-11013一的全部特征值和特征向量，要求A.的所有非主对角线元素的绝对值：:：；=0.1.解 记A =A= (aj)3 3，因为a# =-2是A中所有非主对角线兀素中绝对值最大的元素，取相应的 Give ns 矩阵G12.因为a =1, a20)= T，所以d =(a1a20)/2a1；)= -0.5，t = t an = sg d/(| )d ) = - 0. 6 1 8 0 34cos ) - (1 t2)2=0.850651, sin -1 co - -0.525731,于是因为0： ：

45、：2n(n1)k 1E(Ao)十n)：1，所以当k.E(Ak.i) 0.故k 1E( A).(9.4.15)j =1,2,)|l,n.证毕!為E(A爲(A山37则A2二QTAQ2.中绝对值最大的元素，取相应的Give ns 矩阵G13.因为a：?=2.236068, a32?=3.134730，则AAQ1.cos日sin 0010.850651-0.52573101-s in日cos 60=0.5257310.8506510001一002.2360680-2.2360680.850651A1 = (aj)3 3.L-0.5257310.850651Q1= Q0G：= IG0.850651-0.

46、5257310.5257310.850651010 ,1用A1代替A因为a2；）=0.850651是A1中所有非主对角线元素中绝对值最大的元素，取相应的Give ns 矩阵G23.因为a2； =-2.236068, a)=3,所以d =(a22)（1）-a33 2a：3)= -3.077683；t =tan - - sgncos：-(1 t2)2=0.987688,sin v - tcos：- -0.156434,_100 1_1000cos日sin日=00.987688-0.156434卫-sin Bcos日00.1564340.987688一2.2360680.082241-2.37079

47、8.0.519258记=A2=(aj)3x3,3.134730Q2=Q1GT二0.850651-0.5257310.5192580.840178-0.1564340.0822420.133071 ,0.987688用A2代替A1，重复上述过程:因为a；?-0.519258是A2中所有非主对角线元素0-0.525731记记=G1,12二.1 d2=-0.158384,0.082241-0.519258G23二记G238所以390.0747990记-2.3707880.034190 = A (aj)3.-0.0341903.3720780.8078510.519258Q3=Q2GJ =-0.422

48、828 0.8401780.410602-0.156434则A3二Q3AQ3.因为A3的所有非主对角线元素的绝对值：:：；=0.1，所以迭代停止.此时A的特征值 的近似值分别为A3的对角元素： 1.998721,2：-2.370788,3:3.37208.相应的特征向量的近似值分别为x1(0.807851, -0.422828, 0.410602) k1(1.967479, -1.029776, 1)T,x2(0.519258, 0.840178, -0.156434)T= k2(-3.319342, -5.370814, 1)T,X3：(-0.278834, 0.339584, 0.8982

49、95)T=&(-0.310404, -0.378032, 1)T,其中匕=0.410602, k2=-0.156434, k3=0.898295.A的特征值的精确值为2,2=丄-旦二-2.37228132川，3=133= 3.37228132川.22 2 2相应的特征向量为花=(2, -1, 1)T,x2=(-3.186140|1|, -5.372281川,1)T,x3=(-0.313859丨1|,0.372281川,1)T.从例 9.4.1 可以看出，即使迭代的次数不是很多，迭代矩阵Ak的所有非对角元素的绝对值并不是很小时，用Jacobi 方法求得的结果精度都比较高，因此它是求实对称

50、矩阵的全部特征值和特征向量的一个较好的方法.d =(a1?-a3l)r/2a1(2)=0.865333；t =tan v - sgn(d)d +Jl+d2)=0.457089,COST - (1 t2) J2=0.909493, sin v -tcosv - 0.415720,于是COST0G13=01-sin日0sin0.9094930 = 0cos日-0.4157200.4157200.9094931.998721G3A2GI =0.074799 0-0.2788340.3395840.898295409.4.3改进的Jacobi方法由于每次旋转变换前选非零非主对角线元素的最大值很费时间，

51、为此介绍如下改进方 法.第一种方法：把非主对角线元素按照行的次序印23,川，ain,a23, a24, I u ,a2n,川，ann依次化为零，称为一次扫描.一次扫描后，前面已化为零的元素可能成为非零元素，需要再次扫描.这一方法称为循环 Jacobi 方法，这种方法的缺点是：一些已经足够小的元素也 作化零处理，浪费了时间.第二种方法：首先对实对称矩阵A计算n二nI2V=|2ZZ ai2，(9.4.16)I yj4)设置阀值Vi=v.n，按 盹仔|(,ain,a23,a24,l(),a2n,)l(,an,n的顺序进行扫描.若a0|v,，则选取旋转矩阵Gj作旋转相似变换将aj和aji化为零；否则让

52、a0过关， 即不进行旋转相似变换将其化为零.因为某些绝对值小于V的元素的绝对值可能在后面的旋转变换中增长，所以应进行多次扫描，直到Ai的所有非零非主对角线元素的绝对值都小于Vi.再设置阀值v2vi；n，重复上述过程，直到达到精度要求，即V为止(其中；是指定的精度).这种方法称为 限值 Jacobi 方法.9.5 算法程序9.5.1幕法%幕法% A是方阵， e是精度， 输出向量 V是A的最大特征值 MaxEig对应的特征向量 fun ctio n PowerMethod(A,e)V = on es(size(A,1), 1);for i=1:2V = A * V;M(i) = norm(Y, I

53、 nf);V = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = A * V;M(i) = norm(Y , Inf);V = Y / M(i);endMaxEig = M(i);disp(sprintf(最大特征值 MaxEig %.6fn 相应的特征向量,MaxEig)V=V ,end%反幂法41例 9.5.1 用幕法求方阵_123A= 213.335_按模最大特征值及相应特征向量，要求|mk：10J.解 在 MATLAB 命令窗口中输入：A=1 2 3;2 1 3; 3 3 5; e=10A-3; PowerMethod(A,e) 回车，

54、输出结果：最大特征值 MaxEig =8.358906相应的特征向量V =0.55980.55981.00009.5.2幕法的加速%幕法的加速%A 是方阵，e 是精度，输出向量V 是 A 的最大特征值 MaxEig 对应的特征向量fun ctio n Acc_ PowerMethod (A, e)V = on es(size(A,1), 1);for i=1:2V = A * V;M(i) = norm(Y, I nf);V = Y / M(i);endfor i=3:4Y = A * V;M(i) = norm(Y , Inf);M2(i) = M(i-2) - (M(i-1) - M(i-

55、2)F2 / ( M(i)-2*M(i-1)+M(i-2);V=Y/M2(i);endwhile abs(M2(i)-M2(i-1) ei = i + 1;Y = A * V;M(i) = norm(Y , Inf);M2(i) = M(i-2) - (M(i-1) - M(i-2)F2 / ( M(i)-2*M(i-1)+M(i-2);Y = Y / M2(i);endMaxEig = M2(i);disp(sprintf(最大特征值 MaxEig %.6fn 相应的特征向量,MaxEig)V,end9.5.3反幕法% A 是方阵，e 是精度，输出向量 V 是 A 的最大特征值 MaxEig

56、 对应的特征向量 fun ctio n InversePower(A, e)V = on es(size(A,1), 1);L U = lu(A);for i=1:242Y = U (LV);M(i) = norm(Y, I nf);Vector = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = U (LV);M(i) = norm(Y , Inf);V = Y / M(i);endMi nEig = 1/M(i);disp(sprintf(该矩阵按模最小的特征值是：%.6fn 相应的特征向量是：,MinEig)V,end9.5.4原点平移下

57、的反幕法%原点平移下的反幕法% A 是方阵，e 是精度，%输出向量 V 是 A 的特征值 Eig 对应的特征向量，Appr 是特征值 Eig 的近似值fun ctio n Tranln versePower(A, Appr, e)V = on es(size(A,1), 1);B = A - Appr * eye( size(A,1);L U = lu(B);for i=1:2V = U (LV);M(i) = norm(Y, I nf);V = Y / M(i);endwhile abs(M(i)-M(i-1) ei = i + 1;Y = U (LV);M(i) = norm(Y , In

58、f);V = Y / M(i);endEig = Appr + 1/M(i);disp(sprintf(该矩阵的特征值是：%.6fn 相应的特征向量是：,Eig)V=V,end例 9.5.2 已知特征值的近似值：Y 二0.3589，用原点平移下的反幕法求方阵123A= 21333 5的对应特征值的特征向量.解 在 MATLAB 命令窗口中输入：A=1 2 3;2 1 3; 3 3 5; Appr=-0.3589; e=10A-3; TranlnversePower(A,Appr,e) 输出的结果是：该矩阵的特征值是：-0.358899相应的特征向量是：例 9.5.3 设矩阵43V =0.893

59、10.8931-1.00009.5.5 Householder变换%Householder 变换%参数 X 是列向量function Result = Householder ( X ) d = sig n( X(1) * norm(X, 2);e1 = eye( size(X,1), 1 );u = X + d*e1;p = d * ( d+X(1);H = eye(size(X,1) - pA-1 * (u*u);Result = H;end9.5.6用Householder变换将一个n阶方阵化为上Hessenberg矩阵%用 Householder 变换将一个 n 阶方阵化为上 Hess

60、enberg 矩阵%输出的 Q 是所求的正交矩阵，% 输出的UpH是用Householder变换将n阶方阵 A化成的上 Hessenberg矩阵 function UpHesse nberg( A )N = size(A, 1); |Ai = A;Q = eye(N);for i=2:N-1A = Ai(i-1:N, i-1:N);X = A(2:size(A,1), 1); HH = Householder ( X );Qi = eye(i-1) zeros(i-1, N-i+1); zeros(N-i+1, i-1) H;Q = Q * Qi;Ai = Qi * Ai * Qi;endQ, _UpH=Ai,end4412 122 211A