Cht8矩阵特征值问题计算

1、18.1 8.1 引引 言言第第8 8章章 矩阵特征值问题计算矩阵特征值问题计算物理、力学和工程技术中很多问题在数学上都归结为求矩阵的特征值问题。例如，振动问题(大型桥梁或建筑物的振动、机械的振动、电磁震荡等)，物理学中的某些临界值的确定。它们都归结为下述数学问题。nnnnnnnnijaaaaaaaaaa212222111211)det()( ,)( AIA则称已知定义1定义12.|) 1()( 12211特征多项式特征多项式的为AAaaannnnn.)( .(1.1) 0)det()( 的所有特征值的集合表示的的根称为的特征方程AAAAIA特征值特征值.(1.2) 0)( ,特征值特征向量特

2、征向量的的对应于称为的非零解相应的齐次方程组的为设AxxAIA3.210131012 的特征值及其特征向量求A例1例1.1 ,10 )4( , )3()( )( , )2()0( (1) , 11xxAAAxxAAxxIAIAAxA即的特征值为且非奇异，则设；即的特征值是；即的特征值为；常数的特征值是向量，则是对应的非零特征的特征值是矩阵设kkkknnppppcccR定理1定理14.)det( )2(,1)(tr1 (1) ,), 1( 1niiiinianiniAAA则的特征值是矩阵若定理2定理2).()( , 3AAATnnR则设定理定理5. )()( , , 122211211miiii

3、immmmAAAAAAAAAAA则均为方阵其中每个对角块为分块上三角阵设定理4定理4. , )2(; ) 1 (, , 1的特征向量是则的特征向量是若有相同的特征值与则使非奇异即为相似矩阵与若APyByBAPPPBABA定理5定理56. , 亏损矩阵亏损矩阵定义2定义2为，则称个数少于线性无关的特征向量的且其对应的重特征值有一个如果设AAAkkRnn., ,)( )2(. 1 2121211线性无关对应的特征向量则个不同的特征值有若个线性无关的特征向量具有的充要条件是使非奇异矩阵即可对角化，）（mmnnnnnxxxnmmRnRAAAPPPA定理6定理67则为对称矩阵设对称矩阵的正交约化, )7

