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

数值计算方法矩阵特征值与特征向量的计算一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。21特征值的估计较粗估计(A)||A||欲将复平面上的特征值一个个用圆盘围起来。1.1盖氏圆定义1-1设A=[aij]nn，称由不等式所确定的复区域为A的第i个盖氏圆，记为Gi：i=1，2，…，n。定理3.1-1若为A的特征值，则3定理1-1若为A的特征值，则证明：设Ax=x(x

0)，若k使得因为4例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}。注：定理称A的n个特征值全落在n个盖氏圆上，但未说明每个圆盘内都有一个特征值。G1G2G3G451.2盖氏圆的连通部分称相交盖氏圆之并构成的连通部分为连通部分。孤立的盖氏圆本身也为一个连通部分。定理3.1-2若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。6定理1-2若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。证明:令D=diag(a11，a22，…，ann)，M=A–D，记则显然有A(1)=A，A(0)=D，易知A()的特征多项式的系数是的多项式，从而A()的特征值1()，2()，…，n()为的连续函数。7计算max(x0)，y0=x0/max(x0)求最小模特征值及相应的特征向量幂法是求方阵的最大特征值及对应特征向量的一种迭代法。G4={z：|z+4|0.1)欲缩小Gi，可取bi最大。注：有了R(x(k))，R(x(k+1))，R(x(k+2))，的值，可再用Aitken加速法得到的一个更好的近似值：(1)思想：由定理3.须采用“规范化”的方法xk+1–x*=(c–k)(xk–x*)，其中3盖氏圆与相似变换迭代公式：x(k+1)=A–1x(k)，k=0，1，2，…，0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]若|1|>1则|1ka1|，若|1|<1则|1ka1|0A=[-3,1,0;1,-3,-3;0,-3,4];=…=(Q1…Qk)-1A(Q1…Qk)输入数组x0,eps,A(2)Aitken加速法记Gk=Q1…Qk——正交，if(abs(x(i))>abs(x(k))),k=i;1-2若由A的k个盖氏圆组成的连通部分，含且仅含A的k个特征值。A()的盖氏圆为：因为A(0)=D的n个特征值a11，a22，…，ann，恰为A的盖氏圆圆心，当由0增大到1时，i()画出一条以i(0)=aii为始点，i(1)=i为终点的连续曲线，且始终不会越过Gi；aiii8不失一般性，设A开头的k个圆盘是连通的，其并集为S，它与后n–k个圆盘严格分离，显然，A()的前k个盖氏圆盘与后n–k个圆盘严格分离。当=0时，A(0)=D的前k个特征值刚好落在前k个圆盘G1，…，Gk中，而另n–k个特征值则在区域S之外，从0变到1时，与始终分离（严格）。连续曲线始终在S中，所以S中有且仅有A的k个特征值。9注：1)每个孤立圆中恰有一个特征值。2)例1中G2，G4为仅由一个盖氏圆构成的连通部分，故它们各有一个特征值，而G1，G3构成的连通部分应含有两个特征值。3)因为例1中A为实方阵，所以若为A的特征值，则也是A的特征值，所以G2，G4中各有一个实特征值。101.3盖氏圆与相似变换由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。原理：取对角阵作相似变换阵：P=diag(b1，b2，…，bn)其中bi>0，i=1，2，…，n则与A有相同特征值.而B的第i个盖氏圆为：，

11而B的第i个盖氏圆为：，

适当选取b1，b2，…，bn就有可能使B的某些盖氏圆的半径比A的相应盖氏圆的半径小。121)欲缩小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.03}3

G3，较好。13为了更好地估计另外两个特征值可取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}三个盖氏圆分离，故有1

G1'，2

G2'，3

G3。14幂法是求方阵的最大特征值及对应特征向量的一种迭代法。2.1幂法设An有n个线性无关的特征向量v1，v2，…，vn，对应的特征值1，2，…，n，满足|1|>|2|…|n|(3.2-1)2幂法与反幂法151.基本思想因为{v1，v2，…，vn}为Cn的一组基，所以任给x(0)

0，——线性表示所以有(3.2-2)若a1

0，则因知，当k充分大时A(k)x(0)

