矩阵特征值和特征向量计算

1、 第一节第一节 引言引言第二节第二节 幂法及反幂法幂法及反幂法 幂法幂法 加速方法加速方法 反幂法反幂法第三节第三节 豪斯霍尔德方法豪斯霍尔德方法 正交相似变换正交相似变换(1)(1) 正交相似变换正交相似变换(2)(2)第四节第四节 QR 方法方法主主 要要 内内 容容以下是一些准备知识以下是一些准备知识 工程实践中有多种振动问题，如桥梁或工程实践中有多种振动问题，如桥梁或建筑物的振动，机械机件、飞机机翼的振建筑物的振动，机械机件、飞机机翼的振动，及一些稳定性分析和相关分析可转化动，及一些稳定性分析和相关分析可转化为求矩阵特征值与特征向量的问题。为求矩阵特征值与特征向量的问题。 3. ()0

2、 AIA xAxxxA设 为 的特征值，求相应的齐次方程的非零解（即求的非零解），称为矩阵 对应于 的特征向量。1. (), ()det()0 ()ijnnAaIAA 已 知求 代 数 方 程的 根 。称 为的 特 征 多 项 式 。2.()nA 特 征 多 项 式一 般 有个 零 点 ， 称 为的 特 征 值 。 但高次多项式求根精度低 , 一般不作为求解方法. 目前的方法是针对矩阵不同的特点给出不同的有效方法.其中的特征值及特征向量，求：例A 1210131012A32()d e t () 71 48 0AIA解： 的特征多式421 321的特征值为：求得A机器机器 求解求解1231111

3、 ;0;2111Axxx的 特 征 向 量 ：定理定理1 设为则为为数nnAR的特征值,且向Ax =x，其中x0，(1)ccA的特征值(c常c0);(2) 为-p A -pI的特征值， 即(A -pI)x =(-p)x;为；kk(3)A 的特征值设 为阵为-1-1(4)A非奇异，那么011 且A 的特征值，即A x =x 。设为 阶阵则iijnn (i = 1,2,Ln)n矩A =(a ) 的特定 理2征值， 11(1 )();d e tnniiiiiatrAA12n( 2 )() = .()()nnTARAA设则 定 理 31 11 212 22nnm mAAAAAAA设为块阵定 理 4 A

4、分上 三 角， 即1()()niiiAA个 对块为阵则其 中 每角均方，AByBP yA( 1 )与有 相 同 的 特 征 值 ；( 2 ) 如 果是的 特 征 向 量 ， 则是的 特 征 向 量 。1ABPBPAP设与为相似矩阵即存在非奇异阵 使则定理 (),5注注重要结论重要结论A矩 阵的 特 征 值 不 变 。1ABPAPA一 个 矩 阵经 过 相 似 变 换 ,亏损矩阵亏损矩阵nnARkk设阵个数 为对 应线无 关数矩有 一重的特 征 值 ,但于 它 的性的 特 征 向 量目一个亏损矩阵是一个没有足够特征向量的矩阵，一个亏损矩阵是一个没有足够特征向量的矩阵，亏损矩阵在理论与计算上存在巨

5、大的困难。亏损矩阵在理论与计算上存在巨大的困难。n nAR(1)矩阵可对角化， 即存在非奇异定6 理矩阵使12 nP APAn-1的充分必要条件 是 具有 个线性无关的特征向量。12 , n nmARmmnxxx12n(2)如果有个()不同的特征值,则对应的特征向量,线性无关。7n nAR为对称正定定理 ：矩阵，则：(2)AAnP(1) 的特征值均为实数；有 个线性无关的特征向量；(3)存在一个正交矩阵使得12TnP AP1,2,iinA()为 的特征值12(,)njPuuuu的 列 向 量jA为 的对应于 的特征向量。n nARQ设矩阵则存在正交矩阵使11121222nnTnnRRRRRQ

