第8章矩阵特征值计算

1、上页上页下页下页第第8章章 矩阵特征值计算矩阵特征值计算 8.1 特征值性质和估计特征值性质和估计 8.2 幂法及反幂法幂法及反幂法 8.3 正交变化与矩阵分解正交变化与矩阵分解 8.4 QR方法方法上页上页下页下页8.1 特征值性质和估计特征值性质和估计 工程技术中有多种振动问题，如桥梁或建筑物的工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题征向量的问题. 下面先复习一些矩阵的特征值和特征向量的基础下面

2、先复习一些矩阵的特征值和特征向量的基础知识及性质知识及性质.上页上页下页下页8.1.1 特征值问题及性质特征值问题及性质求求A的特征值问题的特征值问题(1.1)等价于求等价于求A的特征方程的特征方程设矩阵设矩阵A Rnn，特征值问题是求，特征值问题是求 C和非零向和非零向量量x Rn，使，使( )det()0(1.2)pIA的根的根. 在在5.1.3节已经给出特征值的一些性质，下面节已经给出特征值的一些性质，下面再补充一些基本性质再补充一些基本性质. Ax= x (1.1)其中其中 是矩阵是矩阵A的的特征值特征值，x是矩阵是矩阵A属于特征值属于特征值 的的特征向量特征向量. 本章讨论计算矩阵特

3、征值的数值方法本章讨论计算矩阵特征值的数值方法, 在在科学和工程技术中很多问题在数学上都归结为求矩科学和工程技术中很多问题在数学上都归结为求矩阵的特征值问题阵的特征值问题.上页上页下页下页 定理定理1 设设 为为ARnn的特征值的特征值, 且且Ax= x, x 0，则有则有 (2) - - 为为A- - I的特征值，即的特征值，即(A- - I)x=( - - )x ； (1) c 为的为的cA特征值特征值(c 0为常数为常数)； (3) k为为Ak的特征值，即的特征值，即Akx= kx ； (4) 设设A为非奇异矩阵，那么为非奇异矩阵，那么 0 , 且且 - -1为为A- -1的特的特征值，

4、即征值，即A- -1x= - -1x .上页上页下页下页 定理定理2 (1) 设矩阵设矩阵ARnn可对角化，即存在非可对角化，即存在非奇异矩阵奇异矩阵P使使的充分必要条件是的充分必要条件是A具有具有n个线性无关的特征向量个线性无关的特征向量. (2) 如果矩阵如果矩阵ARnn有有m个个(mn)不同的特征不同的特征值值 1, 2, m，则对应的特征向量，则对应的特征向量 x1, x2, , xm 线线性无关性无关.121,nPAP 上页上页下页下页 定理定理3(对称矩阵的正交约化对称矩阵的正交约化) 设设ARnn为对称为对称矩阵，则矩阵，则 (3) 存在一个正交矩阵存在一个正交矩阵P使的使的且且

