矩阵特征值与特征向量的计算

1、.薅羂莄蒅袄羂肄蚁螀羁膆蒄蚆肀艿虿薂聿莁蒂袁肈肁芅袇肇芃薀螃肆莅莃虿肆肅蕿薅肅膇莁袃肄芀薇蝿膃莂莀蚅膂肂薅薁膁芄莈羀膀莆蚃袆膀蒈蒆螂腿膈蚂蚈螅芀蒄薄螄莃蚀袂袃肂蒃螈袃膅蚈蚄袂莇蒁蚀袁葿莄罿袀腿蕿袅衿芁莂螁袈莄薈蚇袇肃莀薃羇膆薆袁羆芈荿螇羅蒀薄螃羄膀蒇虿羃节蚃薅羂莄蒅袄羂肄蚁螀羁膆蒄蚆肀艿虿薂聿莁蒂袁肈肁芅袇肇芃薀螃肆莅莃虿肆肅蕿薅肅膇莁袃肄芀薇蝿膃莂莀蚅膂肂薅薁膁芄莈羀膀莆蚃袆膀蒈蒆螂腿膈蚂蚈螅芀蒄薄螄莃蚀袂袃肂蒃螈袃膅蚈蚄袂莇蒁蚀袁葿莄罿袀腿蕿袅衿芁莂螁袈莄薈蚇袇肃莀薃羇膆薆袁羆芈荿螇羅蒀薄螃羄膀蒇虿羃节蚃薅羂莄蒅袄羂肄蚁螀羁膆蒄蚆肀艿虿薂聿莁蒂袁肈肁芅袇肇芃薀螃肆莅莃虿肆肅蕿薅肅膇莁袃肄芀

2、薇蝿膃莂莀蚅膂肂薅薁膁芄莈羀膀莆蚃袆膀蒈蒆螂腿膈蚂蚈螅芀蒄薄螄莃蚀袂袃肂蒃螈袃膅蚈蚄袂莇蒁蚀袁葿莄罿袀腿蕿袅衿芁莂螁袈莄薈蚇袇肃莀薃羇膆薆袁羆芈荿螇羅蒀薄螃羄膀蒇虿羃节蚃薅羂莄蒅袄羂肄蚁螀羁膆蒄蚆肀艿虿薂聿莁蒂袁肈肁芅袇肇芃薀螃肆莅莃虿肆肅 第3章 矩阵特征值与特征向量的计算一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。3.1 特征值的估计较粗估计r(A) £ |A|欲将复平面上的特征值一个个用圆盘围起来。3.1.1 盖氏图定义3.1-1 设A = aijn´n，称由不等式所确定的复区域为A的第i个盖氏图，记为Gi，i = 1，2，n。定理3.1

3、-1 若l为A的特征值，则证明：设Ax = lx (x ¹ 0)，若k使得因为ÞÞÞ例1 估计方阵特征值的范围解：G1 = z：|z 1|£ 0.6；G2 = z：|z 3|£ 0.8；G3 = z：|z + 1|£ 1.8；G4 = z：|z + 4|£ 0.6。G1G2G3G4注：定理称A的n个特征值全落在n个盖氏圆上，但未说明每个圆盘内都有一个特征值。3.1.2 盖氏圆的连通部分称相交盖氏圆之并构成的连通部分为连通部分。孤立的盖氏圆本身也为一个连通部分。定理3.1-2 若由A的k个盖氏圆组成的连通部分，含且仅

4、含A的k个特征值。证明: 令D = diag(a11，a12，ann)，M = A D，记则显然有A(1) = A，A(0) = D，易知A(e)的特征多项式的系数是e的多项式，从而A(e)的特征值l1(e)，l2(e)，ln(e)为e的连续函数。A(e)的盖氏圆为：因为A(0) = D的n个特征值a11，a12，ann，恰为A的盖氏圆圆心，当e由0增大到1时，li(e)画出一条以li(0) = aii为始点，li(1) = li为终点的连续曲线，且始终不会越过Gi；aiili 不失一般性，设A开头的k个圆盘是连通的，其并集为S，它与后n k个圆盘严格分离，显然，A(e)的前k个盖氏圆盘与后n