6、AQR(实Schur定理10分解)：(1,)iiR im 其中为一阶或二阶方阵，iiRA每个一阶是 的实特征值，iiRA每个二阶 的两个特征值是 的两个共轭复特征值定义1：(, )( )( , )Ax xR xx x,0n nARx x设矩阵为实矩阵，xRayleigh称为对应于 的商n nAR：设为对称矩阵，特征值定理8次序记为2n 1(, ) ();( , )nnAx xxRx x 11.10(,)2.m ax;(,)nxRxA x xx xn0(,)3.m in;( ,)nxRxAx xx x则(),ijn nAa定义2：设令：对于矩阵特征值界如何估计？对于矩阵特征值界如何估计？1(1

7、) | (1, 2 ,) ;nii jjjirain(2)| |,iiiiDzzar zC集合iiiar称复平面上以 为圆心， 为半径的所有圆盘.AGerschgorin的圆 盘称为称为Gerschgorin圆盘定理：(),ijnnAa(1)设则A的每一个特征值都必属于下述某一个圆盘之中:1|niiiijjjiara或者A的特征值都在复平面上n个圆盘的并集并集之中之中. .(2)如果A的m个圆盘组成一个连通的并集S，且S与余下的 则S内包含A的m个特征值个圆盘是分离分离的，n-mn-m特别：若A的一个圆盘iD与其它圆盘是分离的即孤立的 中精确地包含A的一个特征值.,则iDA设 为 的证明：特征

8、值，即12, =()0.TnAxxx xx其中,1 |max | |0,kiinxxx 记Axxk考察的第 方程1=nkjjkjaxx1 (-)nkkkkjjjjkaxax或1|nkkkjkjjkaar即11| |nnkkkkjjkkjjjj kj kaxaxxa于是即结论得证.4 101 01.1 114EAxample估计矩阵的特征值的范围A解 ：有 三 个 圆 盘11111: |4 |1njjjDra33313: |4 |2njjjDra22212: |2njjjDra由上述定理结论可知A的三个特征值位于三个圆盘的并集中，所以D1内恰包含A的一个实特征值1352323.ADD的其它两个特

9、征值，包含在中由于D1是孤立的所以，2323DD问题：如何进一步估计上面两个特征值分别在什么范围？解决途径：若能够改变圆盘的半径，则 有可能将圆盘进行分离，从而可进一步分析特征值的范围.事实上，利用相似矩阵的性质，可使事实上，利用相似矩阵的性质，可使A A的的某些圆盘半径及连通性发生变化某些圆盘半径及连通性发生变化. .具体实施具体实施?11121(1,2, ) ininD1适当选取作=-1ijjin naAD AD对 作相似变换,.即 可 解 决 上 述 提 出 的 问 题()ijjijn nin naAaB即,AB而与有 相 同 的 特 征 值1110.9D若 选 取作相似变换114101

10、010/90.90.94AAD AD对上边同一例题对上边同一例题1A可同上求得 的三个圆盘为：1: |4|1E2:|1 9 / 9E3:|4 |1 .8E123EEE 显 然135 ;219 / 919 / 9;35 .82 .2 幂法幂法是一种计算矩阵的按模最大的特征值与相应的特征向量的迭代方法。适合于大型稀疏矩阵反幂法反幂法是计算Henssenberg阵或对角阵的对应一个给定近似特征值的特征向量的有效方法.一、幂法12(1, 2 ,)innnAin设阶 实 矩 阵的 特 征 值满 足12(1,2, ),ininx xx与相应的特征向量线性无关。且有完全的特征向量组即即构造一个向量序列幂法的

11、基本思想：任意取一个非零的初始向量( 0 )v由矩阵A( 1 )( 0 )( 2 )( 1 )2( 0 )(1 )()1( 0 ),kkkvA vvA vAvvA vAv称为迭代向量。由此计算特征值和特征向量。(0)1 1(1,2, )(1,2, ),iiiniiivinxninvx且为简便，不妨设。因为 线性无关，故必存在 个不全为零的数使得。(1)( )1 (0) kkkvAvA v由11()nkiiiAx11nkiiiix(1)11211122111 () ()kkkknnnvxxx( )11 1()kkkvx21knikiiix其中110,(2,3, )iin设由121lim()0nk

12、iiikix11 lim()0 kikiikx得k故只要 充分大，(1)11111111kkkkvxx必有必有(1)kiv可把作为与相应的特征向量的近似。( )11 1()kkkvx由 lim0 kk及()111limkkkvx ()11kkvAk 即表明序列越来越接近 的对应于的特征向量，也即当时()111kkvx 主特征值主特征值1?()如 何 算()() kkivvi为的 第 个 分 量 。(1)111 1 ,kkvx由( )11 1 kkvx 及(1)1( ) (1,2,) kikivinv即相邻两个迭代向量分量的比值收敛到即相邻两个迭代向量分量的比值收敛到主特征值主特征值. .21