5、 1, 2, n为为A的特征值，而的特征值，而P(u1,u2,un)的列的列向量向量uj为为A的对应于的对应于 j 的单位特征向量的单位特征向量.12,TnP AP (1) A的特征值均为实数；的特征值均为实数； (2) A有有n个线性无关的特征向量；个线性无关的特征向量；上页上页下页下页 定义定义 设设ARnn为对称矩阵，对于任一非零向为对称矩阵，对于任一非零向量量x 0，称，称,),(),()(xxxAxxR 为对应于向量为对应于向量x的的瑞利瑞利(Rayleigh)商商. 定理定理4 设设ARnn为对称矩阵为对称矩阵(其特征值依次记其特征值依次记为为 1 2 n)，则，则 (1). (对

6、任何非零对任何非零xRn)；1(, )( , )nAx xx x (2). ； (1.3)10(, )max( , )nx RxAx xx x 0(, )min.( , )nnx RxAx xx x 上页上页下页下页 证明证明 只证只证(1)，关于，关于(2)自己作练习自己作练习. 由于由于A为实对称矩阵，可将为实对称矩阵，可将 1, 2, n 对应的特对应的特征向量征向量 x1,x2,xn 正交规范化，则有正交规范化，则有 (xi, xj)=ij，设，设x 0为为Rn中任一向量，则有中任一向量，则有. 0,1221 niiniiixxx 于是于是.),(),(11212 niiniiinxx

7、xAx从而从而(1)成立成立. 结论结论(1)说明说明瑞利商瑞利商必位于必位于 n和和 1之间之间.上页上页下页下页8.1.2 特征值估计与扰动特征值估计与扰动 定义定义1 设设n阶矩阵阶矩阵A=(aij)，令，令).,(niarnijjiji211 (1) ； (2) 集合集合 称为复平面上以称为复平面上以aii为圆心，以为圆心，以ri为半径的为半径的n阶矩阵阶矩阵A的的n个个格什戈林格什戈林(Gerschgorin)圆盘圆盘. |,(1,2, )iiiiDzzarzCin 上页上页下页下页 定理定理5 (Gerschgorin圆盘定理圆盘定理) 特别地，如果特别地，如果A的一个圆盘的一个圆

8、盘Di是与其它圆盘分离是与其它圆盘分离(即即孤立圆盘孤立圆盘)，则，则Di中精确地包含中精确地包含A的一个特征值的一个特征值.1(1,2,. ).(1.4)niiiijjj iarain (1) 设设n阶矩阵阶矩阵A(aij)，则，则A的每一个特征值必属的每一个特征值必属于下面某个圆盘之中于下面某个圆盘之中 (2) 如果如果A有有m个个圆盘圆盘组成一个连通的并集组成一个连通的并集S，且，且S与余下与余下n- -m个个圆盘圆盘是分离的，则是分离的，则S内恰包含内恰包含A的的m个个特征值特征值.或者说或者说 A的特征值都在的特征值都在n个圆盘的个圆盘的并集并集中中.上页上页下页下页 证明证明 只就

9、只就(1)给出证明给出证明. 设设 为为A的特征值，即的特征值，即.1knkjjkjkkraa Ax= x，其中，其中 x=(x1,x2, xn)T 0.或或 记记 ，考虑，考虑Ax= x的第的第k个个方程，即方程，即0max1 xxxinik,1knjjkjxxa ,)( kjjkjkkkxaxa 于是于是即即, kjkjkkjjkjkkkaxxaxa 上页上页下页下页这说明，这说明，A的每一个特征值必位于的每一个特征值必位于A的一个圆盘的一个圆盘中，并且相应的特征值中，并且相应的特征值 一定位于第一定位于第k个圆盘中个圆盘中(其中其中k是对应特征向量是对应特征向量x绝对值最大的分量的下标绝

10、对值最大的分量的下标).利用相似矩阵性质，有时可以获得利用相似矩阵性质，有时可以获得A的特征值进的特征值进一步的估计，即适当选取非奇异对角阵一步的估计，即适当选取非奇异对角阵,112111 nD 并做相似变换并做相似变换 .适当选取适当选取 可使某些圆盘半径及连通性发生变化可使某些圆盘半径及连通性发生变化.nnijijaADD 1), 2 , 1(nii 上页上页下页下页.411101014 A 例例1 估计矩阵估计矩阵A的特征值范围，其中的特征值范围，其中 解解 矩阵矩阵A的的3个圆盘为个圆盘为. 24:, 2:, 14:321 DDD 由定理由定理8，可知，可知A的的3个特征值位于个特征值

11、位于3个圆盘的并个圆盘的并集中，由于集中，由于D1是孤立圆盘，所以是孤立圆盘，所以D1内恰好包含内恰好包含A的的一个特征值一个特征值 1(为实特征值为实特征值)，即，即. 531 A的其它两个特征值的其它两个特征值 2, 3包含在包含在D2, D3的并集中的并集中.上页上页下页下页现在取对角阵现在取对角阵,9 . 0000100011 D做相似变换做相似变换.49 . 09 . 00101491011 ADDAA 矩阵矩阵A1的的3个圆盘为个圆盘为. 8 . 14:,919:, 14:321 EEE上页上页下页下页显然，显然，3个圆盘都是孤立圆盘，所以，每一个圆个圆盘都是孤立圆盘，所以，每一个

12、圆盘都包含盘都包含A的一个特征值的一个特征值(为实特征值为实特征值)，且有估计，且有估计 . 2 . 28 . 5,919919, 53321 上页上页下页下页 定理定理6(Bauer-Fike定理定理) 设设 是是A+IRnn的一个的一个特征值，且特征值，且P- -1AP=D=diag( 1, 2, n)，则有，则有1()min,(1.5)pppAPPI 其中其中|p为矩阵的为矩阵的p范数，范数，p=1,2, . 证明证明 由于由于 (A)时显然成立时显然成立,故只考虑故只考虑 (A).这这时时D- - I非奇异，设非奇异，设x是是A+I对应于对应于 的特征向量，由的特征向量，由(A+I-

13、- I)x=0左乘左乘P- -1可得可得111()()()(),DIPxPIPPx 1111() ()(),PxDIPIPPx 而对角矩阵而对角矩阵(D- - I)- -1的范数为的范数为1()1(),min,pADImm 所以有所以有1.pppPIPm 这就得到这就得到(1.5)式式. 这时总有这时总有中的一个中的一个 取到取到m值值.上页上页下页下页 由定理由定理6可知可知|P- -1| |P|=cond(P)是特征值扰动的是特征值扰动的放大系数，但将放大系数，但将A对角化的相似变化矩阵不是唯一的，对角化的相似变化矩阵不是唯一的，所以取所以取cond(P)的下确界的下确界 112()inf

14、()(,) ,(1.6)nAcond P PAPdiag 称为特征值问题的条件数称为特征值问题的条件数. 只要只要v(A)不很大，矩阵微不很大，矩阵微小扰动只带来特征值的微小扰动小扰动只带来特征值的微小扰动. 但是但是v(A)难以计算难以计算,有时只对一个有时只对一个P，用，用cond(P)代替代替v(A).特征值问题的条件数和解线性方程组时的矩阵是特征值问题的条件数和解线性方程组时的矩阵是两个不同的概念，对于一个矩阵两个不同的概念，对于一个矩阵A，两者可能一大一，两者可能一大一小，例如二阶矩阵小，例如二阶矩阵A=diag(1, 10- -10)，有，有v(A) =1，但解，但解线性方程组的矩

15、阵条件数线性方程组的矩阵条件数cond(A)=1010.上页上页下页下页 关于计算矩阵关于计算矩阵A的特征值问题，当的特征值问题，当n=2,3时，我们时，我们还可按行列式展开的办法求特征方程还可按行列式展开的办法求特征方程p( )=0的根的根. 但但当当n较大时，如果按展开行列式的办法较大时，如果按展开行列式的办法, 首先求出首先求出p( )的系数，再求的系数，再求p( )的根，工作量就非常大，用这种办的根，工作量就非常大，用这种办法求矩阵特征值是不切实际的，由此需要研究法求矩阵特征值是不切实际的，由此需要研究A的特的特征值及特征向量的数值方法征值及特征向量的数值方法. 本章将介绍一些计算机上

16、常用的两类方法，一类本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法是幂法及反幂法(迭代法迭代法)，另一类是正交相似变化的，另一类是正交相似变化的方法方法(变化法变化法).上页上页下页下页8.2 幂法及反幂法幂法及反幂法 幂法与反幂法都是求幂法与反幂法都是求实矩阵实矩阵的特征值和特征向的特征值和特征向量的量的向量迭代法向量迭代法，所不同的是幂法是计算矩阵的，所不同的是幂法是计算矩阵的主主特征值特征值(矩阵矩阵按模最大的特征值按模最大的特征值称为称为主特征值主特征值，其模，其模就是该矩阵的就是该矩阵的谱半径谱半径)和和相应特征向量相应特征向量的一种向量迭的一种向量迭代法，特别适用于大型稀

17、疏矩阵代法，特别适用于大型稀疏矩阵. 而反幂法则是计而反幂法则是计算非奇异算非奇异(可逆可逆)矩阵矩阵按模最小的特征值按模最小的特征值和和相应特征相应特征向量向量的一种向量迭代法，特别是计算海森伯格阵或的一种向量迭代法，特别是计算海森伯格阵或三对角阵的对应一个给定近似特征值的特征向量的三对角阵的对应一个给定近似特征值的特征向量的有效方法之一有效方法之一. 下面分别介绍幂法与反幂法下面分别介绍幂法与反幂法. 上页上页下页下页8.2.1 幂法幂法(又称又称乘幂法乘幂法)现讨论求现讨论求 1及及x1的方法的方法.), 2 , 1(nixAxiii 设实矩阵设实矩阵A=(aij)有一个有一个完全的特征

18、向量组完全的特征向量组，即，即A有有n个线性无关的特征向量个线性无关的特征向量，设矩阵，设矩阵A的特征值为的特征值为 1, 2, n, 相应的特征向量为相应的特征向量为x1,x2,xn. 已知已知A的的主特征值主特征值 1是实根，且满足条件是实根，且满足条件)1 . 2(|,|21n 上页上页下页下页)3 . 2(),0(122110 axaxaxavnn设设 幂法的幂法的基本思想基本思想是是: 任取任取非零非零的初始向量的初始向量v0 , 由矩由矩阵阵A构造一向量序列构造一向量序列vk称为迭代向量，由假设，称为迭代向量，由假设，v0可唯一表示为可唯一表示为)2 . 2(.,.,0110212

19、01 vAAvvvAAvvAvvkkk上页上页下页下页)3 . 2(),0(122110 axaxaxavnn设设于是于是).()/(1112111122211101kkniikiiknknnkkkkkxaxaxaxaxaxavAAvv 其中其中.)/(21 niikiikxa 由假设由假设 故故 从而从而), 3 , 2( 1/1nii , 0lim kk )4 . 2(.lim111xavkkk 为为 1的特征向量的特征向量.上页上页下页下页所以当所以当k充分大时，有充分大时，有)5 . 2(,111xavkk 即为矩阵即为矩阵A的的对应特征值对应特征值 1 的一个近似特征向量的一个近似特

20、征向量.用用(vk)i 表示表示vk的第的第i个分量，则当个分量，则当k充分大时，有充分大时，有 )7 . 2(.11 ikikvv即为即为A的的主特征值主特征值 1的近似值的近似值.)6 . 2(,111111kkkkvxaAvv 由于由于 这种由已知非零向量这种由已知非零向量v0及矩阵及矩阵A的乘幂的乘幂Ak构造向构造向量序列量序列vk以计算以计算A的的主特征值主特征值 1(利用利用(2.7)式式)及相应及相应特征向量特征向量(利用利用(2.5)式式)的方法就称为的方法就称为(乘乘)幂法幂法.上页上页下页下页 迭代公式实质上是由矩阵迭代公式实质上是由矩阵A的乘幂的乘幂 Ak与非零向与非零向

21、量量v0相乘来构造向量序列相乘来构造向量序列vk=Akv0，从而计算主特，从而计算主特征值征值 1及其对应的特征向量，这就是及其对应的特征向量，这就是幂法幂法的思想的思想. ).(11 kvvikik 的收敛速度由比值的收敛速度由比值,12 r来确定，来确定，r 越小收敛越快，但当越小收敛越快，但当r 1 1时收敛可能很慢时收敛可能很慢.总结上述讨论，有下面的定理总结上述讨论，有下面的定理.上页上页下页下页 定理定理7 设设ARnn有有n个线性无关的特征向量，个线性无关的特征向量，主特征值主特征值 1满足满足 | 1| 2| | n|，则对任何非零向量则对任何非零向量v0(a1 0)，(2.4

22、)式和式和(2.7)式成立式成立.又设又设A有有n个线性无关的特征向量，个线性无关的特征向量， 1对应的对应的r个线性个线性无关的特征向量为无关的特征向量为x1,x2,xr，则由，则由(2.2)式有式有 如果如果A的主特征值为实的重根的主特征值为实的重根, 即即 1= 2= r, 且且 | r| r+1| | n|，,)/(11110 nriikiiriiikkkxaxavAv 上页上页下页下页).0(lim111 riiiriiikkkxaxav设设 为为 1的特征向量，这说明当的特征向量，这说明当A的主特征值是实的重根的主特征值是实的重根时，定理时，定理7的结论还是正确的的结论还是正确的.

23、 应用应用幂法幂法计算计算A的主特征值的主特征值 1及其对应的特征向及其对应的特征向量时，如果量时，如果| 1|1(或或| 1| 2| | n|，则对任意非零初始，则对任意非零初始向量向量v0=u0(a1 0)，有幂法计算公式为，有幂法计算公式为(uk,vk)则有则有 (1),)max(lim11xxukk (2).lim1 kk上页上页下页下页例例2 用幂法计算矩阵用幂法计算矩阵 0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10 . 1A的主特征值和相应的特征向量的主特征值和相应的特征向量.解解取取 v0=u0=(1,1,1)T, 则则 TTvuv1 ,81

24、82. 0,9091. 01,75. 2,75. 2 ,25. 2 , 5 . 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上页上页下页下页计算结果如下表计算结果如下表k0(1, 1, 1) 1(0.9091, 0.8182, 1) 2.75000005(0.7651, 0.6674, 1)2.558791810(0.7494, 0.6508, 1)2.538002915(0.7483, 0.6497, 1)2.536625616(0.7483, 0.6497, 1)2.536584017(0.7482, 0.6497, 1)2.536559818(0.7482,

25、 0.6497, 1)2.536545619(0.7482, 0.6497, 1)2.536537420(0.7482, 0.6497, 1)2.5365323)(规规范范化化向向量量Tkumaxkkv 上页上页下页下页,5365323. 21 这个结果是用这个结果是用8位浮点数字进行运算得到的位浮点数字进行运算得到的, uk的分量值是舍入值的分量值是舍入值. 于是得到于是得到及相应的特征向量及相应的特征向量(0.7482, 0.6497, 1)T. 1和相应的特和相应的特征向量的真值征向量的真值(8位数字位数字)为为,5365258. 21 .)1,64966116. 0,74822116.

26、 0(1Tx 上页上页下页下页例例 用幂法计算矩阵用幂法计算矩阵 1634310232A的主特征值与其对应的特征向量的主特征值与其对应的特征向量.解解取取 v0=u0=(0,0,1)T , 则则 TTvuv25. 0 , 1, 5 . 01, 4,1 , 4 , 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上页上页下页下页直到直到k=8 时的计算结果见下表时的计算结果见下表k1 2, 4, 1, 4 0.5, 1, 0.252 4.5, 9, 7.75 90.5, 1, 0.86113 5.7222, 11.4444, 8.36111.44440.5, 1, 0.

27、73604 5.4621, 10.9223, 8.2306 10.92230.5, 1, 0.75365 5.5075, 11.0142, 8.2576 11.01420.5, 1, 0.74946 5.4987, 10.9974, 8.2494 10.99740.5, 1, 0.75017 5.5002, 11.0005, 8.2501 11.00050.5, 1, 0.75008 5.5000, 11.0000, 8.2500 11.00000.5, 1, 0.7500TkvTku Tx7500. 0,0 . 1,5 . 0,0000.1111 从而从而k 上页上页下页下页8.2.2 加速

28、方法加速方法1、原点平移法、原点平移法 由前面讨论知道，应用幂法计算由前面讨论知道，应用幂法计算A的主特征值的的主特征值的收敛速度主要由比值收敛速度主要由比值 r=| 2/ 1|来决定，但当来决定，但当r 接近于接近于1时，收敛可能很慢时，收敛可能很慢. 这时，一个补救办法是采用加速这时，一个补救办法是采用加速收敛的方法收敛的方法.其中其中p为参数，设为参数，设A的特征值为的特征值为 i，则对矩阵，则对矩阵B的特征的特征值为值为 i- -p ，而且，而且A, B的特征向量相同的特征向量相同. 引进矩阵引进矩阵 B=A- -pI .上页上页下页下页 如果要计算如果要计算A的主特征值的主特征值 1

29、, 只要只要选择合适的数选择合适的数p，使使 1- -p为矩阵为矩阵B=A- -pI 的主特征值，且的主特征值，且 1212max ppini那么，对矩阵那么，对矩阵B=A- -pI应用幂法求其主特征值应用幂法求其主特征值 1- -p, 收收敛速度将会加快敛速度将会加快. 这种通过求这种通过求B=A- -pI的主特征值和特的主特征值和特征向量，而得到征向量，而得到A的主特征值和特征向量的方法叫的主特征值和特征向量的方法叫原原点平移法点平移法. 对于对于A的特征值的某种分布，它是十分有的特征值的某种分布，它是十分有效的效的.上页上页下页下页例例3 设设AR44有特征值有特征值),4 , 3 ,

30、2 , 1(15 jji 比值比值r=| 2/ 1| 0.9. 做变换做变换B=A- -12I (p=12),则则B的特征值为的特征值为. 1, 0, 1, 24321 应用幂法计算应用幂法计算B的主特征值的主特征值1的收敛速度的比值为的收敛速度的比值为. 9 . 021121212 pp 虽然常常能够选择有利的虽然常常能够选择有利的p值值, 使幂法得到加速使幂法得到加速, 但设计一个自动选择适当参数但设计一个自动选择适当参数p的过程是困难的的过程是困难的. 1212max142 ii上页上页下页下页 下面考虑当下面考虑当A的特征值是实数时，怎样选择的特征值是实数时，怎样选择p使采使采用幂法计

31、算用幂法计算 1得到加速得到加速.ppn 1且使且使收敛速度的比值收敛速度的比值.min,max112 ppppn 设设A的特征值都是实数，且满足的特征值都是实数，且满足)10. 2(,121nn 则对实数则对实数p，使矩阵，使矩阵A- -pI的主特征值为的主特征值为 1- -p或或 n- -p时，时，当当我们计算我们计算 1及及x1时，首先应选取时，首先应选取p使使上页上页下页下页显然，当显然，当 2- -p=- -( n- -p )时，即时，即P=( 2+ n)/2P*时时 为为最小值，这时最小值，这时收敛速度的比值收敛速度的比值为为.2212112nnnpppp 当希望计算当希望计算 n

32、时，应选取时，应选取 p=( 1+ n-1)/2P* 使得应用幂法计算使得应用幂法计算 n得到加速得到加速. 当当A的特征值都是实数，满足的特征值都是实数，满足(2.10)式式且且 2, n能初步估计出来，我们就能确定能初步估计出来，我们就能确定P*的近似值的近似值.nn 121上页上页下页下页 例例4 用原点平移加速法求用原点平移加速法求例例2中矩阵中矩阵A的主特征值的主特征值与其对应的特征向量与其对应的特征向量.25. 125. 05 . 025. 025. 00 . 15 . 00 . 125. 0 B,0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10

33、. 1 A对对B应用幂法，仍应用幂法，仍取取 v0=(1,1,1)T , 则则 .1,75. 0,875. 01, 2,2, 5 . 1,75. 111111TTvuv ), 2 , 1(max1 k/vuvBuvkkkkkkk 解解 取取p=0.75, 做平移变换做平移变换B=A- -pI，则，则上页上页下页下页k0(1, 1, 1) 1(0.8750, 0.7500, 1) 2.00000005(0.7516, 0.6522, 1)1.79140116(0.7491, 0.6511, 1)1.78884437(0.7488, 0.6501, 1)1.78733008(0.7484, 0.6

34、499, 1)1.78691529(0.7483, 0.6497, 1)1.786658710(0.7482, 0.6497, 1)1.7865914)(规规范范化化向向量量Tkumaxkkv 计算结果如下表计算结果如下表由此得由此得B的主特征值为的主特征值为 1 1.7865914，A的主特征值的主特征值 1为为 1 = 1 +0.75 2.5365914.上页上页下页下页 原点位移的加速方法，是一个矩阵变换方法原点位移的加速方法，是一个矩阵变换方法. 这这种变换容易计算，又不破坏矩阵种变换容易计算，又不破坏矩阵A的稀疏性，但的稀疏性，但p的的选择选择依赖依赖对对A的的特征值分布的大致了特征

35、值分布的大致了解解.及相应的特征向量为及相应的特征向量为 x1 (0.7482, 0.6497, 1)T. 与例与例2结果比较，上述结果比例结果比较，上述结果比例2迭代迭代15次还好次还好. 若若迭代迭代15次次 1 1.7865258(相应的相应的 1 2.5365258)就是真就是真值值.上页上页下页下页 设设ARnn为为对称矩阵对称矩阵，称，称.),(),()(xxxAxxR 为向量为向量x的的瑞利商瑞利商，其中，其中(x, x)=xTx为内积为内积. 由定理由定理4知道，实对称矩阵知道，实对称矩阵A的特征值的特征值 1及及 n可用可用瑞利商瑞利商的的极限值表示极限值表示. 下面我们将下

36、面我们将瑞利商瑞利商应用到用幂法计算应用到用幂法计算实对称矩阵实对称矩阵A的主特征值的加速上来的主特征值的加速上来.2、瑞利商、瑞利商(Rayleigh)加速加速上页上页下页下页 定理定理9 设设ARnn为为对称矩阵对称矩阵，特征值满足，特征值满足对应的特征向量对应的特征向量xi满足满足(xi, xj)= ij (单位正交向量单位正交向量), 应应用幂法公式用幂法公式(2.9)计算计算A的主特征值的主特征值 1，则规范化向，则规范化向量量uk的的瑞利商瑞利商给出给出 1的较好的近似值为的较好的近似值为,121nn kkkkkkuuuAuuR2121, 由此可见，由此可见，R(uk)比比 k更快

37、的收敛于更快的收敛于 1.上页上页下页下页 证明证明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(2121122112200001 knjkjjnjkjjkkkkkkkkkaauAuAuAuAuuuAuuR 得得上页上页下页下页 幂法的幂法的瑞利商加速迭代公式瑞利商加速迭代公式可以写为可以写为 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A为为n阶实对称矩阵阶实对称矩阵.,11kkux 对给定的误差限对给定的误差限 ，当，当| k k- -1| 时，取近

38、似值时，取近似值上页上页下页下页8.2.3 反幂法反幂法 反幂法是用于求非奇异矩阵反幂法是用于求非奇异矩阵A的的按模最小的特征按模最小的特征值和对应特征向量值和对应特征向量的方法的方法. 而结合原点平移法的反幂而结合原点平移法的反幂法则可以求矩阵法则可以求矩阵A的任何一个的任何一个给定近似特征值对应的给定近似特征值对应的特征向量特征向量.设矩阵设矩阵A非奇异非奇异,其特征值其特征值 i (i=1,2,n) ,满足满足0121 nn 其相应的特征向量其相应的特征向量x1,x2,xn线性无关，则线性无关，则 A- -1 的特征的特征值为值为1/ i ，对应的特征向量仍为，对应的特征向量仍为xi (

39、i=1,2,n).iiiiiixxAxAx11 上页上页下页下页此时此时, A- -1的特征值满足的特征值满足11111 nn因此因此, 对对A- -1应用幂法应用幂法,可求出其可求出其主特征值主特征值 1/ n k 和和特征向量特征向量 xn uk .从而求得从而求得A的的按模最小特征值按模最小特征值 n 1/ k 和对应的和对应的特征向量特征向量 xn uk ，这种求这种求A- -1的方法就称为的方法就称为反幂法反幂法.上页上页下页下页为了避免求为了避免求A- -1, 可通过解线性方程组可通过解线性方程组Avk=uk- -1得到得到vk，采用采用LU分解法，即先对分解法，即先对A进行进行L

40、U分解分解A=LU, 此时此时反幂反幂法的迭代公式法的迭代公式为为 , 2 , 1/max,1 kvuvvzUvzuLzkkkkkkkkkkk 求求出出解解求求出出解解), 2 , 1()max(11 k/vuvuAvkkkkkkk knknux ,1 反幂法的迭代公式反幂法的迭代公式为为:任取任取v0=u0 0 0，构造，构造上页上页下页下页对给定的误差对给定的误差 ，当，当|kk- -1| n|0，则对任意非零初始向量则对任意非零初始向量u0(an 0) ，由反幂法计算公，由反幂法计算公式构造的向量序列式构造的向量序列vk，uk满足满足 (1)lim,max()nkknxux (2).1)

41、max(limnkkv 收敛速度收敛速度的比值的比值r=| n/ n- -1|. 上页上页下页下页 在反幂法中也可以用在反幂法中也可以用原点平移法原点平移法加速迭代过程加速迭代过程，或或求其它特征值与其对应的特征向量求其它特征值与其对应的特征向量. 如果矩阵如果矩阵(A- -pI)- -1存在，显然其特征值为存在，显然其特征值为,1,1,121pppn 对应的特征向量仍然是对应的特征向量仍然是x1,x2,xn，现对矩阵，现对矩阵(A- -pI)- -1应用幂法，得到反幂法的迭代公式应用幂法，得到反幂法的迭代公式00110,(),(1,2,.).(2.12)/ max()kkkkkuvvApIu

42、kuvv 初初始始向向量量上页上页下页下页 如果如果p是是A的特征值的特征值 j的一个近似值，且设的一个近似值，且设 j与其它与其它特征值是分离的，即特征值是分离的，即),(jippij 就是说就是说1/( j- -p)是矩阵是矩阵 (A- -pI)- -1的主特征值，可用反幂的主特征值，可用反幂法法(2.12)计算特征值计算特征值 j及特征向量及特征向量xj. 设设ARnn有有n个线性无关的特征向量个线性无关的特征向量x1,x2, xn,则则),0(10 jniiiaxau,)max()(0)1(0upIAupIAvkkk 上页上页下页下页,)max()(00upIAupIAukkk 其中其

43、中.)()(10 niikiikxpaupIA 同理可得下面的定理同理可得下面的定理.上页上页下页下页 定理定理11 设设ARnn有有n个线性无关的特征向量，个线性无关的特征向量， 矩阵矩阵A的特征值及对应的特征向量分别记为的特征值及对应的特征向量分别记为 i 及及xi (i=1,2,n)，而，而p为为 j的近似值，的近似值，(A- -pI)- -1存在，且存在，且 (1),)max(limjjkkxxu (2)jkjkkvppv )max(1,1)max(lim即即则对任意非零初始向量则对任意非零初始向量u0(aj 0) ，由反幂法计算公式，由反幂法计算公式(2.12)构造的向量序列构造的向

44、量序列vk，uk满足满足. )( |jippij . |min/ |pprijij 且收敛速度为且收敛速度为上页上页下页下页 由该定理知，对由该定理知，对A- -pI(其中其中p j)应用反幂法，可应用反幂法，可用来计算特征向量用来计算特征向量xj，只要选择，只要选择p是是 j的一个较好的近的一个较好的近似且特征值分离情况较好，一般似且特征值分离情况较好，一般r很小，常常只要迭很小，常常只要迭代一二次就可完成特征向量的计算代一二次就可完成特征向量的计算.反幂法迭代公式中的反幂法迭代公式中的vk是通过解方程组是通过解方程组.)(1 kkuvpIA求得的求得的, 为了节省工作量为了节省工作量, 可

45、以先将可以先将A- -pI进行三角分解进行三角分解.)(LUpIAP 于是求于是求vk相对于解两个三角形方程组相对于解两个三角形方程组.,1kkkkyUvPuLy 上页上页下页下页实验表明实验表明, 按下述方法选择按下述方法选择u0是较好的是较好的: 选选u0使使)13. 2()1 , 1 , 1(011 PuLUv用回代求解三角形方程组用回代求解三角形方程组(2.13)即得即得v1，然后再按公，然后再按公式式(2.12)进行迭代进行迭代.上页上页下页下页反幂法计算公式：反幂法计算公式：1.分解计算分解计算P(A- -pI)=LU,且保留且保留L, U及及P信息信息2.反幂迭代法反幂迭代法(1

46、) 解解Uv1=(1,1,1)T求求v1.),max(kkkkk/vuv (2) 解解k=2,3, 解解Lyk=Puk- -1求求yk 解解Uvk=yk求求vk k=max(vk) 计算计算uk=vk/ k上页上页下页下页例例5 用反幂法求用反幂法求 410131012A的对应于计算特征值的对应于计算特征值 =1.2679(精确特征值为精确特征值为 )的特征向量的特征向量(用用5位浮点数进行计算位浮点数进行计算).333 解解 用部分选主元的三角分解将用部分选主元的三角分解将A- -pI(其中其中p=1.2679)分解为分解为 P(A- -pI)=LU,其中其中,126807. 07321.

47、0010001 L上页上页下页下页,1029405. 0007321. 21017321. 113 U.001100010 P由由Uv1=(1,1,1)T,得得,)8 .3400, 3 .9290,12692(1Tv ,)26796. 0,73206. 0, 1(1Tu 由由LUv2=Pu1 ,得得,)4 .5467,14937,20404(2Tv ,)26796. 0,73206. 0, 1(2Tu 3对应的特征向量是对应的特征向量是,)26795. 0,73205. 0, 1()32,31, 1(3TTx 由此看出由此看出u2是是x3的相当好的近似的相当好的近似.特征值特征值 3的真值为的

48、真值为,26794901. 1/12679. 123 .26794912. 1333 上页上页下页下页8.3 正交变化与矩阵分解正交变化与矩阵分解正交变化是计算矩阵特征值的有力工具，本节介正交变化是计算矩阵特征值的有力工具，本节介绍豪斯绍豪斯霍尔德霍尔德(Householder)变换变换和和吉文斯吉文斯(Givens)变变换换，并利用它们讨论，并利用它们讨论矩阵分解矩阵分解. 主要讨论主要讨论实矩阵实矩阵和和实实向量向量.上页上页下页下页8.3.1 豪斯霍尔德变换豪斯霍尔德变换定义定义2 设向量设向量w Rn，且，且wTw=1，称矩阵，称矩阵)1 . 3(.212222122221)(2212

49、221212121 nnnnnwwwwwwwwwwwwwwwwHTwwIwH2)( 为为初等反射矩阵初等反射矩阵，这个矩阵也称为，这个矩阵也称为豪斯霍尔德变换豪斯霍尔德变换. 如果记如果记w=(w1,w2, ,wn)，则有，则有上页上页下页下页定理定理12 设有初等反射矩阵，其中设有初等反射矩阵，其中H=I- -2wwT，其，其中中wTw=1，则，则(1) H是对称矩阵，即是对称矩阵，即HT=H.(2) H是正交矩阵，即是正交矩阵，即H- -1=H.(3) 设设A为对称矩阵，那么为对称矩阵，那么 A1=H- -1AH=HAH 即亦即亦是对称矩阵是对称矩阵.证明证明 只证只证H的正交性，其中显然

50、的正交性，其中显然. HTH=H2=(I- -2wwT)(I- -2wwT)=I- -4wwT+4w(wTw)wT=I.设向量设向量u 0，则显然，则显然是一个初等反射矩阵是一个初等反射矩阵.222.TuuHIu上页上页下页下页vwSxv 图图 8-18-1Oy下面考察下面考察初等反射矩阵初等反射矩阵的的几何意义几何意义. 参见参见图图8-1，考，考虑以为虑以为w法向量且过原点法向量且过原点O的超平面的超平面S:wTx=0. 设任意向设任意向量量v Rn, 则则v=x+y, 其中其中x S, y ST. 于是于是 Hx=(I- -2wwT)x=x- -2wwTx=x.对于对于y ST，易知，易

51、知Hy=- -y, 从而对从而对任意向量任意向量v Rn，总有，总有 Hv=x- -y=v .其中其中v 为为v关于平面关于平面S的镜面反射的镜面反射(见图见图8-1).初等反射矩阵在计算上的意义是它能用来约化矩初等反射矩阵在计算上的意义是它能用来约化矩阵，例如设向量阵，例如设向量Hx 0，可选择一初等反射矩阵，可选择一初等反射矩阵H，使，使Hx=y.上页上页下页下页定理定理13 设设x,y为两个不相等的为两个不相等的n维向量维向量, |x|2=|y|2,则存在一个初等反射矩阵则存在一个初等反射矩阵H，使，使Hx=y.2222().TTTxyHIwwIxyxy 2xywxy 证明证明 令令 ,

52、 则得到一个初等反射矩阵则得到一个初等反射矩阵而且而且2222()()2()2.TTTTxyxyx xy xHxxxyxxxyxy 因为因为22() ()2().TTTxyxyxyx xy x 所以所以().Hxxxyy 容易说明，容易说明，w是使是使Hx=y成立的唯一长度等于成立的唯一长度等于1的的向量向量(不计符号不计符号).上页上页下页下页定理定理14(约化定理约化定理) 设设x=(x1,x2,xn)T 0, 则存在则存在初等反射矩阵初等反射矩阵H，使，使Hx=- -e1，其中，其中1121212,sgn(),(3.2),1().2THIuuxxuxeux 证明证明 记记y=- -e1，

53、设，设x y，取，取=|x|2，则有，则有|x|2=|y|2，于是由定理，于是由定理13存在存在H变换变换 H=I- -2wwT ，其中其中 ，使，使Hx=y=- -e1. .112xewxe 上页上页下页下页记记u=x+ +e1=(u1,u2,un)T，于是，于是,(0).xdxxdd 1222.TTuuHIIuuu 其中记其中记u=(x1+ +,x2,xn)T， . 显然显然2212u 如果如果和和x1异号，那么计算异号，那么计算x1+ +时有效数字可能时有效数字可能损失，我们取损失，我们取和和x1有相同的符号，即取有相同的符号，即取1/221121sgn()sgn().niixxxx 在

54、计算在计算时，可能上溢或下溢，为了避免溢出，时，可能上溢或下溢，为了避免溢出，将将x规范化规范化上页上页下页下页则有则有H 使使 H x = e1，其中，其中221221(9, 5,1,1) ,108,54,2Tuxeuu 12(),/,/,/,.THIu ud uu ddHH 例例6 设设x=(3, 5, 1, 1)T，则，则|x|2=6. 取取=6,12745994529551,549553195153THIuu 可直接验证可直接验证Hx=- -e1=(- -6, 0, 0, 0)T.上页上页下页下页8.3.2 吉文斯变换吉文斯变换设向量设向量x,y R2，则变换，则变换1122cossi

55、n,sincosyxyx .yPx 或或是平面上向量的一个旋转变换，其中是平面上向量的一个旋转变换，其中cossin( ),sincosP 为正交变换为正交变换.Rn中变换中变换.yPx 其中其中x=(x1,x2,xn)T, y=(y1,y2,yn)T，而，而上页上页下页下页11cossin1( , , )(3.3)1sincos11iPP i jj 第第 行行第第 行行称为称为Rn中平面中平面xi, xj的的旋转变换旋转变换，也称为，也称为吉文斯变换吉文斯变换. P=P(i, j, )=P(i, j)称为平面旋转矩阵称为平面旋转矩阵.上页上页下页下页显然，显然，P=P(i, j, )具有性质

56、：具有性质：(1) P与单位阵与单位阵I只是在只是在(i, i), (i, j), (j, i), (j, j)位置元位置元素不一样，其它相同素不一样，其它相同.(2) P为正交矩阵为正交矩阵(P- -1=PT).(3) P(i, j)A(左乘左乘)只需计算第只需计算第i行与第行与第j行元素，即行元素，即对对A=(aij)mn有有cs,1,2, .scililjljlaalnaa 其中其中c=cos , s=sin .(4) AP(i, j)(右乘右乘)只需计算第只需计算第i列与第列与第j列元素列元素cs(,)(,),1,2,.scliljliljaaaalm 利用平面旋转变换，可得向量利用平

57、面旋转变换，可得向量x中的指定元素变为零中的指定元素变为零.上页上页下页下页定理定理15(约化定理约化定理) 设设x=(x1,xi,xj,xn)T, 其其中不全为零，则可选择平面旋转阵中不全为零，则可选择平面旋转阵P=P(i, j, )，使，使1(,0,) ,TinPxxxxij 其中其中22,arctan(/).iijjixxxxx 于是，由于是，由c, s的取法得的取法得证明证明 取取c=cos =xi/x i , s=sin =xj/x i ，由，由P(i, j, )x =x =(x 1,x i,x j,x n)T, 利用矩阵乘法利用矩阵乘法, 显然有显然有, .iijjijkkxcxs

58、xxsxcxxxki j 22,0.iijjxxxx上页上页下页下页8.3.3 矩阵的矩阵的QR分解与舒尔分解分解与舒尔分解 定理定理16 设设ARnn非奇异，则存在正交矩阵非奇异，则存在正交矩阵P，使使PA=R，其中，其中R为上三角矩阵为上三角矩阵. 证明证明 我们先用吉文斯变换给出构造我们先用吉文斯变换给出构造P的方法的方法 (1) 第第1步约化，由设有步约化，由设有j(j=1,2, ,n)使使aj1 0，则，则可选择吉文斯变换可选择吉文斯变换P(1, j)，将，将aj1处的元素化为零处的元素化为零. 即即若若aj1 0 (j=2,3, ,n)，则存在，则存在P(1, j)使得使得可简记为

59、可简记为P1A=A(2)，其中，其中P1=P(1,n)P(1,2).)2 , 1(), 1()2()2()2(2)2(2)2(2211211AaaaarrrAPnPnnnnn 上页上页下页下页 (2) 第第k步约化，设上述过程已完成第步约化，设上述过程已完成第1步至第步至第k-1步，于是有步，于是有 由设有由设有j(n j k)使使ajk(k) 0 (j=k+1,n)，则可选择，则可选择吉文斯变换吉文斯变换P(k, j) (j=k+1,n)使使1112112222( )221( )( )( )( ).knknkkkkkkknkknknnrrrrrrrPP P AAaaaa PkA(k)=P(k

60、, n)P(k, k+1)A(k)=PkPk- -1P1A=A(k+1),其中其中 Pk=P(k, n)P(k, k+1) .上页上页下页下页 (3) 继续上述约化过程，最后则有继续上述约化过程，最后则有Pn- -1P2P1A=A(n)=R (上三角形矩阵上三角形矩阵). 令令P=Pn- -1P2P1, 它是一个正交矩阵它是一个正交矩阵, 有有PA=R. 也可以用豪斯霍尔德变换构造出正交矩阵也可以用豪斯霍尔德变换构造出正交矩阵P，记，记A(0)=A，它的第一列记为，它的第一列记为u1，不妨设，不妨设a1(0) 0，可按公，可按公式式(3.2)找到矩阵找到矩阵 使使111111,n nTHRHI