5、 k个圆盘严格分离。 当e = 0时，A(0) = D的前k个特征值刚好落在前k个圆盘G1，Gk中，而另n k个特征值则在区域S之外，e从0变到1时，与始终分离（严格）。连续曲线始终在S中，所以S中有且仅有A的k个特征值。注：1) 每个孤立圆中恰有一个特征值。2) 例1中G2，G4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G1，G3构成的连通部分应含有两个特征值。3) 因为例1中A为实方阵，所以若l为A的特征值，则也是A的特征值，所以G2，G4中各有一个实特征值。3.1.3 盖氏圆与相似变换由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换

6、将盖氏圆的半径变小，以得到更好的估计。原理：取对角阵作相似变换阵：P = diag(b1，b2，bn)其中bi > 0，i = 1，2，n则与A有相同特征值.而B的第i个盖氏圆为：， 适当选取b1，b2，bn就有可能使B的某些盖氏圆的半径比A的相应盖氏圆的半径小。1) 欲缩小Gi，可取bi最大。2) 欲缩小除Gi外的圆，可取bi最小。例2，估计的特征值范围。解：A的三个盖氏圆分别为：G1 = z：| z 0.9|£ 0.13；G2 = z：| z 0.8|£ 0.14；G3 = z：| z 0.4|£ 0.03l3 Î G3，较好。为了更好地估计另

7、外两个特征值可取b3最小：取b1 = b2 = 1，b3 = 0.1即，则所以G1' = z：| z 0.9|£ 0.022；G2' = z：| z 0.8|£ 0.023；G3' = z：| z 0.4|£ 0.3三个盖氏圆分离，故有l1 Î G1'，l2 Î G2'，l3 Î G3。3.2 幂法与反幂法幂法是求方阵的最大特征值及对应特征向量的一种迭代法。3.2.1 幂法设An有n个线性相关的特征向量v1，v2，vn，对应的特征值l1，l2，ln，满足|l1| > |l2| ³

8、 ³ |ln| (3.2.1)1. 基本思想因为v1，v2，vn为Cn的一组基，所以任给x(0) ¹ 0， 线性表示所以有 (3.2-2)若a1 ¹ 0，则因知，当k充分大时 A(k)x(0) » l1ka1v1 = cv1 属l1的特征向量另一方面，记max(x) = xi，其中|xi| = |x|¥，则当k充分大时，若a1 = 0，则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0，迭代下去可求得l1及对应特征向量的近似值。2. 规范化在实际计算中，若|l1| > 1则|l1ka1| ®¥，若|l1| &

9、lt; 1则| l1ka1| ® 0都将停机。须采用“规范化”的方法， k = 0，1，2， (3.2-4)定理3.2-1 任给初始向量有， (3.2-5)证明：而注：若的特征值不满足条件(3.2.1)，幂法收敛性的分析较复杂，但若l1 = l2 = = l r且|l1| > |l r +1| ³ ³ |ln|则定理结论仍成立。此时不同初始向量的迭代向量序列一般趋向于l1的不同特征向量。3. 算法求maxa(x)的流程，设数组x(n)数向量x的n个分量数组x = nk = 1for(i = 2 to n, i+)若|xi| > |xk|Tk = ima

