第8章矩阵特征值计算

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

2、特征向量的问题.上页上页上页上页上页上页下页下页下页下页下页下页8.1.1 特征值问题及性质特征值问题及性质求求A的特征值问题等价于求的特征值问题等价于求A的特征方程的特征方程设矩阵设矩阵A Rnn，特征值问题是求，特征值问题是求 C和非零向和非零向量量x Rn，使，使( )det()0pIA的根的根. Ax= x 其中其中 是矩阵是矩阵A的的特征值特征值，x是矩阵是矩阵A属于特征值属于特征值 的的特征向量特征向量. 上页上页上页上页上页上页下页下页下页下页下页下页 定理定理1 设设 为为ARnn的特征值的特征值, 且且Ax= x, x 0，则有则有 (2) - - 为为A- - I的特征值的

3、特征值, 且且 (A- - I)x=( - - )x ; (1) c 为的为的cA特征值特征值(c 0为常数为常数), 且且 (cA)x=(c )x ; (3) k为为Ak的特征值的特征值, 且且Akx= kx ； (4) 若若A为非奇异矩阵为非奇异矩阵, 则则 - -1为为A- -1的特征值的特征值, 且且 A- -1x= - -1x . 特征值的基本性质特征值的基本性质.上页上页上页上页上页上页下页下页下页下页下页下页 定理定理2 (1) 设矩阵设矩阵ARnn可对角化可对角化, 即存在非奇即存在非奇异矩阵异矩阵P 使使的的充分必要条件充分必要条件是是A具有具有n个线性无关的特征向量个线性无

4、关的特征向量. (2) 如果矩阵如果矩阵ARnn有有m个个(mn)不同的特征不同的特征值值 1, 2, m，则对应的特征向量，则对应的特征向量 x1, x2, , xm 线线性无关性无关.121,nPAP 上页上页上页上页上页上页下页下页下页下页下页下页 定理定理3 设设ARnn为对称矩阵为对称矩阵, 则则 (3) 存在一个正交矩阵存在一个正交矩阵P 使的使的且且 1, 2, n为为A的特征值的特征值, 而而P(u1,u2,un)的列向的列向量量uj为为A的对应于的对应于 j 的单位特征向量的单位特征向量.121,TnPAPP AP (1) A的特征值均为实数的特征值均为实数; (2) A有有

