第七章_矩阵特征值与特征向量计算

1、 物理，力学和工程技术中很多问题在数学上归结为求矩阵特征值向量。即下面的数学问题： 1.已知： 求代数方程 的根。 称为A的特征多项式特征多项式。上式展开即为: 的n个解称为A的特征值特征值。 nnija 0det 0.111nnnnccc2.设 为A的特征值，要求相应的齐次方程组 的非零解。即求 的非零解。 此非零解x称为矩阵A的对应于 的特征向量特征向量。0 x xx 回顾几个基本代数结论:定理定理1 若 为A的特征值，则有: (1) (2)nii.1trniiinii11 nA.det21定理定理3. (Gerschgorins定理)设 则A的每一个特征值，必 属于下述某个圆盘之中: 且

2、如果一个特征向量的第i个分量绝对值最大,则对应的特征值一定属 于第i个圆盘中. nnija1,1 . . . . . .ni ii jijiaain定义定义1:设A为n阶实对称矩阵，对于任一非零向量x称称 为对应于向量为对应于向量x的的Rayleigh商商. xxxxxR,(1)A与B有相同的特征值；(2)若x为B 的一个特征向量，则Tx为A的特征向量。定理定理2.设A与B为相似矩阵(即存在T， ， )，则:0detATT1 关于矩阵A的特征值计算问题，当n=2,3时，可以按行列式展开的方法求 的根。但当n较大时，再按行列式展开法先求 的系数，再求 的根，工作量是非常大的，在实际应用中这种方法

3、是不切实际的。 0 (1)(2) (3)1,0,nnx xxRxx x1,0,m axnxRxA x xx xmin,.,0,nnx xxRxx x定理定理4.设 为对称阵，特征值 ,对应的 特征向量 ,组成规范化正交组，即: 则:nnRn.2112,.,nxxx1,0ijijijxxij 在一些工程实际问题中，通常只要求出矩阵按模最大的特征值按模最大的特征值(称之为A的主要特征值主要特征值)和相应的特征向量，对于这种特征值问题应用幂法是较合适的。 幂法幂法是一种计算实矩阵 A 的主特征值的迭代法，最大优点是方法简单，对稀疏矩阵较合适，但有时收敛速度很慢。 设实矩阵 有一个完全的特征向量组,其

4、特征值为 ， 相应的特征向量为 nnija12,.,n12,.,nxxx 已知A的主特征值是实数,且满足:幂法的基本思想基本思想是：任取一个初始向量 , 由矩阵A构造一向量序列（称之为迭代向量序列迭代向量序列。） : (2)123.n00v 211021010,.,.kkkvAv vAvA vvAvAv10111221111112121.kkkkkknnnknkkiiikiknikiiivAvA vaxaaxa xaxa xax于是：其中：111(4)limkkkva x 从 而 ：112, 3,.,)0,ikink 由 假 设 ：， （故 ： 这表明： 越来越接近A的对应于 的特征向量 （数

5、量 不计）或者k充分大时，即迭代向量 为 的特征向量的近似向量(除一个因子外)。1kkv11x1a111(5)kkvaxkv1而由假设, 可表示为特征向量 的线性组合： nxxx.21011221., (0 ),(3)nnva xa xa xa设0v下面考查主特征值 的计算。用 表示 的第 个分量, 则：111111111()( 6 )( 7 )l i miiikkiiikkikikkvaxvaxvv故1( ) ikvkvi 这表明这表明：相邻两次迭代向量的对应分量的比值收敛到主特征值。 这种利用已知非零向量 及矩阵A的乘幂 构造向量序列 以计算A的主特征值 （由(7)式）及相应的特征向量（利

6、用(5)式）的方法称为幂法幂法。 0vkA kv1由(6)及 表达式可见： 的收敛速度由比值 来确定，r 越小收敛越快。而 时收敛可能很慢。将上述讨论总结一下有结论:k11ikikvv21r211r定理定理5： 有n个线性无关特征向量。主特征值 满足 ,则对任何非零初始向量v ,有(4)、(7)成立。nnRA1n.321 01a 若A的主特征值为实的重根，即 ，且又设A有n个线性无关的特征向量, 对应r个线性无关特征向量为： 12.r12.rrrn112,.,(2 ):rxxx则 由有01111krnkkikiiiiiirvA va xax 从而 (设 )相应的（6）（7）式成立。11limr