4、(nnRA定理定理；个线性无关的特征向量有的特征值均为实数；nAA )2( (1).)(,), 2 , 1( )3(21211的特征向量为对应于向量的列而的特征值为且，使得存在正交矩阵jjnin,niuuuuPAAPPP8.), 1( ,| | )2( | 1,)( 圆盘为半径的为圆心以为复平面上以称，）（令设nGerschgoriraDnizrazzDnijaraiiiiiiiiijinnijCA定义3定义3. ), 1( |,| ,)( (1) )( nijniaaanGerschgoriijiinnij某个圆盘之中一个特征值必属于下列的每则设圆盘定理AA定理8定理8.S ,S, 个特征值

5、的中恰有则个圆盘分离与其余且个圆盘构成一个连通域个圆盘中有如果上述的mmnSmnA(2)(2)9.), 2 , 1(. ),(diag11结果质获得特征值的进一步改变，根据相似矩阵性和连通性有时可使某些圆盘半径适当选取，得到选取非奇异对角矩阵niainnijijnADDD.411101014 的特征值的范围估计A例2例2.49 . 09 . 001014,119101910ADDD2|4| 2|1|4|2 . 28 . 55339192919110.), 2 , 1(, )(922211211的特征值为其中使则存在酉矩阵，设定理AnirrrrrrrRSchuriinnnnTnnRAUUUA定理

6、定理.), 2 , 1(, )(1022211211的两个共轭复特征值的两个特征值是二阶时为的实特征值，当是为一阶时其中当使则存在正交矩阵，设分解实ARRARRRRRRRRAQQQAiiiiiiiimmmmTnnmiRSchur定理定理11.)()(),()( , , 商商瑞瑞定义4定义4RayleighRRnn利的为关于向量称对于任一非零向量阶实对称阵是设xxx,xAxxxA.)(),(min 3,)(),(max 2 ,)(),( 1 . , 11111xx,xAxxx,xAxxxx,xAxAAxxxx00nnRnRnnnRn）（）（对于任何非零向量）（则的特征值为阶实对称阵为设定理定理1

7、2 证明证明 只证 1. 由于 为实对称矩阵，可将 对应的特征向量 正交规范化，则有 n,21Anxxx,21.),(ijjixx 设 为 中任一向量，则有展开式 0 xnR,0,211221niiniiixxx于是 .),(),(1212niiniiixxxAx从而1成立. 结论1说明瑞利商必位于 和 之间. n1138.2 8.2 幂法及幂法及反幂反幂法法.1、幂法幂法幂法是一种求实矩阵A的按模最大的特征值1及其对应的特征向量x1的方法。特别适合于大型稀疏矩阵。 ., , ,)(2121nnnnnnijRaxxxA对应的特征向量为为其特征值有一个完全特征向量组设(2.1)

8、, 21n满足的主特征值是实根，且并设A.11的基本方法及现在讨论求x14 幂法的基本思想是任取一个非零的初始向量 ，由矩阵 构造一向量序列 0vA,011021201vAAvvvAAvvAvvkkk（2.2）称为迭代向量. 由假设， 可表示为 0v),0(122110设nnxxxv（2.3）于是 15),()/(1112111122211101kkniikiiknknnkkkkkxxxxxxvAAvv其中.)/(21niikiikx由假设),3 ,2(1/1nii故 ,0limkk（2.4）.lim111xvkkk从而 这说明序列 越来越接近 的对应于 的特征向量，或者说当 充分大时 kkv

9、1A1k16,111xvkk（2.5）即迭代向量 为 的特征向量的近似向量（除一个因子外）. kv1 再考虑主特征值 的计算，用 表示 的第 个分量，则 1ikv )(kvi,)()()()()()(1111111ikiikiikikxxvv（2.6）故 ,)()(lim11ikikkvv（2.7）也就是说两相邻迭代向量分量的比值收敛到主特征值. 17 这种由已知非零向量 及矩阵 的乘幂 构造向量序列 以计算 的主特征值 （利用（2.7）式）及相应特征向量（利用（2.5）式）的方法称为幂法幂法. 0vAkAkvA1 由（2.6）式知， 的收敛速度由比值 来确定， 越小收敛越快，但当 时收敛可能

10、就很慢. 11)()(ikikvv12rr112r 定理定理12 12 设 有 个线性无关的特征向量，主特征值 满足 nnA Rn1,321n则对任何非零初始向量 , (2.4), (2.7) 式成立. )0(1v18 如果 的主特征值为实的重根，即 ，且 Ar21,1nrr又设 有 个线性无关的特征向量， 对应的 个线性无关特征向量为 ，则由（2.2）式 An1rrxxx,21,)/(11110nriikiiriiikkkxxvAv).0(lim111riiiriiikkkxxv设 这说明当 的主特征值是实的重根时，定理5的结论还是正确的. A 应用幂法计算 的主特征值 及对应的特征向量时，

11、如果 （或 ），迭代向量 的各个不等于零A11111kv的分量将随 而趋向于无穷（或趋于零），这样在计算机实现时就可能“溢出”. k19就有得值最大的分量，规范化的绝对为向量记做改进为了避免“溢出”下面 ).()max( )max( .0 0vvvuvv(2.9) )1,2,( ./),max( , , )0( , , 131001021kanRkkkkkkknnnvuvAuvvuvA计算，对任何非零初始向量其特征值个线性无关的特征向量有设定理定理.lim ,)max(lim 111kkkkxxu则20 ,)max()max( , 对于任0011100100AvAvvvuAvAuvvu，给非零

12、向量事实上，,)max()max( ,)max( 020222200212vAvAvvuAvvAAuv.)max()max( ,)max( , 00010vAvAvvuvAvAvkkkkkkkk21,121221110nknnkkkaaaxxxvA)( )max( max)max(111212211121221100kaaaaaanknnknknnkkkkxxxxxxxxvAvAu22)( maxmax )max(max)max(111211221112122111010kaaaaaanknnknknnkkkkxxxxxxvAvAv.12确定收敛速度由比值r23A=1 1 0.5;1 1 .2

13、5;.5 .25 2u=1,1,1v=A*u,v1=max(v),u=v/v1 例例3 3 用幂法计算 0.225.05.025.00.10.1A的主特征值和相应的特征向量. 计算过程如表8-1. 245365323.2)16497.07482.0(205365374.2)16497.07482.0(195365456.2)16497.07482.0(185365598.2)16497.07482.0(175365840.2)16497.07483.0(165366256.2)16497.07483.0(155380029.2)16508.07494.0(105587918.

14、2)16674.07651.0(57500000.2)18182.09091.0(1)111(0)max(18kTkvuk（规范化向量）表25 8.2.2 8.2.2 加速方法加速方法 原点平移法原点平移法 由前面讨论知道，应用幂法计算 的主特征值的收敛速度主要由比值 来决定，但当 接近于1时，收敛可能很慢. A12rr 一个补救的办法是采用加速收敛的方法. 引进矩阵 ,pIAB其中 为选择参数. 设 的特征值为 ，则 的相应特征值为 ，而且 的特征向量相同. pAn,21Bpppn,21BA,26 如果要计算 的主特征值 ，就要适当选择 使 仍然是 的主特征值，且使 A1pp1B.1212p

15、p对 应用幂法，使得在计算 的主特征值 的过程中得到加速. 这种方法通常称为原点平移法. BBp1 例例4 4 设 有特征值 44RA),4,3 ,2, 1(15jjj比值 . 作变换 9.0/12r),12(ppIAB则 的特征值为 B.1,0, 1,2432127应用幂法计算 的主特征值 的收敛速度的比值为 B1.9.021121212pp 选择有利的 值，虽然能够使幂法得到加速，但问题在于如何选择适当的参数 . pp 设 的特征值满足 A,121nn（2.10）则不管 如何， 的主特征值为 或 .当希望计算 及 时，首先应选择 使 ppIABp1pn11xp,1ppn28且使收敛速度的比

16、值 .min,max112ppppn显然，当 , 即 时 为最小，这时收敛速度的比值为 ppppn112*22ppn.2*212112nnnpppp 当 的特征值满足（2.10）且 能初步估计时，就能确定 的近似值. An,2*p 当希望计算 时，应选择n*,211ppn29 例例5 5 计算矩阵0.225.05.025.00.10.1A的主特征值. 使得应用幂法计算 得到加速. n 作变换 取 ，则 ,pIAB75.0p.25.125.05.025.025.00.15.00.125.0B30对 应用幂法，计算结果如表8-2. B7865914.1)16497.07482.0