5、n个线性无关的特征向量个线性无关的特征向量;上页上页上页上页上页上页下页下页下页下页下页下页 定义定义 设设ARnn为对称矩阵为对称矩阵, 对于任一非零向量对于任一非零向量x 0, 称称,),(),()(xxxAxxR 为对应于向量为对应于向量x的的瑞利瑞利(Rayleigh)商商. 定理定理4 设设ARnn为对称矩阵为对称矩阵(其特征值依次记其特征值依次记为为 1 2 n), 则则 (1). (对任何非零对任何非零xRn);1(, )( , )nAx xx x (2). ;10(, )max( , )nx RxAx xx x 0(, )min.( , )nnx RxAx xx x 上页上页上

6、页上页上页上页下页下页下页下页下页下页 证明证明 只证只证(1). 由于由于A为实对称矩阵为实对称矩阵,可将可将 1, 2, , n 对应的特征向量对应的特征向量 x1, x2, , xn 正交规范化正交规范化, 即有即有 (xi, xj)=ij , 设设x 0为为Rn中任一向量中任一向量, 则则. 0,1221 niiniiixxx 于是于是2121(, )( , )niiiniiAx xx x 瑞利商瑞利商必位于必位于 n和和 1之间之间, (1)成立成立.n 1. 上页上页上页上页上页上页下页下页下页下页下页下页8.1.2 特征值估计与扰动特征值估计与扰动 定义定义1 设设n阶矩阵阶矩阵

7、A=(aij)，令，令).,(niarnijjiji211 (1) ; (2) 集合集合 称为复平面上以称为复平面上以aii为圆心，以为圆心，以ri为半径的为半径的n阶矩阵阶矩阵A的的n个个格什戈林格什戈林(Gerschgorin)圆盘圆盘. |,(1,2, )iiiiDzzar zCin上页上页上页上页上页上页下页下页下页下页下页下页 定理定理5 (Gerschgorin圆盘定理圆盘定理) 特别地，如果特别地，如果A的一个圆盘的一个圆盘Di是与其它圆盘分离是与其它圆盘分离(即即孤立圆盘孤立圆盘)，则，则Di中精确地包含中精确地包含A的一个特征值的一个特征值.1(1,2,. ).niiiijj

8、j iarain (1) 设设n阶矩阵阶矩阵A(aij)，则，则A的每一个特征值必属的每一个特征值必属于下面某个圆盘之中于下面某个圆盘之中 (2) 如果如果A有有m个个圆盘圆盘组成一个连通的并集组成一个连通的并集S，且，且S与余下与余下n- -m个个圆盘圆盘是分离的，则是分离的，则S内恰包含内恰包含A的的m个个特征值特征值.或者说或者说 A的特征值都在的特征值都在n个圆盘的个圆盘的并集并集中中.上页上页上页上页上页上页下页下页下页下页下页下页 证明证明 只就只就(1)给出证明给出证明. .1knkjjkjkkraa Ax= x，其中，其中 x=(x1,x2, xn)T 0.或或记记 , 考虑考

9、虑Ax= x的第的第k个方程个方程, 即即1maxkii nxxx ,1knjjkjxxa ,)( kjjkjkkkxaxa 于是于是即即,kkkkjjkkjj kj kaxaxxa 设设 为为A的特征值的特征值, 上页上页上页上页上页上页下页下页下页下页下页下页A的每一个特征值必位于的每一个特征值必位于A的一个圆盘中的一个圆盘中, 并且相并且相应的特征值应的特征值 一定位于第一定位于第k个圆盘中个圆盘中(其中其中k是对应特征是对应特征向量向量x绝对值最大的分量的下标绝对值最大的分量的下标).利用相似矩阵性质利用相似矩阵性质,有时可以获得有时可以获得A的特征值进一的特征值进一步的估计步的估计,

10、 即适当选取非奇异对角阵即适当选取非奇异对角阵,112111 nD 并做相似变换并做相似变换 .适当选取适当选取 可使某些圆盘半径及连通性发生变化可使某些圆盘半径及连通性发生变化.nnijijaADD 1), 2 , 1(nii 上页上页上页上页上页上页下页下页下页下页下页下页410101114A 例例1 估计矩阵估计矩阵 的特征值范围的特征值范围.解解 矩阵矩阵A的的3个圆盘为个圆盘为. 24:, 2:, 14:321 DDD由圆盘定理可知由圆盘定理可知, A 的的3个特征值位于个特征值位于3个圆盘的并集个圆盘的并集中中, 由于由于D1是孤立圆盘是孤立圆盘, 所以所以D1内恰好包含内恰好包含

11、A的一个的一个特征值特征值 1 (为实特征值为实特征值), 即即. 531 A的其它两个特征值的其它两个特征值 2, 3包含在包含在D2, D3的并集中的并集中.上页上页上页上页上页上页下页下页下页下页下页下页取对角阵取对角阵1100010000.9D 做相似变换做相似变换.49 . 09 . 00101491011 ADDAA矩阵矩阵A1的的3个圆盘为个圆盘为. 8 . 14:,919:, 14:321 EEE3个圆盘都是孤立圆盘个圆盘都是孤立圆盘, 所以所以, 每一个圆盘都包含每一个圆盘都包含A的的一个特征值一个特征值(为实特征值为实特征值), 且有估计且有估计上页上页上页上页上页上页下页

12、下页下页下页下页下页 . 2 . 28 . 5,919919, 53321 定理定理6(Bauer-Fike定理定理) 设设 是是A+IRnn的一个的一个特征值特征值, 且且P- -1AP=D=diag( 1, 2, n), 则有则有其中其中|p为矩阵的为矩阵的p范数，范数，p=1,2, .1()minpppAPPI cond (),pP . 8 . 14:,919:, 14:321 EEE上页上页上页上页上页上页下页下页下页下页下页下页1()minpppAPPI 证明证明 (A)时显然成立时显然成立,故只考虑故只考虑 (A).这时这时D- - I非奇异，设非奇异，设x是是A+I对应于对应于

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

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

15、页下页下页下页下页 计算矩阵计算矩阵A的特征值的特征值. 当当n=2, 3时时, 可按行列式展可按行列式展开的办法求特征方程开的办法求特征方程p( )=0的根的根. 但当但当n较大时较大时, 如果如果按展开行列式的办法按展开行列式的办法, 首先求出首先求出p( )的系数的系数, 再求再求p( )的根的根, 工作量就非常大工作量就非常大,用这种办法求矩阵特征值是不用这种办法求矩阵特征值是不切实际的切实际的,由此需要研究由此需要研究A的特征值及特征向量的数值的特征值及特征向量的数值方法方法. 介绍介绍(计算机上计算机上)两类方法两类方法:1. 幂法及反幂法幂法及反幂法(迭代法迭代法)，2. 正交相

16、似变换的方法正交相似变换的方法(变化法变化法).上页上页上页上页上页上页下页下页下页下页下页下页8.2 幂法及反幂法幂法及反幂法 幂法与反幂法都是求幂法与反幂法都是求实矩阵实矩阵的特征值和特征向量的特征值和特征向量的的向量迭代法向量迭代法. . 幂法幂法是计算矩阵的是计算矩阵的按模最大的特征值按模最大的特征值和和相应特征相应特征向量向量的一种向量迭代法的一种向量迭代法, ,特别适用于大型稀疏矩阵特别适用于大型稀疏矩阵. . 反幂法反幂法是计算非奇异是计算非奇异( (可逆可逆) )矩阵矩阵按模最小的特征按模最小的特征值值和和相应特征向量相应特征向量的一种向量迭代法的一种向量迭代法, ,特别是计算

17、海特别是计算海森伯格阵或三对角阵的对应一个给定近似特征值的森伯格阵或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一特征向量的有效方法之一. . 上页上页上页上页上页上页下页下页下页下页下页下页8.2.1 幂法幂法(又称又称乘幂法乘幂法)现讨论求现讨论求 1及及x1的方法的方法.), 2 , 1(nixAxiii 设实矩阵设实矩阵A=(aij)有一个有一个完全的特征向量组完全的特征向量组, 即即A有有n个线性无关的特征向量个线性无关的特征向量, 设矩阵设矩阵A的特征值为的特征值为 1, 2, n, 相应的特征向量为相应的特征向量为x1,x2,xn. 已知已知A的主特的主特征值征值 1是

18、实根是实根, 且满足条件且满足条件12| |,n 上页上页上页上页上页上页下页下页下页下页下页下页011221(0),nnva xa xa xa设设 幂法的幂法的基本思想基本思想是是: 任取任取非零非零的初始向量的初始向量v0 , 由矩由矩阵阵A构造一向量序列构造一向量序列vk称为迭代向量，由假设称为迭代向量，由假设, v0可唯一表示为可唯一表示为102210110,.,.kkkvAvvAvA vvAvAv 上页上页上页上页上页上页下页下页下页下页下页下页011221(0),nnva xa xa xa设设于是于是1011122211111121()().kkkkkknnnnkkkiiikivA

19、vA vaxaxaxa xaxa x 其中其中21().nkikiiiax 由假设由假设 故故 从而从而), 3 , 2( 1/1nii , 0lim kk 111limkkkva x 为为 1的特征向量的特征向量.上页上页上页上页上页上页下页下页下页下页下页下页所以当所以当k充分大时充分大时, 有有111,kkva x 则则vk为矩阵为矩阵A的的对应特征值对应特征值 1 的一个近似特征向量的一个近似特征向量.用用(vk)i 表示表示vk的第的第i个分量个分量, 则当则当k充分大时充分大时, 有有 11.kikivv 即为即为A的的主特征值主特征值 1的近似值的近似值.111111,kkkkv

20、Ava xv 由于由于 这种由已知非零向量这种由已知非零向量v0及矩阵及矩阵A的乘幂的乘幂Ak构造向量序列构造向量序列vk以计算以计算A的的主特征主特征值值 1及相应特征向量的方法就称为及相应特征向量的方法就称为(乘乘)幂法幂法.上页上页上页上页上页上页下页下页下页下页下页下页 幂法幂法的思想的思想: 由矩阵由矩阵A的乘幂的乘幂 Ak与非零向量与非零向量v0相相乘来构造向量序列乘来构造向量序列vk=Akv0，从而计算主特征值，从而计算主特征值 1及其对应的特征向量及其对应的特征向量. ).(11 kvvikik 的收敛速度由比值的收敛速度由比值,12 r来确定，来确定，r 越小收敛越快，但当越

21、小收敛越快，但当r 1 1时收敛可能很慢时收敛可能很慢.上页上页上页上页上页上页下页下页下页下页下页下页 定理定理7 设设ARnn有有n个线性无关的特征向量个线性无关的特征向量,主主特征值特征值 1满足满足 | 1| 2| | n|, 则对任何非零向量则对任何非零向量v0 ,又设又设A有有n个线性无关的特征向量个线性无关的特征向量, 1对应的对应的r个线性个线性无关的特征向量为无关的特征向量为x1, x2, , xr , 则则如果如果A的主特征值为实的的主特征值为实的r重根重根, 即即 1= 2= r , 且且 | r| r+1| | n|，01111(),rnkkkikiiiiii rvA

22、va xax 111,kkva x 11.kikivv 上页上页上页上页上页上页下页下页下页下页下页下页01111(),rnkkkikiiiiii rvA va xax ).0(lim111 riiiriiikkkxaxav设设 为为 1的特征向量的特征向量, 这说明当这说明当A的主特征值是实的重根时的主特征值是实的重根时,定理定理7的结论还是正确的的结论还是正确的.上页上页上页上页上页上页下页下页下页下页下页下页 应用应用幂法幂法计算计算A的主特征值的主特征值 1及其对应的特征向及其对应的特征向量时量时, 如果如果| 1|1(或或| 1| 2| | n|, 则对任意非零初始向则对任意非零初始

23、向量量v0=u0(a1 0), 有幂法计算公式为有幂法计算公式为(uk,vk)则有则有,)max(lim).1(11xxukk .lim).2(1 kk上页上页上页上页上页上页下页下页下页下页下页下页例例2 用幂法计算矩阵用幂法计算矩阵 0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10 . 1A的主特征值和相应的特征向量的主特征值和相应的特征向量.解解取取 v0=u0=(1,1,1)T, 则则 TTvuv1 ,8182. 0,9091. 01,75. 2,75. 2 ,25. 2 , 5 . 211111 ), 2 , 1(max1 k/vuvAuvkkkkk

24、kk 上页上页上页上页上页上页下页下页下页下页下页下页计算结果如下表计算结果如下表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, 0.6497, 1)2.536545619(0.7482, 0.6497, 1)2.536537420(0.748

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

26、幂法计算矩阵 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.73604 5.46

27、21, 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.

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

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

30、下页下页下页下页下页例例3 设设AR44有特征值有特征值),4 , 3 , 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上页上页上页上页上页上页下页下页下

31、页下页下页下页 下面考虑当下面考虑当A的特征值是实数时，怎样选择的特征值是实数时，怎样选择p使采使采用幂法计算用幂法计算 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=(

32、2+ n)/2P*时时 为为最小值，这时最小值，这时收敛速度的比值收敛速度的比值为为.2212112nnnpppp 当希望计算当希望计算 n时，应选取时，应选取 p=( 1+ n-1)/2P* 使得应用幂法计算使得应用幂法计算 n得到加速得到加速. 当当A的特征值都是实数，满足的特征值都是实数，满足(2.10)式式且且 2, n能初步估计出来，我们就能确定能初步估计出来，我们就能确定P*的近似值的近似值.nn 121上页上页上页上页上页上页下页下页下页下页下页下页 例例4 用原点平移加速法求用原点平移加速法求例例2中矩阵中矩阵A的主特征值的主特征值与其对应的特征向量与其对应的特征向量.25.

33、125. 05 . 025. 025. 00 . 15 . 00 . 125. 0 B,0 . 225. 05 . 025. 00 . 10 . 15 . 00 . 10 . 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.00000

34、005(0.7516, 0.6522, 1)1.79140116(0.7491, 0.6511, 1)1.78884437(0.7488, 0.6501, 1)1.78733008(0.7484, 0.6499, 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.上页上页上页上页上页上页下页下页下

35、页下页下页下页 原点位移的加速方法，是一个矩阵变换方法原点位移的加速方法，是一个矩阵变换方法. 这这种变换容易计算，又不破坏矩阵种变换容易计算，又不破坏矩阵A的稀疏性，但的稀疏性，但p的的选择选择依赖依赖对对A的的特征值分布的大致了特征值分布的大致了解解.及相应的特征向量为及相应的特征向量为 x1 (0.7482, 0.6497, 1)T. 与例与例2结果比较，上述结果比例结果比较，上述结果比例2迭代迭代15次还好次还好. 若若迭代迭代15次次 1 1.7865258(相应的相应的 1 2.5365258)就是真就是真值值.上页上页上页上页上页上页下页下页下页下页下页下页 设设ARnn为为对称

36、矩阵对称矩阵，称，称.),(),()(xxxAxxR 为向量为向量x的的瑞利商瑞利商，其中，其中(x, x)=xTx为内积为内积. 由定理由定理4知道，实对称矩阵知道，实对称矩阵A的特征值的特征值 1及及 n可用可用瑞利商瑞利商的的极限值表示极限值表示. 下面我们将下面我们将瑞利商瑞利商应用到用幂法计算应用到用幂法计算实对称矩阵实对称矩阵A的主特征值的加速上来的主特征值的加速上来.2、瑞利商、瑞利商(Rayleigh)加速加速上页上页上页上页上页上页下页下页下页下页下页下页 定理定理9 设设ARnn为为对称矩阵对称矩阵，特征值满足，特征值满足对应的特征向量对应的特征向量xi满足满足(xi, x

37、j)= ij (单位正交向量单位正交向量), 应应用幂法公式用幂法公式(2.9)计算计算A的主特征值的主特征值 1，则规范化向，则规范化向量量uk的的瑞利商瑞利商给出给出 1的较好的近似值为的较好的近似值为,121nn kkkkkkuuuAuuR2121, 由此可见，由此可见，R(uk)比比 k更快的收敛于更快的收敛于 1.上页上页上页上页上页上页下页下页下页下页下页下页 证明证明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(2121122112200001 knjkjjnjkjjkkkkkkkkk

38、aauAuAuAuAuuuAuuR 得得上页上页上页上页上页上页下页下页下页下页下页下页 幂法的幂法的瑞利商加速迭代公式瑞利商加速迭代公式可以写为可以写为 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A为为n阶实对称矩阵阶实对称矩阵.,11kkux 对给定的误差限对给定的误差限 ，当，当| k k- -1| 时，取近似值时，取近似值上页上页上页上页上页上页下页下页下页下页下页下页8.2.3 反幂法反幂法 反幂法是用于求非奇异矩阵反幂法是用于求非奇异矩阵A的的按模最小的特征按模最小的特征值和对应特征向量值和对应特征向量的方法的方法. 而结合原点平移法

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

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

41、k 求求出出解解求求出出解解), 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)max(limnkkv 收敛速度收敛速度的比值的比值r=| n/ n- -1|. 上页上页上页上页上页上页下

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

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

44、 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)构造的

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

46、得的求得的, 为了节省工作量为了节省工作量, 可以先将可以先将A- -pI进行三角分解进行三角分解.)(LUpIAP 于是求于是求vk相对于解两个三角形方程组相对于解两个三角形方程组.,1kkkkyUvPuLy 上页上页上页上页上页上页下页下页下页下页下页下页实验表明实验表明, 按下述方法选择按下述方法选择u0是较好的是较好的: 选选u0使使)13. 2()1 , 1 , 1(011 PuLUv用回代求解三角形方程组用回代求解三角形方程组(2.13)即得即得v1，然后再按公，然后再按公式式(2.12)进行迭代进行迭代.上页上页上页上页上页上页下页下页下页下页下页下页反幂法计算公式：反幂法计算公

47、式：1.分解计算分解计算P(A- -pI)=LU,且保留且保留L, U及及P信息信息2.反幂迭代法反幂迭代法(1) 解解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 解解 用部分选主元的三角分解将

48、用部分选主元的三角分解将A- -pI(其中其中p=1.2679)分解为分解为 P(A- -pI)=LU,其中其中,126807. 07321. 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对应的特征向量是

49、对应的特征向量是,)26795. 0,73205. 0, 1()32,31, 1(3TTx 由此看出由此看出u2是是x3的相当好的近似的相当好的近似.特征值特征值 3的真值为的真值为,26794901. 1/12679. 123 .26794912. 1333 上页上页上页上页上页上页下页下页下页下页下页下页8.3 正交变换与矩阵分解正交变换与矩阵分解8.3.1 豪斯霍尔德变换豪斯霍尔德变换定义定义 设向量设向量w Rn，且，且wTw=1，称矩阵，称矩阵211212212221212222122( ).2212nnnnnww ww ww www wH ww ww ww TwwIwH2)( 为为

50、初等反射矩阵初等反射矩阵，这个矩阵也称为，这个矩阵也称为豪斯霍尔德变换豪斯霍尔德变换. 如果记如果记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的正交性，其中显然的正交性，其中显然. HTH=H2=(I- -2wwT)

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

52、.初等反射矩阵在计算上的意义是它能用来约化矩初等反射矩阵在计算上的意义是它能用来约化矩阵阵. 设向量设向量Hx 0, 可选择一初等反射矩阵可选择一初等反射矩阵H, 使使Hx=y. Hy=(I- -2wwT)y=y- -2wwTy=-y.上页上页上页上页上页上页下页下页下页下页下页下页定理定理13 设设x,y为两个不相等的为两个不相等的n维向量维向量, |x|2=|y|2,则存在一个初等反射矩阵则存在一个初等反射矩阵H，使，使Hx=y.2222().TTTxyHIwwIxyxy 2xywxy 证明证明 令令 , 则得到一个初等反射矩阵则得到一个初等反射矩阵而且而且2222()()2()2.TTT

53、Txyxy x xy xHxxxyxxxyxy 因为因为22() ()2().TTTxyxyxyx xy x所以所以().Hxxxyyw是使是使Hx=y成立的唯一长度等于成立的唯一长度等于1的向量的向量(不计符号不计符号).上页上页上页上页上页上页下页下页下页下页下页下页定理定理14(约化定理约化定理) 设设x=(x1,x2,xn)T 0, 则存在则存在初等反射矩阵初等反射矩阵H，使，使Hx=- -e1，其中，其中1121212,sgn(),1().2THIuuxxuxeux 证明证明 记记y=- -e1, 设设x y, 则有则有|x|2=|y|2, 由定理由定理13,存在变换存在变换 H=I

54、- -2wwT , 使使Hx=y=- -e1, ,其中其中112.xewxe 上页上页上页上页上页上页下页下页下页下页下页下页,(0).xdxxdd 1222.TTuuHIIuuu 其中记其中记u=(x1+ +,x2,xn)T， . 2212u 如果如果和和x1异号异号, 那么计算那么计算x1+ +时有效数字可能损时有效数字可能损失失, 我们取我们取和和x1有相同的符号有相同的符号, 即取即取1/221121sgn()sgn().niixxxx 在计算在计算时时, 可能上溢或下溢可能上溢或下溢, 为了避免溢出为了避免溢出, 将将x规范化规范化记记u=x+ +e1=(u1,u2,un)T，于是，

55、于是上页上页上页上页上页上页下页下页下页下页下页下页则有则有H 使使 H x = e1，其中，其中221221(9, 5,1,1) ,108,54,2Tuxeuu12(),/,/,/,.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，则变换，则变换1122cossin,sincosyxyx .y

56、Px 或或是平面上向量的一个旋转变换，其中是平面上向量的一个旋转变换，其中cossin( ),sincosP 为正交变换为正交变换.R2中变换中变换.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)称为平面旋转矩阵称为平面旋转矩阵.上页上页上页上页上页上页下页下页下页下页下页下页显然，显

57、然，P=P(i, j, )具有性质：具有性质：(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

58、 利用平面旋转变换，可得向量利用平面旋转变换，可得向量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, 利用

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

60、零. 即即若若aj1 0 (j=2,3, ,n)，则存在，则存在P(1, j)使得使得可简记为可简记为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( )