7、kiikkiva x 10riiia x 这表明：当A的主特征值为实的重根时，上一定理结论仍正确。 利用幂法计算A的主特征值 及对应的特征向量时,如果 (或 )，迭代向量 的各个不等于零的分量将随 11111kvk而趋向于无穷(或趋于零)，这在计算机上计算时可能产生“溢出”。为克服这个缺点，我们将迭代向量规范化。设 ，其规范化向量为 ，其中max(v)表示向量v的绝对值最大的分量绝对值最大的分量。0v maxvuv 在定理5的条件下幂法可这样进行：任取一初始向量 构造向量序列：（ 为向量） 01aiv01100110 ,m axm ax A vvvAuAvuvA v220022122020,m

8、axmaxmaxA vA vvvAuuAvvA v00100,.m axm axkkkkkkA vA vvuAvA v00v（下面 为向量。） 由(3)式:ix0111121(8)knnkkkiiiiiiiiA vaxa xax111210011121:maxmaxknkiiikikkknkiiiia xaxA vuA va xax则1121111121m axm axkniiiikniiiia xaxxkxa xax 这表明：规范化向量序列收敛到主特征值对应的特征向量。同理可证:111210110111121maxmaxknkiiikikkknkiiiia xaxA vvAva xax则收敛

9、速度取决于比值11112111112maxmaxmaxkniiiikkniiiia xaxvka xax 12r (9) 有: .0010maxkkkkkvuvAuvuv,.2,1k111,maxmaxlimlimkkkkxuvx 定理定理6: 有n个线性无关的特征向量，主特征值 满足： , 则对任意非零初始向量： , 按下述方式构造向量序列：nnRA1n.321001,0vua总结以上讨论有下面结论：解： 取初始向量 ,按(9)迭代5次得到数据如下 表: k (规范化向量) 0 1 1 1 1 1 1 1 0.2143 0.4821 1 12.00 27.00 56.00 2 0.1875