13、幂 法 的 收 敛 速 度 依 赖 于 比 值，A按 上 面 式 子 计 算 矩 阵按 模 最 大的 特 征 值 与 相 应 的 特 征 向 量 的方 法 称 为 幂 法 。11 lim()0 kikiikx由比值越小，收敛越快。两点说明：(0)110,v） 如 果的 选 取 恰 恰 使 得幂 法 计 算 仍 能 进 行 。1x它 在方 向 上 的 分 量 不 为 零 ， 于 是以 后 的 计 算 就 满 足 所 设 条 件 。( ),kv因为计算过程中舍入误差的影响，迭代若干次后，必然会产生一个向量因此，幂法实际使用的计算公式是( )11111 2,(1)0(1)kkvx ）因计算过程中可能

14、会出现溢出或成为的情形。解决方法：每次迭代所求的向量都要归一化。( )( )( )( )( )1( 1)( )( 1)1 m ax(0,1,2, ) kkkkkrrii nkkkruvvvvkvAuv 注意注意( )( )1 max.kkrii nrvv 是满足的最小的下标(1)(0)(1)(1)(0)max()max()vAvuvAv(2)2(0)(2)(2)2(0)max()max()vA vuvA v (0 )()(0 )m ax()kkkA vuA v(1)(0)(0)vAuAv2(0)(2)(1)(0)max()A vvAuAv(0)( )(1)1(0)max()kkkkA vvAu

15、Av( )(1)(0) kkkvAvA v由1nkiiiix()11121 ()nkkkiiiivxx( 0 )()( 0 )m a x ()kkkAvuAv1112111121() max ()nkkiiiinkkiiiixxxx()1111 max (kkkxux11()m a x ()xkx同上讨论有1112()11111121() max()nkkiiiiknkkiiiixxvxx1112()111121() max()max()nkiiiiknkiiiixxvxx1k 13n nARn：设有 个线性无关定理的特征向量，112| |,n主特征值 满足(0)(0)vu则对任意的非零初始向

16、量( )( ):kkvu按下述方法构造的向量序列及( 0)( 0)()(1)()()()0ma x()kkkkkkkvuvAuvvu()111( 1 ) l i mma x ()( 2 ) l i mkkxux则有11 1.(),(,),2 .1,03.max4. 5., ,66.,1,3ijnriri nrAavvvNkrvvvuvAuvvkNkk ：输入初始向量误差限 ， 最大迭代次数。置求整数 ，算法使，v计算置若输出停机；否则，转若置转 ；否则， 输出失败信息，停机。210021012A例：用幂法求矩阵的按模最大的特征值和相应的特征向量。( 0 )3( 0 , 0 , 1 ),1 0.

17、Tv取( 0 )( 0 ) (0 , 0 ,1),Tuv解 ：(1)(0)(0, 1,2) ,2,TvAu( 1 )( 1 )( 0 ,0 . 5 , 1 ),Tvu(2 )(1) (0.5,2, 2.5) ,2.5, TvAu (8)(7)(2.7650948, 2.9981848,2.9990924)vAu2.9990924(8)(9)(8)(0.9219772, 0.9996973,1)(2.8436517, 2.9993946,2.9996973)uvAu3 2.99969732.99909240.000604910由12 .9 9 9 6 9 7 3 . 故 (2.8436517,

18、2.9993946,2.9996973)x 。相应特征向量为1233,2,1A事实上， 的特征值，11 -1,1 T与对 应 的 特 征 向 量 为 （ ，）。212 .3此 例 中 比 值 为两种特殊情况两种特殊情况 .12 121mnm前面假定如果按模最大的特征值有多个，即幂法是否有效？1,112mAmn（ ） 是 重根，即矩阵仍有 个线性无关的特征向量。(1)111 1111 ()()1111kkvxxm mkkmnxxn nmm此时有 ,1(1)1 ()11 1kmkkvxxm m显然，只要不全为零，当 充分大时，就有1 11(1) 12( )(1)xxAm mkvimkvikv因也是

19、矩阵 相应于的特征向量，故有为相应的特征向量，即对这种情况幂法仍然有效。12132,An （）且 矩 阵有个 线 性 无 关 的 特 征 向 量 。(1)111311 1223 3111( 1)() ()kkkkknnnvxxxx (1)(21)2111122(2 )211122 () ()kkkkkvkvxxvxx由上式可知，是个摆动序列，当 充分大时，有(2)2(2)( )11,2( )/kkkiiikivvvv(1)1111 122( )11 122 ( 1) ( 1)kkkkkkvxxvxx 又由(1)( )1111 1(1)( )11112 22 2( 1)kkkkkkkvvxvvx

20、 故在这种情况下，仍可按幂法产生向量序列。121211 nmmnA综上可知，当 的特征值分布为或时，用幂法可以计算出及相应的特征向量。An此外，当矩阵 无 个线性无关的特征量时，幂法收敛很慢，亦应考虑改用其他方法。幂法计算简便易行，它是求大型稀疏矩阵按模最大特征值的常用方法。(1)( )( )12kkkvAvv 如果按迭代所得向量序列呈有规律的摆动，则可能为的情况。否则应考虑用别的方法求解。 因为幂法的收敛速度是线性的，而且依赖因为幂法的收敛速度是线性的，而且依赖于比值于比值 ，当比值接近于，当比值接近于1 1时，幂法收敛时，幂法收敛很慢幂法加速有多种，以下介绍两种。很慢幂法加速有多种，以下介

21、绍两种。210 AAI（ 一 ） 原 点 移 位 法矩 阵与的 特 征 值 有 以 下 关 系 ：00iiAAI若 是 的特征值，则就是的特征值，而且相应的特征向量不变。(1)120101122101010()() ()kkknnnuuu010002101 (2,3,)iiin适 当 地 选 取， 使 得且(1)( )0(1)( )0 ()kkkkAIxAxxAI x如果对矩阵按计算，则有010A- I -A这样，用幂法计算的最大模特征值及相应特征向量的收敛速度比对 用幂法计算要快。这种加速收敛的方法称为原点移位法。0012312300)nnA原点移位法使用简便，但选取困难。在一些简单情形，可