17、(107866587.1)16497.07483.0(97869152.1)16499.07484.0(87873300.1)16501.07488.0(77888443.1)16511.07491.0(67914011.1)16522.07516.0(5)1 1 1(0)max(28kTkvuk（规范化向量）表由此得 的主特征值为 的主特征值 为 BA,7865914.111,5365914.275.01131与例3结果比较，上述结果比例3迭代15次还好. 若迭代15次， （相应的 ). 7865258.115365258.21 原点位移的加速方法，是一个矩阵变换方法. 这种变换容易计算，又

18、不破坏矩阵 的稀疏性，但 的选择依赖于对 的特征值分布的大致了解. ApA 瑞利商加速瑞利商加速 定理定理14 14 设 为对称矩阵，特征值满足 nnA R,321n对应的特征向量满足 ，应用幂法计算 的主特征值 ，则规范化向量 的瑞利商给出 的较好的近似 ijjixx),(A1ku1.),(),(2121kkkkkOuuuAu32 证明证明 由（2.8）式及 ,)max(,)max(001100uAuAAuvuAuAukkkkkkk得 .),(),(),(),(2121122112200001knjkjjnjkjjkkkkkkkkOuAuAuAuAuuuAu（2.11）).( 1 ,333a

19、P：作作业业33三、反幂三、反幂法法反幂法可求非奇异实矩阵的按模最小特征值及特征向量。 .)max( ),max( , 2 , 1 ,)max(, ,11100kkknknkkkkknkRvvuxvvvuuAvuv计算任取非零向量. , , 11kkkkkkkyUvuLyLUuAvv求解即分解求得利用经常来求得求解可通过主元高斯消去法34 ).1,2,( ,)max( , )0( , 0 , 151100121kanRkkkkknnnnnvvuuAvvuA计算，对任何非零初始向量其特征值满足向量个线性无关的特征为非奇异矩阵且有设定理定理.1)max(lim ,)max(lim nkknnkkv

20、xxu则35特征值.移来加速迭代或求其他在反幂法中可用原点位.,1 ,1 ,1 )(21211nnppppxxxIA对应的特征向量仍然为存在，则特征值为若(2.12) ).1,2,( ,)max( ,)( 1100kpkkkkkvvuuIAvvu，计算对任何非零初始向量)( , 0 ijpppijj的近似，并且满足的特征值是如果A36.)()()(111及其对应的特征向量特征值法可求的主特征值，用上述算是则pppjjIA).( , 0 , ), 2 , 1( , 16ijpppninRijjiinn且的近似是而和值和特征向量记为的特征个线性无关的特征向量有设xAA定定理理(2.12) ).1,