10、x = xk幂法流程：输入数组x0, eps, Ax1 = x0y = x1/maxa(x1)x0 = Ay|maxa(x1) maxa(x0)| < eps输出y, maxa(x0)例1，用幂法求的最大模特征值及对应特征向量见P312function y = maxa(x)k=1;n=length(x);for i=2:n if (abs(x(i)>abs(x(k),k=i; end;end;y=x(k);A=2,4,6;3,9,15;4,16,36;x0=1;1;1;y=x0/maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)>0.001

11、x0=x1; y=x0/maxa(x0) x1=A*yend;ymaxa(x1)3.2.2 加速方法幂法的迭代公式：当k ®¥时，max(x(k) ® l1，其中|l1| > |l2| ³ ³ |ln|注：幂法的收敛速度取决于比值|l2| / |l1|，考虑收敛加速1. 特征值的Aitken加速法(1) 思想：由定理3.2.1的证明知Þ(3.2.6)解之得 (3.2.7)使用l1(k+2)作为l1的近似值的算法称为Aitken加速法。(2) Aitken加速法设xk线性收敛到x*，即存在c，|c| < 1，满足xk+1 x

12、* = (c dk)( xk x*)，其中令则算法：计算流程图输入x0计算max(x0)，y0 = x0/max(x0)计算x1=A y0，max(x1)，y1= x1/max(x1)x2 = A y1，l1 = l0计算max(x2)y2= x2/max(x2)l0=max(x2)-max(x2)-max(x1)2/max(x2)-2max(x1)+max(x0)x0 = x1，x1 = x2|l1 - l0| > eps输出l0例2 用幂法求方阵A的最大模特征值，并用Aitkem加速法解：见（P314）x0=1;1;1;y0=x0/maxa(x0)x1=A*y0;y1=x1/maxa

13、(x1)x2=A*y1;y2=x2/maxa(x2)l0=maxa(x2)-(maxa(x2)-maxa(x1)2/(maxa(x2)-2*maxa(x1)+maxa(x0)while (abs(l1-l0)>0.01 x0=x1;x1=x2;l1=l0; x2=A*y2 maxk=maxa(x2) y2=x2/maxk l0=maxa(x2)-(maxa(x2)-maxa(x1)2/(maxa(x2)-2*maxa(x1)+maxa(x0)end;2. 原点平移法思想：由矩阵论知，若l为A的特征值则l a为A aI的特征值，且特征向量相同。若l1 a为A aI的最大模特征值，且。（lk

14、 a是A aI的次最大模特征值），则对A aI计算l1 a及对应的特征向量比对计算收敛得快，此即为原点平移法。计算l1 a及特征向量的迭代公式特征向量：，max(x(k) ® l1 a，Þ a + max(x(k) ® l1。注：a的选取较为困难。例3 设，求最大模特征值及特征向量。解：（P315）幂法：A=-3,1,0;1,-3,-3;0,-3,4;x0=0;0;1;k=1;y=x0/maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)>0.001 x0=x1; y=x0/maxa(x0) k=k+1 x1=A*yend;ym

15、axa(x1)原点平移法：A=-3,1,0;1,-3,-3;0,-3,4;x0=0;0;1;k=1;y=x0/maxa(x0)x1=(A+4*eye(3)*ywhile(abs(maxa(x1)-maxa(x0)>0.001 x0=x1; y=x0/maxa(x0) k=k+1 x1=(A+4*eye(3)*yend;ymaxa(x1)-43. 对称矩阵的Rayleigh商加速法定义 设A对称，x ¹ 0，则称为关于的Rayleigh商思想：A对称特征值l1，l2，ln均为实数，且存在特征向量v1，v2，vn为标准正交基。设，a1 ¹ 0，则当k充分大时，M'

16、与k无关)注；此比Aitken加速中的(3.2-6)更快公式称为Rayleigh商加速法。其中注：有了R(x(k)，R(x(k+1)，R(x(k+2)，的值，可再用Aitken加速法得到的一个更好的近似值：因为所以例4 设，用Rayleigh商加速法求的最大模特征值及特征向量，并与幂法相比较。解：（P317）幂法：A=6,2,1;2,3,1;1,1,1;x0=1;1;1;k=1;y=x0/maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)>0.001 x0=x1; y=x0/maxa(x0) x1=A*y k=k+1end;ymaxa(x1)Rayleig

17、h商加速法：A=6,2,1;2,3,1;1,1,1;x0=1;1;1;k=1;r=0;y=x0/maxa(x0)x1=A*ywhile(abs(r1-r)>0.001 x0=x1;r1=r; y=x0/maxa(x0) x1=A*y r = y'*x1/(y'*y) k=k+1end;ymaxa(x1)r3.2.3 反幂法用代替作幂法，即反幂法1. 求最小模特征值及相应的特征向量若可逆，|l1| > |l2| ³ ³ |ln|为其特征值，则为A-1的最大模特征值。迭代公式：x(k+1) = A1x(k)，k = 0，1，2，但A1不易求，通常可解

18、方程组Ax(k+1) = x(k)来求x(k+1)即有 (3.2.12)当k ® ¥时有注：为解(3.2-12)中的方程组。对作LR分解（带行交模）PA = LR则有2. 求任一特征值及相应特征向量反幂法结合原点平移法思想：若已知为lj的近似值，则的特征值是而显然非常大（最大）比值很小迭代公式：当k ® ¥时有注：(1) 若有分LR解则迭代公式 (3.2-16)(2) 在(3.2-16)中直接取z(1) = (1，1)T作初值开始迭代称为半次迭代法例5 设的一个特征值的l的近似值，用带原点平移的反幂法求l及相应的特征向量见P320format long;

19、A=-1,2,1;2,-4,1;1,1,-6;x0=1;1;1;B=A+6.42*eye(3);C=lu(B);R = triu(C,0);L =eye(3)+ tril(C,-1);y=x0/maxa(x0);z=1,1,1'x1=inv(R)*zwhile(abs(maxa(x1)-maxa(x0)>0.001 x0=x1; y=x0/maxa(x0) z=inv(L)*y x1=inv(R)*zend;-6.42+1/maxa(x1)预备知识：矩阵论1. 矩阵QR分解定理 设可逆，则存在正交阵Q与上三角阵R使A=QR注：方法 1) 使用史密斯正交变换 2) 使用Househ

20、older变换（反） 3) 使Givens变换2. 矩阵Schur分解定理 设，则存在正交阵Q使实Schur型其中至多2阶。若1阶，其元素即A的特征值若2阶其特征值为A的一对共轭复特征值。注：想加快迭代速度通常先将A化为上Hessenberg阵3. 正交相似于一个n阶上Hessenberg矩阵（）证明：见（P125）§3.3 QR方法 QR方法即使用QR分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛使用的方法之一。3.3.1 QR方法的计算公式思想：从A1 = A出发用正交相似变换得序列Ak使当k ® ¥时，Ak本质收敛到块上三角阵方法：设A1 = Q

21、1R1（QR分解），令A2 = R1Q1，设A2 = Q2R2，令A3 = R2Q2，即 k = 1，2， (3.3-1)Ak的性质： Ak A：Ak+1 = RkQk = (Qk-1Ak)Qk= Rk = Qk-1Ak= (Q1Qk)-1A(Q1Qk)记Gk = Q1Qk 正交，故有Ak A，且有A1Gk = GkAk+1 记Hk = RkR1，则Ak = GkHkQR分解GkHk = (Q1Qk)(RkR1) = Gk-1QkRk Hk-1 = Gk-1AkHk-1= A1Gk-1Hk-1 = = A1k = Ak注：为求得A的特征值，只须Ak能趋于块上三角阵。定义： 矩阵列Ak，当k &

22、#174; ¥时，若其对角元均收敛，且严格下三角部分元素收敛到0，则称Ak本质收敛到上三角阵。 矩阵列Ak，当k ® ¥时，若其对角子块收敛到1阶或2阶的方阵，其下部收敛到0，则称Ak本质收敛到块上三角阵。定理3.3-1设A的特征值满足条件：|l1| > |l2| > > |ln| > 0，vi为li对应的特征向量，i = 1，2，n。记X = (v1，v2，vn)，若有直接三角分解X-1 = LU（杜利特尔分解），则(3.3-1)序列Ak本质收敛于上三角阵，其主对角元素均为的特征值。例1 用QR方法求的A特征值，其中见（P322）注：若A

23、不满足定理条件，Ak不一定本质收敛于上三角矩阵。3.3.2 上Hessenberg矩阵的QR方法及带原点平移的QR方法 在使用QR方法之前，先A将作正交相似变换化为上Hessenberg矩阵H，然后对H作QR迭代，可大量节省运算量。 Givens变换记s = sinq，c = cosq，则为旋转变换正交阵。推广到n维：称为Givens矩阵或Givens变换（旋转变换）。易知J(i，k，q)为正交阵。 对上Hessenberg矩阵用Givens变换作QR分解令hi+1,i* = si hii + ci hi+1,i = 0，即选择qi使右边第i+1行第i列元素为0，而H的第i行与第i+1行零元素位置上左乘J(i，i+1，qi)后仍为0，其他行则不变。（可以证明） 这样i = 1，2，n-1共n-1次左乘J后H变为上三角阵R。即U定理 = R，其中UT= J(n-1，n，qn-1)J(1，2，q1)正交，且为下Hessenberg阵，Þ U为上Hessenberg阵 Þ H = UR (QR分解) 记H1 = H，设H1= U1R1 令H