22、估。如当矩阵 的特征值满足（或时，20210122221121112 nnnnnn且 因此，用原点移位法求 可使收敛速度加快。02020101(),2 (2,3, )niin取(4)4(5)5454 (3.1000568,2.214326, 0.9687661) 3.1000568(3.0999984,2.2142846, 0.9687501) 3.09999840.000058410 xx 0414051302.9,102.8AA -4例：,用原点移位法求矩阵 的按模最大的特征值，要求误差不超过10 。(0)(1)( )0(1,1,1) ,()TkkxxAI x解：取按进行计算06.9140

23、510.10100.1AI1 3.09999842.95.9999984A所以，矩阵 的按模最大的特征值为 (3.0999984,2.2142846, 0.9687501) .Tx 相应的特征向量为 1236,3,2.8,A不难求出， 的特征值为 1211222112121 lim0 ()22kkkkkkkkkkkkkkkkkkkkkkkkaaaacaaaaaakaaaaaaaaaaaaaaaaaaaaaAitkena如果序列线性收敛到 ，即则当 充分大时，有序列比更快地收敛到 ，这就是加速法。将这一方法用于幂法所产生的序列，可加快幂法的收敛速度。（二）、幂法的埃特肯(Aitken)加速101

24、012210021001021 1.(),( ,),2 .1,0,0,13.max4. ()5. 26., , ,77., ,ijnriri nrAaxxxNkrxxxxyxAyxaaaaaaxkN 输入初始向量误差限 ，最大迭代次数 。置求整数，使，计算置计算若输出停机；否则，转若置算法：0,1,3kk 转 ； 否则，输出失败信息，停机。(也可采用幂法迭代两步或三步，加速一次的方法） 反幂法是计算矩阵按模最小的特征值及特征向反幂法是计算矩阵按模最小的特征值及特征向 量的方法，也是修正特征值、求相应特征向量的量的方法，也是修正特征值、求相应特征向量的 最有效的方法。最有效的方法。11 , 1

25、AnnuAAuuuA uA uu设 为阶非奇异矩阵，为 的特征值 与相应的特征向量，即111 AAAAA 此式表明，的特征值是 的特征值的倒数，而相应的特征向量不变。因此，若对矩阵用幂法，即可计算出的按模最大的特征值，其倒数恰为 的按模最小的特征值。这就是反幂法的基本思想。1(1)( ) kkAAAxx因为的计算比较麻烦，而且往往不能保持矩阵的一些好性质（如稀疏性），因此，反幂法在实际运算时以求解方程组 1.AALU反幂法计算的主要步骤：对 进行三角分解(1)1( )(1) kkkxA xxA代替幂法迭代求得，每迭代一次要解一个线性方程组。由于矩阵在迭代过程中不变，故可对 先进行三角分解，每次

26、迭代只要解两个三角形方程组。( )( )1( )( )( )2.max, kkrii nkrkkrxxxxy 求整数 ，使得 计算 ( )(1)3. kkLzyUxz解方程组* 0 ()iiAAI用带原点移位的反幂法来修正特征值， 并求相应的特征向量是非常有效的。设已知 的一个特征值 的近似值为， 因接近 ，一般有 故是矩阵的按模最小的特征值，反幂法的一个应用* ()iiAI且由上式可知，比值较小。因此，对用反幂法求一般收敛很快，通常只要经过二、三次迭代就能达到较高的精度。*11.(),( ,),2 . 1,1ijnAaxxxNku输入 近似值 ，初始向量 误差限 ， 最大迭代次数 。置算法：

27、*1*3.()max5. 1116., ,77.,1,4riri nrAILUrxxxxyLzyUxzxxukNkku 作三角分解 4.求整数 ，使，计算置若则置,输出停机； 否则，转若置转 ； 否则，输出失败信息，停机。(0) (0,0,1) .210021012TAxA：用反幂法求矩阵 接近2.93的特征值， 并求相应的特征向量，取 这里矩阵 A 为 例 ，2.930.9310 2.9300.931010.93AIAI：对作三角分解得解41000.9310 01000.93101/0.931000.931/0.933 3.0000954,3 10(1, 0.9992431,0.999147

28、8)(1,-1,1)0.001TTur 按算法迭代 次，与准确值 的误差小于 ，与准确值比较， 残差2 Jacobi 方法12 Jacobi,:(1), (,) (1,2, ),2(),Tniijn nTAQQ AQdiaginAQAaQBQ AQ 方法是求实对称矩阵的全部特征值及相应特征向量的一种方法 它基于以下两个结论任意实对称矩阵 可通过正交相似变换化成对角阵 即存在正交矩阵使得其中是 的特征值中各列即为相应的特征向量。（ ）在正交相似变换下，矩阵元素的平方和不变。设为正交矩阵，记22,1,1(), ,nnijn nijiji ji jbabJacobiA则方法的基本思想是通过一次正交变

29、换 将 中一对非零的非对角元素化成零 并且使得非对角元素的平方和减少。反复进行上述过程，使变换后的矩阵的非对角元素的平方和趋于零，从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。一、矩阵的旋转变换11cossin1( )1sincos11ijAnV设 为 阶实对称矩阵，考虑矩阵(1)(1)(1)22(1)22(1)(1)(1 ( ) ()cossinsin2 sincossin2 cossin (, ) TijijijijiiiijjijjjiijjijikkiikjkjkVAV AVaaaaaaaaaaaaaki ja是正交矩阵，若记则有)(1)(1)(1)(1)(1)sincos (

30、 , )1()sin2cos2 20,22/() kjikjkkllkklijjijjiiijijijiijjaaaaaak li jaaaaaatgaaa 如果取 使得(1)(1)(/4)0,.ijjiijjiaaAaa则有得到一个使 中非零的非对角元素变成零的正交相似变换(1)(2)( )(1)2(1)(1) 2(1)(1)(1)(1) , ( ),()() ( , )co ssinkklklk lk lklklikkiikjkjkkAAAAAE AaE Aaaak li jaaaaaa对重复上述过程得到一个矩阵序列。可证，虽然这种变换不一定能使矩阵中非对角元中零元素的个数单调增加，但可保

31、证非对角元的平方和递减。以 与为例：设由(1)(1)(1) 2(1) 2(1) 2,2222,sincos (, ) ()()2()() 2()( )2( ) jikjkklikjkk l i jk i jklikjkijk l i jk i jaaki jE AaaaaaaE AaE A 上式表明，在上述旋转变换下，非对角元的平方和严格单调递减，因而对角元的平方和单调递增。正利用这点，导出了Jacobi方法。(1)( )( )( )( )1,( ,0,2., ,max3. 2kkijkkkijlml m n l miiAAAJacobiAaJacobikAAi jaaaactg通过一系列旋转

32、相似变换将 变成求得 的全部特征值与特征向量的方法称为方法。如果在对 作相似变换的过程中，每一步都选绝对值最大的非对角元素以此确定旋转矩阵，这种方法称为古典方法。其计算过程如下：1.令求整数使得计算旋转矩阵)( )2( )( )2,( )(1),21 cos,sin(1 ,)kkjjkijkijabtgsign aaaacdbcVVb (1)(1)2( )2( )( )(1)2( )2( )( )(1)(1)( )( )(1)(1)( )( )(1)(1)4.2, 2, , (1,2, , )kkkkkkkkkiiiijjijjjiijjijkkkkkkkkilliiljljlljiljlkk

33、lmmlAac ad acdaad ac acdaaacadaaadacaln li jaaa 计算( )(1)(1)(1)(1)2 ( ,1,2, , ) 0 5.()() klmkkijjikklml mn l mi jaaE Aa置计算(1)(1)(1)(1)(0)(1)1122( ) 6.(), 12 l mkkkknnnE AaaaQVVVkkA 若则为特征值，的各列为相应特征向量；否则，返回 ，重复上述过程。一般地，古典Jacobi法不能在有限步内将 化成对角阵，但有以下收敛性结果。( )(0)( )( )( )( )2( )(1)(1)2(1)2, lim()0, max1 ()

34、()(1) () ()( kkkkkijlml mkkijkkkiijjAnAJacobiAAAE AJacobiJacobiaaaE An nAaaa设 为 阶实对称阵，对 用古典法得到序列其中则即古典法收敛。由古典法计算过程故有又由计算的公式可得定理：证：( )2( )2( )2(1)2(1)2( )2( )2(1)( )( )2(1)( )1) ()2()() ()() () ()()2() ,22 ()1 ()1( )(1)(1)2 11,lim(1)kkkiijjijkkkkiljliljlkkkijkkkkaaaaaaE AE AaE AE AE An nn nEn n于是有因(

35、)2()lim1( )0(1)kkkAE An n JacobiJacobiJacobi定理表明，古典法是收敛的，进一步分析还可得出法收敛较快。另外，这种方法对舍入误差有较强的稳定性，因而解的精度高，且所求得的特征向量正交性很好。它的不足之处是运算量大，且不能保持矩阵的特殊形状（如稀疏性），因此法是求中小型稠密实对称矩阵的全部特征值与特征向量的较好方法。(0)12210 121012 1,2,20,41cossin,21102211 ()02 20 01AAijctgVV用Jacobi方法求 的特征值。首先取因故有于是例：解：(1)(0 )(0 )(0 ) 111100222221011110

36、121022220120010011102103211222TAVAV(1)(2)(1)(1)(1) 1,3,22,cos0.88807,sin0.459700.8880700.45970 0100.4597000.888070.633970.3250600.3250630.6279600.627962.36603 2,3,TijtgVAVA Vij 再取下面应取重复上述(5)(5)22.,0.585790.00203830 0.00203833.414010.01675800.0167582.00020()2(0.00203830.016758 )0.00056997AE A过程 如此继续下

37、去 可得123123 3.41401,2.00020,0.58579223.4142136,2,220.58578640.0002036 :,0(),( ),kkkijkijijjiAkaVaa 所以 的特征值为准确值最大误差为古典Jacobi法的改进 取定正序列以为限 逐行检查非对角元 若就跳过 否则以消去元和反复进行上述过程 直到+1,1kkkkA所有非对角元的绝对值均小于再以为限 进行第轮循环消元。当充分小时，所得到的矩阵的对角元即为 的全部特征值。一、基本一、基本QRQR方法方法 60年代出现的QR算法是目前计算中小型矩阵的全部特征值与特征向量的最有效方法。 实矩阵、非奇异。 理论依据

38、：任一非奇异实矩阵都可分解成一个正交矩阵Q和一个上三角矩阵R的乘积，而且当R的对角元符号取定时，分解是唯一的。()(1)(1)(1)11111(2)1(2)111 QRQR (1,2,). ,kkkkkkAQ RkAR QAAAAAQ RQARAR QQAQAA基本方法的基本思想是利用矩阵的分解 通过迭代格式将化成相似的上三角阵（或分块上三角阵），从而求出矩阵 的全部特征值与特征向量。由 即 。于是 即与相似。() (2,3,)kAAk 同理可得，。故它们有相同的特征值。定理（QR方法的收敛性）,n nRij设 A= a(1)0;12n如果 A的特征值满足: 11211(2)(,),nAAXD

39、XDdiagXXLU 有标准型 其中 且设 有三角分解 （L 为单位下三角矩阵，U 为上三角矩阵）， 12*()knARk k本质上则由 QR算法产生的 A本质上收敛于上三角矩阵，即当 时( )( )( ),(1)lim;(2),lim0;(2),;kkijkiiikkijkkijAaaijaija若记 则当时当时的极限不一定存在 可证，在一定条件下，基本QR方法产生的矩阵序列A(k) “基本”收敛于一个上三角阵（或分块上三角阵）。即主对角线（或主对角线子块）及其以下元素均收敛，主对角线（或主对角线子块）以上元素可以不收敛。特别的，如果A是实对称阵，则A(k) “基本”收敛于对角矩阵。 因为上

40、三角阵的主对角元（或分块上三角阵中，主对角线子块的特征值）即为该矩阵的特征值，故当k充分大时， A(k)的主对角元（或主对角线子块的特征值）就可以作为A的特征值的近似。基本的QR方法的主要运算是对矩阵QR分解，分解的方法有多种。介绍一种Schmit正交化方法。12122111122211212121212111211211111112222121,(,) (1,2, )., /. , ,0, /,1, nTjjjnjAnAa aaaaaajna abaabaa b baaaa aa aa aa aaab ba aaa aabbbbbbbb设 为 阶非奇异实矩阵，记 其中取取则2,0.b1112

41、1111 , /(2,3, ),1(1,2,),kkkkiiikkknkkkkkkkkbaa b bbbbknb bbbknaa b ba bbb b一般地取 则向量组正交，且即1111121212112211, ,QRkkkkkkknnnnnnnnaabbabbbbAaaabbbaababbabQRbabbSchmit即 于是这就是用正交化方法对矩阵进行的分解。 基本QR方法每次迭代都需作一次QR分解与矩阵乘法，计算量大，而且收敛速度慢。因此实际使用的QR方法是先用一系列矩阵相似变换将 A 约化成拟上三角矩阵(称为上Hessenberg 矩阵)，然后对此矩阵用基本QR方法。因为拟上三角矩阵具

42、有较多零元素，故可减少运算量。化A为相似的拟上三角阵的方法有多种。123111222211 210121012.: ( 2 ,1, 0 ),(1, 2 ,1),( 0 ,1, 2 ).21 (, 0 ),55,42136(1, 2 ,1)(, 0 )(,1),5555: 5TTTTTTTS c h m itAQ Raaaababaabb用正 交 化 方 法 对 进 行分 解解因 而 有例2222333113223332365(,),707070,246(,)777 123 (,)141 1 44TTTbbbbaabbabbbbb121212112211 , , ,253701145451515

43、670214070516700570314002 40 7nnnnnnnnAa aab bbaa ba bba bba bb 所以1HouseholderAHHouseholderA（ ）用变换将矩阵 化成拟上三角矩阵或上海森柏格矩阵。(2) 用变换将对称矩阵 化成对称三对角矩阵。为了求矩阵特征值先进行初等变换把矩阵变成较简单形式8.3、豪斯豪尔德(Householder)方法221122112122122212 (,)1,12222122 22212TnnnTnnnnwwwwwwww ww ww www wHIwww ww wwHouseholder设向量 满足 则称为 矩阵或反射矩阵。1、

44、豪斯豪尔德(Householder)变换11;THHHH 可证其具有以下性质：（ ）是实对称的正交矩阵，即 (2) det()1;H (3)11n1,1,;Hw仅有两个不等的特征值，其中 是重特征值是 单重特征值为其相应的特征向量:0,().TwoS w xRH xwxw (4) 考虑以 为法向量过原点 的超平面为任意 的数 有( )22,1,nx yRyHouseholderHHxxy 定理 : 设为中任意非零向量 且则存在矩阵使得。222222222222222222() 2, (2)2( ) .2 () () 22() .TTTTTTTTTxxywHIwwxxyxxyHxIwwxxxxy

45、xxxyxxyxxyxxyx xxy xxxxyxHxxy 证：令于是 由 范数的定义. 代入上式得222,()()iiiiixHouseholderxyexHouseholderxxe sign xwxxe sign x此定理表明，对任一非零向量 都可以构造一个变换，它将 变成事先给定的单位向量的倍数。特别地取 则 经过 变换后可变成只有一个分量不为零。实际计算时，为避免误差取 。111211121222123233311(2,3, ),nnnnnnnnniihhhhhhhhHhhhhhhin称形如 的矩阵为拟上三角阵，也称为上海森堡（Hessenberg)阵。如果此对角线元全不为零 则称该