1ka1v1=cv1(属1的特征向量)另一方面，记max(x)=xi，其中|xi|=||x||，则当k充分大时，16若a1=0，则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0，迭代下去可求得1及对应特征向量的近似值。2.规范化在实际计算中，若|1|>1则|1ka1|，若|1|<1则|1ka1|0都将停机。须采用“规范化”的方法，k=0，1，2，…(3.2-4)17规范化：，k=0，1，2，…(3.2-4)定理3.2-1任给初始向量x(0)0，有(3.2-5)18易知A()的特征多项式的系数是的多项式，A=[-3,1,0;1,-3,-3;0,-3,4];解：Matlab代码如下注：想加快迭代速度通常先将A化为上Hessenberg阵注：1)每个孤立圆中恰有一个特征值。G4={z：|z+4|0.欲将复平面上的特征值一个个用圆盘围起来。求任一特征值及相应特征向量while(abs(maxa(x1)-maxa(x0)))>0.定理设ARnn，则存在正交阵Q使迭代公式：x(k+1)=A–1x(k)，k=0，1，2，…，因为{v1，v2，…，vn}为Cn的一组基，所以而B的第i个盖氏圆为：，矩阵特征值与特征向量的计算解：Matlab代码如下输入数组x0,eps,A②矩阵列{Ak}，当k时，若其对角子块收敛到1阶或2阶的方阵，其下部收敛到0，则称{Ak}本质收敛到块上三角阵。A=[6,2,1;2,3,1;1,1,1];=…=(Q1…Qk)-1A(Q1…Qk)注：若A的特征值不满足条件(3.|1-0|>eps证明：19而注：若A的特征值不满足条件(3.2-1)，幂法收敛性的分析较复杂，但若1=2=…=r，且|1|>|r+1|…|n|，则定理结论仍成立。此时不同初始向量的迭代向量序列一般趋向于1的不同特征向量。203.算法求maxa(x)的流程，设数组x(n)存放向量x的n个分量幂法流程：数组x=[n]k=1for(i=2ton,i++)若|x[i]|>|x[k]|Tk=imax=x[k]输入数组x0,eps,Ax1=x0y=x1/maxa(x1)x0=Ay|maxa(x1)–maxa(x0)|>eps输出y,maxa(x0)21Matlab幂法流程：输入数组x0,eps,Ay=x0/maxa(x0)x1=Ay|maxa(x1)–maxa(x0)|>epsx0=x1y=x0/maxa(x0)x1=Ay输出y,maxa(x1)22例1，用幂法求的最大模特征值及对应特征向量（见P312）解：首先给出函数代码：functiony=maxa(x)k=1;n=length(x);fori=2:nif(abs(x(i))>abs(x(k))),k=i;end;end;y=x(k);23幂法代码：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.001x0=x1;y=x0/maxa(x0)x1=A*yend;ymaxa(x1)输入数组x0,eps,Ay=x0/maxa(x0)x1=Ay|maxa(x1)–maxa(x0)|>epsx0=x1y=x0/maxa(x0)x1=Ay输出y,maxa(x1)242.2加速方法幂法的迭代公式：当k

时，

max(x(k))

1，其中|1|>|2|…|n|注：幂法的收敛速度取决于比值|2|/|1|，考虑收敛加速251.特征值的Aitken加速法(1)思想：由定理3.2-1的证明知26(3.2-6)27解之得(3.2-7)使用1(k+2)作为1的近似值的算法称为Aitken加速法。28为了更好地估计另外两个特征值可取b3最小：(1)思想：由定理3.k=1;n=length(x);k=1，2，…计算max(x0)，y0=x0/max(x0)2-16)中直接取z(1)=(1，…，1)T作初值开始迭代称为半次迭代法记Gk=Q1…Qk——正交，适当选取b1，b2，…，bn就有可能使B的某些盖氏圆的半径比A的相应盖氏圆的半径小。y=x0/maxa(x0)一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。y=x1/maxa(x1)l0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0))注：(1)若有LR分解，则迭代公式|1|>|2|…|n|(3.方法2)使用Householder变换（反射）y=x0/maxa(x0)对作LR分解（带行交换）PA=LR则有若a1=0，则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0，迭代下去可求得1及对应特征向量的近似值。x0=x1;k=k+1|1-0|>eps2-12)中的方程组。若|1|>1则|1ka1|，若|1|<1则|1ka1|03盖氏圆与相似变换注:此比Aitken加速中的(3.2上Hessenberg矩阵的QR方法及带原点平移的QR方法由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。解：A的三个盖氏圆分别为：y=x0/maxa(x0)x0=[1;1;1];k=1设{xk}线性收敛到x*，即存在c，|c|<1，满足注：幂法的收敛速度取决于比值|2|/|1|，考虑收敛加速得上Hessenberg阵列{Hk}。输入数组x0,eps,A取平移量对H1–I作QR分解：H1–I=Q1R1，令H2=R1Q1+I，G4={z：|z+4|0.令hi+1,i*=sihii+cihi+1,i=0，即选择i使右边第i+1行第i列元素为0，方法2)使用Householder变换（反射）x1=inv(R)*z即UTH=R，UT=J(n-1，n，n-1)…J(1，2，1)，其中UT正交，且为下Hessenberg阵记X=(v1，v2，…，vn)，若有直接三角分解X-1=LU（杜利特尔分解），则(3.(2)Aitken加速法设{xk}线性收敛到x*，即存在c，|c|<1，满足xk+1–x*=(c–k)(xk–x*)，其中令则步骤：计算29流程图输入x0计算max(x0)，y0=x0/max(x0)计算x1=Ay0，max(x1)，y1=x1/max(x1)x2=Ay1，1=0计算max(x2)y2=x2/max(x2)0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]x0=x1，x1=x2|1-0|>eps输出030Matlab流程图输入A，x0计算max(x0)，y0=x0/max(x0)计算x1=Ay0，max(x1)，y1=x1/max(x1)计算x2=Ay1，max(x2)，y2=x2/max(x2)0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]|1-0|>epsx0=x1，x1=x2，1=0x2=Ay1计算max(x2)y2=x2/max(x2)0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]输出031例2用幂法求方阵A的最大模特征值，并用Aitkem加速法解：见（P314）32幂法A=[-4,14,0;-5,13,0;-1,0,2];x0=[1;1;1];k=1y=x0/maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)))>0.01x0=x1;k=k+1maxa(x0)y=x0/maxa(x0)x1=A*yend;33Aitkem加速A=[-4,14,0;-5,13,0;-1,0,2];l1=0;k=1x0=[1;1;1];y0=x0/maxa(x0)x1=A*y0;y1=x1/maxa(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.01x0=x1;x1=x2;l1=l0;k=k+1x2=A*y2maxk=maxa(x2)y2=x2/maxkl0=maxa(x2)-(maxa(x2)-maxa(x1))^2/(maxa(x2)-2*maxa(x1)+maxa(x0))end;342.原点平移法思想：由矩阵论知，若为A的特征值则–a为A–aI的特征值，且特征向量相同。若1–a为A–aI的最大模特征值，且（k–a是A–aI的次最大模特征值）则对A–aI计算1–a及对应的特征向量比对A计算收敛得快，此即为原点平移法。计算1–a及特征向量的迭代公式特征向量：特征值：max(x(k))

1–a，

a+max(x(k))

1。35注：a的选取较为困难。例3设，，求最大模特征值及特征向量。解：（P315）36幂法：A=[-3,1,0;1,-3,-3;0,-3,4];x0=[0;0;1];k=1y=x0/maxa(x0)maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)))>0.01x0=x1;k=k+1y=x0/maxa(x0)maxa(x0)x1=A*yend37原点平移法：A=[-3,1,0;1,-3,-3;0,-3,4];x0=[0;0;1];k=1y=x0/maxa(x0)maxa(x0)x1=(A+4*eye(3))*ywhile(abs(maxa(x1)-maxa(x0)))>0.01x0=x1;k=k+1y=x0/maxa(x0)maxa(x0)x1=(A+4*eye(3))*yend;ymaxa(x1)-4383.对称矩阵的Rayleigh商加速法定义设A对称，x