10、0.4483 1 8.357 19.98 44.57 3 0.1860 0.4463 1 8.168 19.60 43.92 4 0.1895 0.4460 1 8.157 19.57 43.88 5 0.1859 0.4460 1 8.156 19.57 43.88 Tuv11100TkuTkv361641593642A例1：用幂法计算下面矩阵的主特征值及对应的特征向量。88.431对应的特征向量为:Tu0000. 14460. 01859. 01故按模特征值为: 应用幂法计算A的主特征值的收敛速度主要由比值 来决定。12r1r但当 接近 时，收敛可能很慢，这时可采用加速方法补救。.)(,为

11、先定的参数引进矩阵ppIAB.,.,.,2121特征向量相同而且的相应的特征值为则的特征设BApppBAnn且使的主特征值为仍使就要选择的主特征值如果需要计算,11BpPA1212pp1BBpA对 应用幂法，使得在计算 的主特征值的过程中得到加速。这种方法称为原点平移法。 对于 的特征值的某种分布，它是十分有效的。:例如. 9 . 0) , 4 , 3 , 2 , 1( ,15)(12rjjaAjnnij比值有特征值1,0,1,2124321的特征值为则取作变换BppIAB1222111:10.92Bprp应用幂法计算 的主特征值的收敛速度1ppAp可见的确能通过选择有利的参数 使幂法得到加速

12、，但设计一个自动选择适当参数 的过程是困难的。下面考虑当 的特征值是实数时，怎样选择 ，而使用幂法计算主特征值得以加速。.,)10(.:1121pppIABpAnnn或的主特征值为为何值则不论的特征值满足设pppxn111:,)(使首先应选择时特征向量及当我们所希望计算211:max,minnppwpp且使收敛速的比值为nnnnnpppprwpppp212*1*1*2*222:,2,)(这时收敛速度的比值为为最小时时则当2n2n002p0pn0pn02ppy.,),10(*2的近似值我们就可以确定能初步估计时且的特征值满足因此当pAn.2,*11得到加速使得应用幂法计算应选择时当希望计算nnn

13、pp的主特征值计算例0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10 . 1:A,15(8):A若直接用幂法求 的主特征值 则迭代次位浮点数 得到Tx) 1 ,6497. 0 ,7483. 0(,5366256. 211:0.75p 下面用原点平移加速法，计算结果如下 取，则：11110.251.000.501.000.250.25:0.50865914 ():2.5365914,:(0.7482,0.6497,1)TpBApxB=A- I=，对B应用幂法，迭代仅十次即有的主特征值 。从而 的主特征值为特征向量为11112.5365258

14、15:1.786258,2.53652588Ap与矩阵 的主特征值比较可知：后者比前者迭代次数少，但精度更高，即提高了收敛速度。若用后一方法也迭代次，则有对应有，已达 位有效数了。ApA原点平移加速方法是一个矩阵变换方法，这种变换容易计算而且又不破坏矩阵 的稀疏性。但 的选择必须依赖于对的特征分布已有大致了解，这在实际中是相当困难的。00112(, )(, ),( , )( , )(.)maxminn nn nxxnx Rx RnAx xAx xx xx xRayleighA由于，则我们有可能将商应用到用幂法计算对称矩阵 的特征值的加速收敛上来。12117:.(,)(6(9),:n nniji

15、jkARxxAuRayleigh定理设为对称矩阵，特征值满足。对应的特征向量满足，应用幂法 公式见定理 中式计算 的主特征值则规范化向量的商给出的较好近似2211(,)() )()(,)kkkkkAuuOuu。 证明见教材:商加速计算过程为因此Rayleigh001( )10(11)max()(,)(,)kkkkkkkkkkvuvAuvuvAuuuu.,近似特征值的特征向量及计算对应于一给定的其特征向量模最小的特征值及反幂法用来计算矩阵按12121,:.:, .,:n nnnARAxxxA设为非奇异矩阵的特征值为相应的特征向量为则的特征值为12111111.,:,.,nnnnxxxx对应的特征

16、向量为.1的特征值问题的按模最大的问题就是计算的按模最小的特征值因此计算AAn:.,1),(11反幂法迭代公式为的按模最小的特征值从而求得的主特征值可求得矩阵称为反幂法应用幂法迭代对于nnAAA.,.2 , 1)max(. :. 011100求得可以通过解方程迭代向量构造向量序列任取初始向量kkkkkkkkuAvvkvvuuAvuv1210018:.0(0) ,:1(1)lim,(2)lim max()max():.nnnkknkkkknnnnAnAuvavuxuvx定理设 有 个线性无关的特征向量， 为非奇异阵且特征值满足，则对任何初始非零向量，由反幂法构造的向量序列满足收敛速度比值为.及其

17、特征向量或求其它特征值点平移来加速迭代过程在反幂法中也可以用原:)(,1的特征值为则逆存在的若pIApIAnxxx,.,:.21对应的特征向量仍然是:,)(1得到反幂法的迭代公式应用幂法现对 pIA00110 ()1,2,. 13.11,2.lim.max()max()kkkknkkkuvvApIukvukpvv初始向量（ ）pppn1.,1,1211()1()(13)jjijjpAppijApIp若 是 的特征值的一个近似值，且即为的主特征值。则可用反幂法计算该特征值及其特征向量。:)()(,)max()()max()(),0(,.,10000)1(01021由此可得结论其中则个线性无关的特

18、征向量有inikiikkkkkkkjniiinnnxpaupIAupIAupIAuupIAupIAvaxauxxxnRA1019:(1,2. )()()0(0):(13) :n niijjiniijikkARnAx inpApIppijua xavu定理设有 个线性无关特征向量， 的特征值及对应的特征向量记为及， 为的近似值，存在，且，并且为给定初始向量， 则 由反幂法迭代公式构造的向量序列，满足11)lim,(2)lim max(),:max()1()maxmax()jkkkkjjjjijkixuvxpppkrvp 即， 收敛速度比值为:()jjjApIpxpr由此定理知 对其中应用反幂法，

19、可用来计算特征向量，只要选择的 是的一个较好的近似且特征值分离情况比较好，一般 很小，常常只须迭代一、二次就能完成特征向量的计算。1()():kkkvApI vuApI反幂法迭代公式中的是通过解方程组求得的，为节省计算工作量，可先得进行三角分解1():kkkkkP ApILUPvLyPuPUvy为某个排列矩阵，于是求 相当于两个三角形方程组为某个排列矩阵。11,2.(14)max()kkkkkkkLyPuUvykvuv从而反幂法公式可写成.:00是较好的按下述方式选择经验表明uv 110(1,1,1,.,1)(15)TUvL Pu1(15)(14)v由回代求解可得， 然后再由迭代求解。3:1.