46、矩阵为不可约的上Hessenberg矩阵。 Householder A讨论用 变换将一般矩阵 相似变换 成 Hessenberg阵.11111111,1000 01HouseholderHHH AHHHHHnHouseholder首先，选取 矩阵 使得经 相似变换后的矩阵 的第一列中有尽可能多的零元素。为此，应取 为如下形式其中为阶 矩阵。11121112131111122122221213122211111(,) ,(,) ,.(1,0,0) ,2TTnnTnnnnTaaHH AHaaaaH aH AHaaaaaaAaaHH aH AHn 于 是 有 其 中由 上 节 定 理 ， 只 要 取

47、使 得 就 会使 得 变 换 后 的 矩 阵 的 第 一 列 出 现 个 零 元 。111111211111122221122,() ()1000*0100*00* *00waae sign awHHaae sign aHouseholderHH H AH HH为避免在计算 时会产生较大的误差 取 。同理，可构造如下列形式 矩阵 使得*12222112222, .nnnnnHouseholderH HHHH H AH HHHHHessenbergAH如此进行 次，可以构造 个 矩阵 使得其中 为上 矩阵。特别地，当 为实对称矩阵，则经过上述正交变换后，变为三对角阵。:5222321052 22

48、 2 02100241HouseholderAA例 用 变 换 将 矩 阵 化 成 拟 上 三 角 阵 。1111122(1,0,0) ,(1,00) .21,02TTaH aHIHIHouseholderHH 解：因 由 为使 矩阵 满足 222(), () (2,2)2(1,0)(22,2) 222 ,22 22iiiiTTTTxxe sign xwxxe sign xwwww 由 公式 ，。22 222222104422 2012222244221000010022 0022220022THIwwH2112 100052223 201001052 22 222 0002102222024