0，则称为x关于A的Rayleigh商思想：A对称，特征值1，2，…，n均为实数，且存在特征向量v1，v2，…，vn为标准正交基。设，a1

0，则39

40当k充分大时，M'与k无关)41注:此比Aitken加速中的(3.2-6)更快公式称为Rayleigh商加速法。其中42注：有了R(x(k))，R(x(k+1))，R(x(k+2))，的值，可再用Aitken加速法得到的一个更好的近似值：因为所以43G2={z：|z–3|0.y=x0/maxa(x0)则k有Hk~H1~A。0=max(x2)-[max(x2)-max(x1)]^2/[max(x2)-2max(x1)+max(x0)]定理设ARnn，则存在正交阵Q使while(abs(maxa(x1)-maxa(x0)))>0.x1=inv(R)*z思想：A对称，特征值1，2，…，n均为实数，且存在特征向量v1，v2，…，vn为标准正交基。G3'={z：|z–0.|1-0|>epsy=x0/maxa(x0)适当选取b1，b2，…，bn就有可能使B的某些盖氏圆的半径比A的相应盖氏圆的半径小。例2用幂法求方阵A的最大模特征值，并用Aitkem加速法例4设，用Rayleigh商加速法求的最大模特征值及特征向量，并与幂法相比较。y2=x2/max(x2)设An有n个线性无关的特征向量v1，v2，…，vn，对应的特征值1，2，…，n，满足xk+1–x*=(c–k)(xk–x*)，其中1)欲缩小Gi，可取bi最大。一些工程技术问题需要用数值方法求得矩阵的全部或部分特征值及相关的特征向量。由于特征值是相似不变量，所以代数上常用相似变换将矩阵化简以得到特征向量，这里也可用相似变换将盖氏圆的半径变小，以得到更好的估计。QR方法即使用QR分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛使用的方法之一。G1={z：|z–0.例4设，用Rayleigh商加速法求的最大模特征值及特征向量，并与幂法相比较。解：（P317）幂法：A=[6,2,1;2,3,1;1,1,1];x0=[1;1;1];k=1y=x0/maxa(x0)maxa(x0)x1=A*ywhile(abs(maxa(x1)-maxa(x0)))>0.001x0=x1;k=k+1y=x0/maxa(x0)maxa(x0)x1=A*yend;44Rayleigh商加速法：A=[6,2,1;2,3,1;1,1,1];x0=[1;1;1];r=0;k=1y=x0/maxa(x0)maxa(x0)x1=A*ywhile(abs(r1-r))>0.001x0=x1;r1=r;k=k+1y=x0/maxa(x0)maxa(x0)x1=A*yr=y'*x1/(y'*y)end452.3反幂法——用A–1代替A作幂法，即反幂法1.求最小模特征值及相应的特征向量若A可逆，|1|>|2|…|n|为其特征值，则为A-1的最大模特征值。迭代公式：x(k+1)=A–1x(k)，k=0，1，2，…，但A–1不易求，通常可解方程组Ax(k+1)=x(k)来求x(k+1)即有(3.2-12)46