21、2,( ,)max( ,)( ),0(110kpakkkkkjvvuuIAvu计算对任何非零初始向量.)max(1,1)max( ,)max( jkjkjjkppvvxxu则37 .min确定收敛速度由比值pprijij.算迭代一两次就可完成计很小，离较好，一般的较好近似且特征值分是只要rpj. ) 1 , 1 , 1 (:10110vPuLUvu，即得使得实际选取T. )( :.)( 11kkkkkkppyUvPuLyLUIAPLUuvIA，组：，需解两个三角形方程分解借助反幂法需要解方程组：38反幂法计算公式反幂法计算公式: .)max( )( 1kkkkkpvvuuvIA., 3 , 2

22、 ,/ ),max( , (2) ;/),max( , ) 1 , 1 , 1 (1 . 2.,)( : . 111111111kpkkkkkkkkkTvuvyUvPuLyvuvvUvULPLUIAPLU得）解（反幂法迭代，保存分解.2679. 133410131012 3的特征向量的特征值求A例例6 6format long;A=2 1 0;1 3 1;0 1 4,p=1.2679,B=A-p*eye(3);L U P=lu(B);L,U,P,v=U1 1 1, mu=max(v);u=v/mu,v=U(L(P*u), mu=max(v);u=v/mu,lamda=p+1/mu393 Hou

23、seholder3 Householder方法方法一、引言一、引言本节讨论两个问题:. ) 1 ( 伯格矩阵海森似变换约化实矩阵为上用初等反射阵作正交相. )2( 对角矩阵的问题对称似变换约化对称矩阵为用初等反射阵作正交相. 值问题或对称对角矩阵的特征格矩阵题就转化为求上海森伯于是，原矩阵特征值问已经学过的知识：40介绍初等反射阵和平面旋转阵，一、初等反射阵一、初等反射阵阶矩阵则且设nRTn , 1 , www定义9定义9TwwIH2 .rHouseholde,)(矩阵也称矩阵或镜面反射称为初等反射.212222122221)( ,),(221222121212121nnnnnTnwwwwww

24、wwwwwwwwwwwwwHw则若记41.2)( , 022uuuIuHuT都有初等反射阵对于因此 , ,., (4) ),( , (3) , , (2) , , (1) , 1 ,2 112亦是对称矩阵则是对称矩阵若或即是对合矩阵即是正交矩阵即是对称矩阵则其中设初等反射矩阵HAHAHHAHHIHIHHHHHHwwwwIHTTTT定定理理2 25 5. , , 02222xyHxyx都有令因此，对于42, ,0| SySxyxRRSnTn其中有则对设超平面：vvxwx几几何何意意义义vSv. ).|(2)2(,02)2(,yxyyyyyyyxxxxTTHvwwwIHwwwIH于是. , , ,

25、 yHxHxyyx使得反射阵则总存在初等且互异设向量nR定定理理2 26 6xyS43 , ,)0 , 0 , 1 (0),( 11121，使得则总存在初等反射阵，单位向量设约化定理eHxuuIHexTTTnxxx) )定定理理2 27 7( (其中 . )(,),( ,)sgn(1222121121xxxxxTnuexux）（则设*,)(),(,121uauIHaHaHaHaHAAiTinnmR.1)( ,1次乘法次除法用次乘法用计算mmiTiTuauau 2次乘法用 mnHA44. , ,) 1 , 2 , 2(的后两个分量为零使得确定初等反射阵给定向量HxHxT : :练练习习.151

26、,15,) 1 , 2 , 5(,)0 , 0 , 3(, 3 :151415231152151132313232TTTuuIHuHx答45.规范化：).218(P算算法法6 646.cossinsincos)( ,cossinsincos , ,21212为正交矩阵换，其中是平面上向量的旋转变或则变换设向量PPxyyxxxyyR二、平面旋转矩阵二、平面旋转矩阵47.)(11cossin11sincos11),( , 变换或平面上的旋转为为平面旋转矩阵，则称设GivensxxxjiRjinPyPyx定义定义48);,( , ,),( 2)(1),(),(1jikxycxsxysxcxyjiPP

27、jijikkjijjiiT，则设）（；）正交（的性质和作用：xPyPPP. 3APPA，）（./,/,), 0 ,( ),( ,),( 2211ijiijiiTnijiTnjixxsxxcxxxxxxjixxxxxx其中使则可选择平面旋转矩阵不全为零其中，设约化定理PxPx) )定定理理2 28 8( (49阶矩阵则且设nRTn , 1 , www定义定义TwwIH2 .rHouseholde,)(矩阵也称矩阵或镜面反射称为初等反射都有初等反射阵对于因此 , 0u , ,.2)( 22uuuIuHT , ,)0 , 0 , 1 (0),( 11121，使得则总存在初等反射阵，设约化定理eHxu

28、uIHexTTTnxxx) )定定理理( (其中 . )(,),( ,)sgn(1222121121xxxxxTnuexux50 般矩阵为上海森伯格阵用正交相似变换约化一二二、 .),(,1131211) 1 (221) 1 (12112122221112111TnnnnnnnaaaaaaaaaaaaacAcAAA其中步约化：设第 ),1否则这步不需做不妨0(c ,111111111，使得取初等反射阵ecRuuIRT其中 . )(,),( ,)sgn(2111221211131121111121211aaaaaTnuecuc51,00 )2(222)2(12)2(11)2()2(3)2(2)2

29、(3)2(33)2(32)2(2)2(23)2(221)2(1)2(13)2(12111) 1 (221111) 1 (12111112AcAARARcRRAUAUA 0nnnnnnnaaaaaaaaaaaaaa，则令111RU.),()2(2)2(322Tnaac其中52)()(1,)(,)(, 1)(1, 1)(, 1)(,1) 1(1, 12)(2)(1, 2)(2) 1(1, 2)2(221)(1)(1, 1)(1) 1(1, 1)2(12) 1 (1111111111121 ,1, 1knnkknkknknkkkkkkkkkkkkkkknkkkkkkknkkkkkkkkkkkkkaa

30、aaaaaaaaaaaaaaaaakkUUAUUUAUAUUU使得步，即有步约化：假定已完成第第, )(22)(12)(11knkkkkkAcAA 0.)(11阶上海森伯格矩阵是其中kkA53 ,11，使取初等反射阵不妨ecRuuIRckkkTkkkkk0其中 . )(,),( ,)sgn()(, 12221)(,)(, 2)(, 112)(, 1kkkkkkkTkknkkkkkkkkkkkkkkkaaaaauecuc. ) 1(221) 1(12) 1(11)(22)(12)(111kkkkkkkkkkkkkkkkkkkAcAARARcRRAAUAUARIU0000，则令.1) 1(11阶上

31、海森伯格矩阵是其中kkA. )(22)(22)(12）为对称阵时只需计算（当及步约化只需计算第RARARARRAkkkkkkk54.* 2) 1(1)2(1, 1)3(332)2(221) 1 (1121121nnnnnnnnnnaaaaanUAUUUA步约化：第).( , 002112221上海森伯格矩阵使得，则存在初等反射阵设HAUUUAUUUUUUATnnnnnR定定理理1 17 755.32 12,1, . 2; . 10010kkkkkkkkkknkIUUURIUAUUAecRRU）（；，其中）约化计算（；，使）计算初等反射阵（对：00算法1.724232734 1约化为上海森伯格矩

32、阵用豪斯霍尔德方法将 AA例例7 7 ,)4 , 2(111111，使得，求初等反射阵解：ecRRcT56.0.4472 0.8944- 0.8944- 0.4472- 16)522(4)522(4)522()522(52111 )522(52)( ,)4 , 522( , 5220)sgn(21111121111111121211TTaauuIRecuc，.2.2000 0.4000- 0.0000 0.4000- 7.8000 4.4721-0.4472- 7.6026 4.0000-1112HUAUA，111RU57 称阵为对称三对角阵用正交相似变换约化对三、. , 8111222111

33、2112221CUAUUUUUUAnnnnnnnnncbbcbbcbbcR使得反射阵为对称阵，则存在初等设定定理理1 1.17知：由定理证明58, .11) 1(221) 1(12) 1(11)(22)(12)(111knkkkkkknkkkkkkkkkkkkkkAcAARARcRRAAUAUA00步约化中，在第下面考虑实现过程).()(22即可实际上对角线以下元素和对称，故只需计算因kkkkRARRA,)(21,1)(221knkkTkkkkknkkkkRRururtuAr记)( )(221)(221)(22TkkkkkTkkkkkkuuAAuuIRAR则TkkTkkkTkkkTkkkTkk

34、TkkkkkTkkTkkkTkkTkkkTkkTkkkkTkkkTkkkuttuAuururururuAuruuruurAuruuAuuurA)(2211)(221)(221)(221)(22 )(2)(2 )( 59 算法算法2 2（豪斯霍尔德约化对称阵为对称三对角阵） 设 为对称阵，本算法确定初等反射阵 使 （为对称三对角阵）. nnA R)2, 1(nkUkCUAUUUnn2112 的对角元 存放在数组 内， 的次对角元素 存放在数组 内. C), 1(nici)(ncC) 1, 1(nibi)(nb 数组 最初可用来存放 及 ，确定 中向量 的分量存放在 的相应位置. 冲掉 ，约化 的

35、结果冲)(nbkrktkRkuAkkkaA掉 ，数组 的上部分元素不变. 如果第 步不需要变换则置 为零. AAkk60dbuauuuunkidauabadadRacnkkkkkkkkkkkknkiikkkikikikkkkkiknikkkkk)7)6)5)sgn()4), 1(/)3)4,0,0,0)2max)1)2()1(2,2, 1.1,1,1,112,11计算转则如果计算确定变换对于61jkijikijijkkkkkikiikiknijjkjiikjjkijikTkkkubbuaaikjnkiRARnkiusbbtuhssuauahbnkiruuAs*, 1, 1)4), 1(/)2/()()3)b()a(, 1)201)3()(2211)(22对于对角线以下部分计算计算对于及计算）应用变换621,11,11.2)5nnnnnnnnnabacack继续循环 对对称阵 用初等反射阵正交相似约化为对称三对角阵大约需要 次乘法. A332n634 QR4 QR方法方法Rutishauser(1958)利用矩阵的三角分解提出计算矩阵特征值的LR算法，Francis(1961,1962)利用矩阵的QR分解建立计算矩阵特征值的QR方法.QR方法是一种变换方法，是