20、2679(33)(5)210131014AA例 用反幂法求 的对应于近似特征值精确解的特征向量 用 位浮点数进行计算LUpIAPppIA)(:)2679. 1(:分解为其中解将用部分选主元的三角分解0011000101029405. 0007321. 2101732. 11,126807. 07321. 00100013PUL其中:) 1,1 ,1 (1得由TUv TTvvuv)26795. 0,73198. 0, 1 ()max(,)8 .3400,3 .9290,12692(111121223323:(20404,14937,5467.4) ,(1,0.73206, 0.26796):(1

21、,13, 23)(1,0.73205, 0.26796).TTTTLUvPuvuxux由得而所对应的特征向量是可见已是 的相当好的近似了:210121012AA 例 利用逆幂法计算下面矩阵 的按模最小特征值:1002101/21003 2102/31004 3ALUALU 解对 作分解得)()0, 0, 1 (10本例中未进行标准化解取kkkkTyUvvLyv计算结果如下表：111/45/85/819/1619/1667/3205/8121/167/475/32301/37/67/87/421/1623/869/321v2v0v0y1y4v3y2y3v1,2,3)(i /)(3)(4iivv由

22、于 的平均值之倒数为而A的按模最小特征值精确值为可见已获得了较好的特征值。5858. 0221330.585922767 12 23227,)387 14133(三个比值： ，均值为 于是求实对称阵A的特征向量问题就等于寻找正交矩阵P使为对角阵，其主要困难在于寻找正交阵。DTpAp Jacobi方法是用来求实对称矩阵的全部特征值及对应特征向量的一种变换方法。其基本思想是基本思想是：通过一组平面旋转变换（正交相似变换）将对称阵A化为对角阵，得其全部特征值。i)1(ni nnR12(,.)TnPAPdiagD 由代数理论：为对称矩阵，则存在一正交阵P，使且D的对角元就是A的特征值。jvTP的列向量

23、就是A的对应于 的特征向量。j先考察平面上正交变换： ，其中为平面旋转矩阵。则 ，其中：1122PyxyxcossinsincosP11122122TccPAPcc 221111222122221122212112221121cossinsin2sincoscos21()sin2cos22caaacaaaccaaa2cos2sin)(210 2111221221aaaCC使故可选角先看2阶对称阵 能否找到对称阵P使A经过正交相似变换化为对角阵？aaaaA22122111 即 （当 时，可选取 ）从而可使aaatg22111222aa22114),(21diagPAPTcossinsincos(

24、, )ijiijjkkyPxki jyxxyxxyx 记为将上面方法推广，首先在 Rn 中引进平面旋转变换：( , )11cossin11sincos11i jijiPjP平面内的一个旋转变换中称之为,.21.21),(),(xxRyyyyyxxxxxjinTTnjiynjix12, ) ( , ) ( , ) ( , )3TPPPi ii jj jj iPAAijAijAijijPAPT有如下性质：为正交矩阵。和单位矩阵只在（，四个位置上的元素不一样。只改变 第 行元素与第 行元素；AP 只改变 的第列与第 列元素；只改变 的第 行，第 行，第列，第 列元素。2210n nTFFACPRPA

25、PCA定理：为对称阵， 为正交阵，则：( )2222111.2sincos(2)2sincos(3)cossinsincosn nijTijiijjijiijjjjjjijACRPcPAPcaaacaaa定理 ：为对称阵，为一平面旋转阵，则的元素的计算公式为：12.()sin2cos2(4)2ijjijjiiijccaaa 5.cossin ,(, .1,., )(7)cossin ,(, .1,., )(8)6.,( , .,1,., )(9)kikjkikjkjkilklkijki jknki jkncal ki jl kncaacaa第、列元素：由对称性知：（7）、（8）的计算已经由（5