(3.2-12)当k

时有注：为解(3.2-12)中的方程组。对作LR分解（带行交换）PA=LR则有472.求任一特征值及相应特征向量——反幂法结合原点平移法思想：若已知为j的近似值，则的特征值是而显然非常大（最大），比值很小迭代公式：48迭代公式：当k

时有注：(1)若有LR分解，则迭代公式

(3.2-16)49(2)在(3.2-16)中直接取z(1)=(1，…，1)T作初值开始迭代称为半次迭代法50例5设的一个特征值的近似值，用带原点平移的反幂法求及相应的特征向量（见[P320]）51解：Matlab代码如下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.001x0=x1;y=x0/maxa(x0)z=inv(L)*yx1=inv(R)*zend;-6.42+1/maxa(x1)52预备知识：矩阵论1.矩阵QR分解定理设A

Rnn可逆，则存在正交阵Q与上三角阵R使A=QR注：方法1)使用史密斯正交变换方法2)使用Householder变换（反射）方法3)使Givens变换（旋转）532.矩阵Schur分解定理设A

Rnn，则存在正交阵Q使——实Schur型其中Rii至多2阶。若1阶，其元素即A的特征值，若2阶其特征值为A的一对共轭复特征值。注：想加快迭代速度通常先将A化为上Hessenberg阵543.设A

Rnn，则A正交相似于一个n阶上Hessenberg矩阵（i>j+1

hij=0）证明：见（P125）55QR方法即使用QR分解构造迭代序列，是目前求一般矩阵全部特征值的最有效并广泛使用的方法之一。3.3.1QR方法的计算公式定义：①矩阵列{Ak}，当k

时，若其对角元均收敛，且严格下三角部分元素收敛到0，则称{Ak}本质收敛到上三角阵。②矩阵列{Ak}，当k

时，若其对角子块收敛到1阶或2阶的方阵，其下部收敛到0，则称{Ak}本质收敛到块上三角阵。思想：从A1=A出发用正交相似变换得序列{Ak}，使当k

时，Ak本质收敛到块上三角阵3.3QR方法56方法：设A1=Q1R1（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）=…=(Q1…Qk)-1A(Q1…Qk)记Gk=Q1…Qk——正交，故有Ak~A，且A1Gk=GkAk+1②记Hk=Rk…R1，则Ak=GkHk——QR分解GkHk=(Q1…Qk)(Rk…R1)=Gk-1QkRkHk-1=Gk-1AkHk-1=A1Gk-1Hk-1=…=A1k=Ak57注：为求得A的特征值，只须{Ak}能趋于块上三角阵。定理3.3-1设A的特征值满足条件：|1|>|2|>…>|n|>0，vi为i对应的特征向量，i=1，2，…，n。记X=(v1，v2，…，vn)，若有直接三角分解X-1=LU（杜利特尔分解），则(3.3-1)序列{Ak}本质收敛于上三角阵，其主对角元素均为A的特征值。例1用QR方法求的A特征值，其中见（P322）注：若A不满足定理条件，{Ak}不一定本质收敛于上三角矩阵。583.2上Hessenberg矩阵的QR方法及带原点平移的QR方法在使用QR方法之前，先A将作正交相似变换化为上Hessenberg矩阵H，然后对H作QR迭代，可大量节省运算量。①Givens变换记s=sin，c=cos，则为旋转变换——正交阵。59推广到n维：称为Givens矩阵或Givens变换（旋转变换）。易知J(i，k，)为正交阵。60②对上Hessenberg矩阵用Givens变换作QR分解令hi+1,i*=sihii+c