49、10022100052510100103222 000223220012220022HH H AH H 于是有3、拟上三角矩阵的QR分解2121(2)(2)(2)11112131(2)(2)(112223221 1 0(cossin00sincos000011nnHnGivensHQRhVrhhhhhhV H因为拟上三角阵 的特殊形状，通常用个旋转变换（又称变换）可将它化成上三角矩阵，从而得到 的分解式。具体步骤为：设否则进行下一步），取旋转矩阵 则2)(2)(2)(2)32333(2)(2)1(2)221121111112111 cos, sin, .nnnnnhhhhhhhHrhhrr其中

50、232222232(3 )(3 )(3 )(3 )11213111(3 )(3 )(3 )223212(3 )(3 )(3 )23331332(3 )(3 )(43414 0(10cossinsincos 11 nnnnnnnnhVrhhhhrhhhhhhVHhhh（）（）设否 则 进 行 下 一 步 ） ， 再 取 旋 转 矩 阵 则(3 )3 )(3 )(3 )1( 2 )( 2 )( 2 )2( 2 )23222222223222 cos, sin, ()() .nnnnHhhhhrhhrr其 中( )(1)1( )( )( )( )1111111( )( )( )11111( )( )

51、( )1( )( )( )1111( )( )1 1 kkkkkkkkkknnkkkkkkknknkkkkkknknkkkkkknknkknnnnkHVHrhhhhrhhhhhhhhhhh假设上述过程已进行了步，有()11()()1()2()21 0,11 cossinsincos1 cos, sin, ()() . kkkkkkkkkkkkkkkkkkkkkkkkkkhVhhrrrhh设取其中(1)(1)(1)11111(1)(1)1( )(1)(1)(1)(1)111111(1)(1)(1)21212(1)(1)1 1kkkkknkkkkkknkkkkkkkkkknknkkkkkknknk

52、knnnnrhhhrhhVHHhhhhhhhhn于是因此，最多做次旋转变( )( )( )112131( )( )2232( )( )1122133 nnnnnnnnnnnnnnnrhhhrhhHVVV HRrhr换，即得121 32121 32123 (2,3, ) 4 ,( ) iiTTTnnTTTnnVinHV VVRQRQ V VVnQRO nHRQQRQR因为均为正交矩阵，故其中仍为正交矩阵。可算出完成这一过程的运算量约为比一般矩阵的分解的运算量少一个数量级。可证明仍是拟上三角阵，于是可按上述步骤一直迭代下去，这样得到的方法的运算量比基本方法大为减少。需要说明QR的是，通常用方法计算特征值，然后用反幂法求其相应的特征向量。222532644445 (6,4)64 (1,0)(652,4) (0.957092,0.289784) 2100.9160250.27 73500.8320500.55 2010.2773500.0839747TTTTTQRAAwwwwIww用方法求矩阵的全部特征值。首先将 化成拟上三角阵解：，取例：147000.5547000.83