26、）、（6）得到。3.:cossin ,(, .1,)(5)4.cossin ,(, .1,)(6)kiikjkikkjjkikjkiki jknjki jknccaaccaa第 行元素第 行元素：),(. 32. 2),.(. 10,01222222222222),(jikljikACCAAacaaaccaaccPAPccPAPPaRlklkijjjiijjiijkikjkikTjiijTjiijnn的元素有如下关系：与且的非对角元素为：使则可选择一平面旋转阵的一个非对角元素，为为对称阵，：定理 aaijijACAC2222,少了的非对角元素平方和减比的非对角线元素平方和大了的对角线元素平方和

27、增的对角线元素平方和比即)10(22122)()(2)()(aSSaDDijAcijAc知：由定理:Jacobi下面介绍方法 AAA继续这个过程，连续对 实行一系列平面旋转变换，消除非对角线绝对值最大元素，直到 的非对角线元素全化为充分小为止，从而求得的全部近似特征值。222222(1)(1)(1)1,2221222(2)(2),()max012(,),0lkijl kl kTjjiiaaaAjiPAP A Pcc再选的非对角元素中的主元素，如：由定理，又可选择一平面旋转变换使的非对角元素。（前一次消除了的这次可能又变为非零了。）11111111,1111,max,012,)10(lklkAj

28、jiiAPjijiaaaP i jccT1先对非对角线元素中选绝对值最大的元素（称为主元素），如：可设，否则 已对角化了。由定理知：可选择一平面旋转阵，使=AP 的非对角元素113.(1,2,)lim(n nTmmmmmmJacobiAmDRAP A PA定理：（方法的收敛性）为对称阵，对A实行一系列平面旋转变换：则对角阵） , )21200sin2cos202221,2()1541(15)Tijjjiii jijjiijjiijijiijjCtgtgdaPAPaaccccaPaaatg（关于sin 、cos 的计算：由定理知：当时，可选 满足：使的非对角线元素即：选使：即选择使：其中（ ）c

29、ossin,sincos,(1,2,)(14)TmliljlililjljlnRPrrRrrR的计算公式为：，则的两列元素，若记只需计算右乘用初等正交阵可采用累积计算方法对),(,.21jiRRRRRRIPPPRPRPRPRmmTmTmTTTm的近似特征向量。的列向量就是记法收敛性定理可知：由关于特征向量的计算则ADAJacobiRPPPRPPPPPPPTmTmTTTmTmTTmm.,:2121121 1.),(2),(3)(15),2(15)ijijiilliijjtc sdaaaaaa（若则可取将带入并利用可知：)16()(cossin11cos0101)(.1)(2,max. 022sc

30、ttcdddStddStgdtdaaaaaaijjjiilkklijij时当时当则：设.cos,sin,011) 1(01122来计算可利用元素时当由此可见时当选取aaaddijjjiiddtgddtg.11012),(22)(dtgaadtgdtgjjii即时当的特征值。法计算对称阵例：用210121012AJacobi23()()mmlkijARmlkASa 的第 列，第 列两列元素，再利用对称性传递元素。因而不用做 个矩阵乘法。在计算机上，需两组存贮单元以存贮（或）和可用控制迭代终止。)17(0 ccaacaacjiijijjjjjijiiiittA1m中的公式，只需计的11后，再利用定

31、理cossini但如果先 计如PT1mAP1mA1m元素与 计素的非 对非对角线元素中Am于选每选代一次主要工作在 110011d2dS(d)t0a122a22a11d则2)j1,(i1a12选非对角元素：A,A0步第一27071068. 07071068. 07071068. 0207071068. 001.7071068. 021.7071068. 021111112PPAtTctSC366025. 26297630. 006279630. 033250567. 003250567. 06339746. 04597008. 09880738. 05176381. 07071068. 070

32、71068. 06.5031)3, 1(7071068. 0212223131332233113211212222212111113PAPAaacccccccaacaacaTsctdsscttji则类似计算可得：对）计算如（）由（元素的主元素第二步，在中选非对角 JacobiA法是一个求 的全部特征值、特征向量的有效迭代方法，但其计算过程由于平面旋转变换将破坏原矩阵的稀疏带状特征。这是其缺点。2T55453A20.58578790.203825201020.20382523.4140130.0167579010AP A P00.016757902.000198再类似地以为基础迭代 步有：的精

33、确精确特征值A0.585787932.00019823.4140131的特征 值特征值近A于是得近似特征向量。的A的列向量，即得PT5PT4PT3PT2PT1RT5且可逐 步可逐0.585786）221（ 2322，414214。3）221（ 21一对共轭复特征值。的是阶对角块的两个特征值的实特征值，每一个二角块即为一个一阶对为一阶或二阶方阵，每其中的对角块AAmiRii), 2 , 1(我们称这种分块上三交阵为矩阵A的Schur分块上三角阵，上三角阵和对角阵是它的特殊情形。定理7.9并没有解决如何计算全部特征的问题。7.4 QR算法算法7.4.1 化矩阵为化矩阵为Hessenberg形形对于

34、实对称矩阵，可通过正交相似变换约化为对角矩阵。那么，对于一般的实矩阵，通过正交相似变换可约化到什么程度呢？线性代数中有如下结果。，使得存在正交矩阵对于任何矩阵QRAnn,mmmmTRRRRRRAQQ22211211定理定理7.9 （实（实Schur分解定理）分解定理）定义定义7.20ijnnijbRbB的次对角线以下的元素）（若矩阵(ij+1), 则称B为上Hessenberg矩阵，简称Hessenberg形，即B的形状为B是可约的，则（有一个次对角元如若BnkbBkk),110, 1矩阵的特征值问题。问题约简为求解较小的形，可把求解特征值的约的否则是不可约。对于可Hessenberg可以用平

35、面旋转变换化矩阵为Hessenberg形,下面介绍另一种正交变换。为了节省运算工作量，实用的方法是先将矩阵约化为与Schur分块上三角阵很近似的Hessenberg形。定义定义7.3则称设向量, 1,2wRwnTwwIwH2)(（7.4.1）为（初等）镜面反射矩阵，（初等）镜面反射矩阵，或Householder变换矩阵变换矩阵。Houholder矩阵H=H（w）有如下性质： ，。事实上，显然有是对称正交阵，即HHHHHHTT1(1)得知又由12 wwwT。IwwwwwwIHHHTTTT)(442(2)22,xyHxyRxn有记对任何(3) 记S为与w垂直的平面,则几何上x与y=Hx关于平面S对

36、称。事实上,由得知xwwIHxyT)2( 。wxwyxT)(2上式表明向量x-y与w平行，注意到y与x的长度相等，于是x经变换后的象y=Hx是x关于s对称的向量，如图7-1所示。xwyx-y图图7-1对应于性质（2），有下面的定理。，使则有镜面反射矩阵且设HyxyxRyxn,22定理定理7.10得Hx=y。则有令12,)(22wwwIHyxyxwT证证知由yyxxTT。22)()(2)(2yxyxyxyyyxxxxyxTTTTT。yyxxyxxyxyxxxwwxHxTT)()(2222由此可得定理得证。，使得有镜面反射矩阵是对该定理的一个重要应用HxxxxTn0),(21。1eHx（7.4.2

37、）的计算公式为。矩阵其中HexxsignT)0 , 0 , 1 (,)(121。，TuuIHxexu111)(,（7.4.3）值计算的尽量大，从而有利于数使分母的的符号的选取，是为了关于稳定性。（7.4.2）的意义是对向量作消元运算。与平面旋转不同的是，镜面反射变换可成批的消去向量的非零元。，使构造镜面反射对于向量HxT,) 1 , 1 , 5 , 3(。TxxsignHx)0 , 0 , 0 , 1 ()(21例例7.4 解解。Texuxxsignx) 1 , 1 , 5 , 9(, 6)(, 61212）得。按（3 . 4 . 754,1082u。531591535955294599452

38、75411TuuIH（7.4.4），使得存在正交阵对于任何矩阵，QRAnn,。AQQBT定理定理7.11为Hessenberg形。维向量。不含对角线）的第一列对角线以下（为，设1111nAaAA证证111111,12 .4 .7eeaHHn其中，使得阶对称正交阵）可构造根据（。用是对称正交阵，显然。记1111111), 1 ()0 , 0 , 1 (PPPHdiagPRnT作相似变换，易知对11AP。11112PAPA造维向量，那么同利可构角线）列对角线以下（不含对的第为记2222nAa22112222) , 0 , 0 , 1 (,2IReeaHHnnT。记其中，使得阶对称正交阵做相似对为对

39、称正交阵。用显然阶单位向量，为222222),(2APPHIdiagP 变换，有。212223PAPA如此类推，经n-2步对称正交相似变换，得到Hessenberg形矩阵。2211222221nnnnnnPPAPPPPPAPA）。定理得证。，则有（若记4 . 4 . 7,2211nnPPPQAB为，使得存在正交阵对于任何对称矩阵AQQBQRATnn,推论推论7.1对称三对角阵。这时，可以的维数大于的向量的矩阵，第一列被消元形。对于阶数大于231a上述定理7.11的证明是构造的，即可以用镜面反射化矩阵Hessenberg形。此定理可用平面旋转变换来证明，即也可用平面旋转变换化矩阵为Hesenbe

40、rg逐个化为零。个分量开始的非零元素的从第，把连续使用平面旋转变换21a如此类推，最后得到的正交矩阵Q，是平面旋转矩阵的乘积。7.4.2 QR算法及其收敛性算法及其收敛性 QR算法可以用来求任意的非奇异矩阵的全部特征值，是目前计算这类问题最有效的方法之一。它基于对任何实的非奇异矩阵都可以分解为正交阵Q和上三角矩阵R的乘积。定理定理7.2(QR分定理分定理)正交阵与为非奇异矩阵，则存在设nnRA上三角阵R,使得A=QR,且当R的对角元素均取正时,分解是唯一的。证证类似于定理7.11的证明，对矩阵A的左乘一系列正交变换矩阵，可以将A化为上三角形矩阵，因此，可得A的QR分解。下面证明分解的唯一性。设

41、有两种分解2211RQRQA此可得的对角元均为正数。由，而且21RR。11212RRQQT上式左边为正交阵，即。）（）（1112112RRRRT这个式子左边是下三角阵，则右边是上三角阵，所以只能是对角阵。设),(21112ndddRRD，进而，从而故有且则有122, 2 , 1, 0,RRIDnidIDDDiT定理得证。，21QQ 一般按平面旋转变换或镜面反射变换作出的分解A=QR，R的对角元不只要令一定是正的。设),(ijrR ),(22221111nnnnrrrrrrdiagD就是符合的上三角阵，这样，为对角元是为正交阵，QRArRDRQDQii1_定理7.12的唯一QR分解。BAQQBR

42、QBQRAQRAT。这说明，则有，那么令分解，即的设有，得阵。令分解，又可得一新的矩继续作有相同的特征值。对与AAQRBA1如下的算法：。, 2 , 1,1kQRARQAkkkkkk（7.4.5）的方法称为）得到矩阵序列由（kA5 . 4 . 7或称为基本基本QR算法。 QR算法，证证 容易证(1)从它递推得,)()(11121121111kkTkkTkkkTkkQAQQQQAQQQQAQA1_1_1221_kkkkkkkRAQRRRQQQRQ。1_1_1_1_1_1kkkkTkkRQARQAQQ）。定理得证。，即证得（由此推及21111_1_AARQRQkA定理定理7.13QR算法产生的序列

43、 满足:;1kkTkkQAQA。，其中12_21_,RRRRQQQQRQAkkkkkk（1）（2）kA一般情形下，QR算法的收敛性比较复杂。若矩阵序列 的角元均收敛，且严格下三角部分元素均收敛到零，则对求A的特征值而言已经足够了。此时，我们称 基本收敛到上三角阵基本收敛到上三角阵。下面对最简单的情性给出收敛性定理。kA设矩阵 的特征值满足kA定理定理7.14,nnRA, 021n121),(, 2 , 1,XxxxXnixnii的逆可分解为。若矩阵对应特征向量=LU，其中L为单位下三角阵U为上三角阵，则QR算法产生的序列 基本收敛到上三角阵，其对角极限为。niaikiin, 2 , 1,lim

44、)(kAkA更一般地，在一定条件下，由QR算法生成的序列 收敛为Schur分块上三角形，对角块按特征值的模从大到小排列，上述定理是它的特殊情形。当收敛结果为Schur分块上三角形时，序列 的对角块以上的元素以及2阶块的元素不一定收敛，但这不影响求全部特征值。例例7.5 用QR方法求下列矩阵的全部特征值。112112314)2( ,151311313153) 1 (解解先用镜面反射变换化矩阵A为Hessenberg形矩阵 ，然后用平面旋转变换作QR分解进行迭代，生成序列 。（1）的计算结果为1AkA，2062. 14639. 05361. 59738.129284.137282. 23077.

45、40000. 31A，5229. 14656. 09381. 06362. 05511. 18265. 07711.171133.102A，9980. 10012. 06317. 10019. 30000. 02165. 99820.160001. 616A。9999. 10000. 06329. 100001. 30000. 02364. 99712.160000. 623A该矩阵A非对称，从计算结果看，收敛于上三角阵。（2）的计算结果为,0000. 10000. 10000. 10000. 18284. 24142. 18284. 20000. 41A,4000. 04899. 03266

46、. 02667. 17454. 01121. 59379. 13333. 22A,3333. 27456. 07263. 33336. 00002. 06516. 38171. 00003. 225A。9998. 02349. 22366. 29996. 20001. 02374. 29999. 20002. 226A 从计算结果来看，迭代收敛于Schur分块上三角形，对角块分别是1阶和2阶子 ，迭代已接近收敛。阶子矩阵的特征值都是的右下角的和矩阵。事实上，矩阵iAA0000.19999.022625 一般在实际使用QR方法之前，先用镜面反射变换将A化为Hessenberg形矩阵H，然后对H作

47、QR迭代，这样可以大大节省运算工作量。因为上 Hessenberg阵H的 次对角线以下元素均为零，所以用平面旋转变换作QR分解较为方便。 对i = 1,2,.n - 1,依次用平面旋转矩阵J（i , i+1）左乘H，使J（i , i+1）H的第i +1行第i列元素为零。左乘J(i,i+1)后，矩阵H的第i行与第i+1行零元素位置上仍为零，其他行不变。这样，共n-1次左乘正交矩阵后得到上三角阵R。即 =R， =J（n-1,n）J(n-2,n-1)J(1,2)。可以验证 是一个下Hessenberg阵，即U是一个上Hessenberg阵。这样，得到H的QR分解H=UR。在作QR迭代时，下一步计算R

48、U，容易验证RU是一个上Hessenberg阵。以上说明了QR算法保持了H的上Hessenberg结构形式。 HUTUTUT例 7.6求Hessenberg形矩阵 2100322023011525H的特征值。 解解），使得，（），（），（有平面旋转矩阵的次对角线非零元素，根据令43322111JJJHHH 7822. 00002736. 35242. 2005288. 25852. 10381. 203922. 04912. 59612. 10992. 58989. 03962. 0000740. 01761. 09813. 004192. 08804. 01887. 01961. 01038

49、. 01923. 00377. 09806. 0111RUH，逆序相乘，求出和）。然后将求得的，（），（），（其中2111213243HRUJJJUT7031. 03099. 0001294. 38525. 04770. 205361. 15171. 29401. 13997. 02390. 05922. 19508. 56157. 4112URH重复上面的过程，计算11次得0000. 1*1211. 03290. 1*5910. 38789. 1*0000. 412H至此，不难看出，一个特征值是4，另一个特征值是-1，其他两个特征值是方程01211. 03290. 15910. 38789.

50、 1的特征方程为阵。事实上，可以求得矩的根，求得为Hi 21, 020775234上述用QR方法求得的特征值是该特征方程的准确解。7.4.3 带原点位移的带原点位移的QR算法算法 前面我们介绍了在反幂法中应用原点位移的策略，这种思想方法也可用于QR算法。一般我们针对上Hessenberg矩阵讨论QR算法，并且假设每次QR迭代中产生的 都是不可约的，否则，可以将问题分解为较小型的问题。这样，带原点位移的带原点位移的QR算法算法可以描述为：。（正交相似变换），分解），（形（初始化），的为取,2, 111kQRAQRRQIsAHessenbergAAkkkkkkk变换的变换称为原点位移的到这里QRA

51、Akk1,即每个所以，由于kkTkkkkkkTkkTkkkkkTkkkkQAQAQIsRQQQQsQRQQIsQR1,)( 特征值。的其余近似算法，就可逐步求出应用方法，即对近似特征值。采用收缩的作为且基本收敛于三角阵，并则在一定的条件下，序主子矩阵，若选取阶顺的是。设阵的序列相似于上述变换就产生一正交反复应用不同的位移相似。实际计算时，用相似，从而与原矩阵都与AQRRBAaAasnABaAAHessenbergsssAAAnnkknnkknnkkknnkijkkkk11211,1,1,11,1,11,1,knnknnknnknnknnknnaaaaAaa则与相邻元素进行，取准或者将充分小的准

52、则可以是判别 近似特征值。的一个就作为足上述准则时，可认为为给定的精度要求。满其中Aaaknnknn, 01,根据QR算法的收敛性质，位移量有下列两种取法： 最接近的一个。的特征值中与取为矩阵。knnknnknnknnknnkknnkaaaaasas1, 11, 121变换可表达为进行原点位移的矩阵转变换对具体计算时，用平面旋QRAHessenberg1。IsPPPRARIsAPPPTnnTTnn1, 12312121111223, 1,)(矩阵。仍为上容易验证，HessenbergA2例7.7 用带原点位移的QR算法求下列矩阵的特征值：611142121A解 先用镜面反射变换把A化为上Hes

53、senberg矩阵。按（7.4.3）式有4 .62 .002 .06 .35051,525/105/15/2000121AHHAwwIHTT,114997409. 0,392590761. 0, 5 . 6211s角元素，则有量，即取位移量为右下若按第一种方法取位移,cossin0sincos0001,1000cossin0sincos2222211111PP同理可得421062856. 6002432908. 00002432908. 0779171336. 4666846149. 00666846149. 0200234192. 04 . 6,021202899. 0183563711.

54、0743000879. 1076516671. 0137183510. 3844655679. 5)4 . 6(21121121IPPRAIAPPRTT421066615.6000866668465.4036401350.00036401350.0287735078.0,421066615.6000000006.00000000006.0862139274.4157002612.00157002612.0283205888.043AA。，矩阵的特征值，得。求左上角故有特征值866925525. 4287992138. 022421066615. 6213类似于上面的计算可得，矩阵的特征值，则有

55、量，即取右下角若按第二种方法取位移469693846. 6221s。421066615. 6000866649445. 4037723983. 00037723983. 0287716058. 0,421066615. 6000011074. 00000011074. 0861785493. 4162696110. 00162696110. 0282852106. 0,421048287. 6005374767. 00005374767. 0773105873. 4689146437. 00689146437. 0194154158. 0432AAA。，由此可得特征值421066615. 6866925525. 4287992139. 0321 该问题如果不用带原点位移的QR算法，而是用基本QR算法，则收敛速度很慢，计算结果为。28799