矩陣特征值與特征向量的計算課件

數(shù)值計算方法矩陣特征值與特征向量的計算一些工程技術(shù)問題需要用數(shù)值方法求得矩陣的全部或部分特征值及相關(guān)的特征向量。21特征值的估計較粗估計(A)||A||欲將復(fù)平面上的特征值一個個用圓盤圍起來。1.1蓋氏圓定義1-1設(shè)A=[aij]nn，稱由不等式所確定的復(fù)區(qū)域為A的第i個蓋氏圓，記為Gi：i=1，2，…，n。定理3.1-1若為A的特征值，則3定理1-1若為A的特征值，則證明：設(shè)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個蓋氏圓上，但未說明每個圓盤內(nèi)都有一個特征值。G1G2G3G451.2蓋氏圓的連通部分稱相交蓋氏圓之并構(gòu)成的連通部分為連通部分。孤立的蓋氏圓本身也為一個連通部分。定理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()的特征多項式的系數(shù)是的多項式，從而A()的特征值1()，2()，…，n()為的連續(xù)函數(shù)。7計算max(x0)，y0=x0/max(x0)求最小模特征值及相應(yīng)的特征向量冪法是求方陣的最大特征值及對應(yīng)特征向量的一種迭代法。G4={z：|z+4|0.1)欲縮小Gi，可取bi最大。注：有了R(x(k))，R(x(k+1))，R(x(k+2))，的值，可再用Aitken加速法得到的一個更好的近似值：(1)思想：由定理3.須采用“規(guī)范化”的方法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)輸入數(shù)組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為終點的連續(xù)曲線，且始終不會越過Gi；aiii8不失一般性，設(shè)A開頭的k個圓盤是連通的，其并集為S，它與后n–k個圓盤嚴格分離，顯然，A()的前k個蓋氏圓盤與后n–k個圓盤嚴格分離。當=0時，A(0)=D的前k個特征值剛好落在前k個圓盤G1，…，Gk中，而另n–k個特征值則在區(qū)域S之外，從0變到1時，與始終分離（嚴格）。連續(xù)曲線始終在S中，所以S中有且僅有A的k個特征值。9注：1)每個孤立圓中恰有一個特征值。2)例1中G2，G4為僅由一個蓋氏圓構(gòu)成的連通部分，故它們各有一個特征值，而G1，G3構(gòu)成的連通部分應(yīng)含有兩個特征值。3)因為例1中A為實方陣，所以若為A的特征值，則也是A的特征值，所以G2，G4中各有一個實特征值。101.3蓋氏圓與相似變換由于特征值是相似不變量，所以代數(shù)上常用相似變換將矩陣化簡以得到特征向量，這里也可用相似變換將蓋氏圓的半徑變小，以得到更好的估計。原理：取對角陣作相似變換陣：P=diag(b1，b2，…，bn)其中bi>0，i=1，2，…，n則與A有相同特征值.而B的第i個蓋氏圓為：，

11而B的第i個蓋氏圓為：，

適當選取b1，b2，…，bn就有可能使B的某些蓋氏圓的半徑比A的相應(yīng)蓋氏圓的半徑小。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最?。喝1=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冪法是求方陣的最大特征值及對應(yīng)特征向量的一種迭代法。2.1冪法設(shè)An有n個線性無關(guān)的特征向量v1，v2，…，vn，對應(yīng)的特征值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及對應(yīng)特征向量的近似值。2.規(guī)范化在實際計算中，若|1|>1則|1ka1|，若|1|<1則|1ka1|0都將停機。須采用“規(guī)范化”的方法，k=0，1，2，…(3.2-4)17規(guī)范化：，k=0，1，2，…(3.2-4)定理3.2-1任給初始向量x(0)0，有(3.2-5)18易知A()的特征多項式的系數(shù)是的多項式，A=[-3,1,0;1,-3,-3;0,-3,4];解：Matlab代碼如下注：想加快迭代速度通常先將A化為上Hessenberg陣注：1)每個孤立圓中恰有一個特征值。G4={z：|z+4|0.欲將復(fù)平面上的特征值一個個用圓盤圍起來。求任一特征值及相應(yīng)特征向量while(abs(maxa(x1)-maxa(x0)))>0.定理設(shè)ARnn，則存在正交陣Q使迭代公式：x(k+1)=A–1x(k)，k=0，1，2，…，因為{v1，v2，…，vn}為Cn的一組基，所以而B的第i個蓋氏圓為：，矩陣特征值與特征向量的計算解：Matlab代碼如下輸入數(shù)組x0,eps,A②矩陣列{Ak}，當k時，若其對角子塊收斂到1階或2階的方陣，其下部收斂到0，則稱{Ak}本質(zhì)收斂到塊上三角陣。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)，冪法收斂性的分析較復(fù)雜，但若1=2=…=r，且|1|>|r+1|…|n|，則定理結(jié)論仍成立。此時不同初始向量的迭代向量序列一般趨向于1的不同特征向量。203.算法求maxa(x)的流程，設(shè)數(shù)組x(n)存放向量x的n個分量冪法流程：數(shù)組x=[n]k=1for(i=2ton,i++)若|x[i]|>|x[k]|Tk=imax=x[k]輸入數(shù)組x0,eps,Ax1=x0y=x1/maxa(x1)x0=Ay|maxa(x1)–maxa(x0)|>eps輸出y,maxa(x0)21Matlab冪法流程：輸入數(shù)組x0,eps,Ay=x0/maxa(x0)x1=Ay|maxa(x1)–maxa(x0)|>epsx0=x1y=x0/maxa(x0)x1=Ay輸出y,maxa(x1)22例1，用冪法求的最大模特征值及對應(yīng)特征向量（見P312）解：首先給出函數(shù)代碼：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)輸入數(shù)組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īng)蓋氏圓的半徑小。y=x0/maxa(x0)一些工程技術(shù)問題需要用數(shù)值方法求得矩陣的全部或部分特征值及相關(guān)的特征向量。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及對應(yīng)特征向量的近似值。x0=x1;k=k+1|1-0|>eps2-12)中的方程組。若|1|>1則|1ka1|，若|1|<1則|1ka1|03蓋氏圓與相似變換注:此比Aitken加速中的(3.2上Hessenberg矩陣的QR方法及帶原點平移的QR方法由于特征值是相似不變量，所以代數(shù)上常用相似變換將矩陣化簡以得到特征向量，這里也可用相似變換將蓋氏圓的半徑變小，以得到更好的估計。解：A的三個蓋氏圓分別為：y=x0/maxa(x0)x0=[1;1;1];k=1設(shè){xk}線性收斂到x*，即存在c，|c|<1，滿足注：冪法的收斂速度取決于比值|2|/|1|，考慮收斂加速得上Hessenberg陣列{Hk}。輸入數(shù)組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加速法設(shè){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及對應(yīng)的特征向量比對A計算收斂得快，此即為原點平移法。計算1–a及特征向量的迭代公式特征向量：特征值：max(x(k))

1–a，

a+max(x(k))

1。35注：a的選取較為困難。例3設(shè)，，求最大模特征值及特征向量。解：（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商加速法定義設(shè)A對稱，x

0，則稱為x關(guān)于A的Rayleigh商思想：A對稱，特征值1，2，…，n均為實數(shù)，且存在特征向量v1，v2，…，vn為標準正交基。設(shè)，a1

0，則39

40當k充分大時，M'與k無關(guān))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)]定理設(shè)ARnn，則存在正交陣Q使while(abs(maxa(x1)-maxa(x0)))>0.x1=inv(R)*z思想：A對稱，特征值1，2，…，n均為實數(shù)，且存在特征向量v1，v2，…，vn為標準正交基。G3'={z：|z–0.|1-0|>epsy=x0/maxa(x0)適當選取b1，b2，…，bn就有可能使B的某些蓋氏圓的半徑比A的相應(yīng)蓋氏圓的半徑小。例2用冪法求方陣A的最大模特征值，并用Aitkem加速法例4設(shè)，用Rayleigh商加速法求的最大模特征值及特征向量，并與冪法相比較。y2=x2/max(x2)設(shè)An有n個線性無關(guān)的特征向量v1，v2，…，vn，對應(yīng)的特征值1，2，…，n，滿足xk+1–x*=(c–k)(xk–x*)，其中1)欲縮小Gi，可取bi最大。一些工程技術(shù)問題需要用數(shù)值方法求得矩陣的全部或部分特征值及相關(guān)的特征向量。由于特征值是相似不變量，所以代數(shù)上常用相似變換將矩陣化簡以得到特征向量，這里也可用相似變換將蓋氏圓的半徑變小，以得到更好的估計。QR方法即使用QR分解構(gòu)造迭代序列，是目前求一般矩陣全部特征值的最有效并廣泛使用的方法之一。G1={z：|z–0.例4設(shè)，用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.求最小模特征值及相應(yīng)的特征向量若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.求任一特征值及相應(yīng)特征向量——反冪法結(jié)合原點平移法思想：若已知為j的近似值，則的特征值是而顯然非常大（最大），比值很小迭代公式：48迭代公式：當k

時有注：(1)若有LR分解，則迭代公式

(3.2-16)49(2)在(3.2-16)中直接取z(1)=(1，…，1)T作初值開始迭代稱為半次迭代法50例5設(shè)的一個特征值的近似值，用帶原點平移的反冪法求及相應(yīng)的特征向量（見[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預(yù)備知識：矩陣論1.矩陣QR分解定理設(shè)A

Rnn可逆，則存在正交陣Q與上三角陣R使A=QR注：方法1)使用史密斯正交變換方法2)使用Householder變換（反射）方法3)使Givens變換（旋轉(zhuǎn)）532.矩陣Schur分解定理設(shè)A

Rnn，則存在正交陣Q使——實Schur型其中Rii至多2階。若1階，其元素即A的特征值，若2階其特征值為A的一對共軛復(fù)特征值。注：想加快迭代速度通常先將A化為上Hessenberg陣543.設(shè)A

Rnn，則A正交相似于一個n階上Hessenberg矩陣（i>j+1

hij=0）證明：見（P125）55QR方法即使用QR分解構(gòu)造迭代序列，是目前求一般矩陣全部特征值的最有效并廣泛使用的方法之一。3.3.1QR方法的計算公式定義：①矩陣列{Ak}，當k

時，若其對角元均收斂，且嚴格下三角部分元素收斂到0，則稱{Ak}本質(zhì)收斂到上三角陣。②矩陣列{Ak}，當k

時，若其對角子塊收斂到1階或2階的方陣，其下部收斂到0，則稱{Ak}本質(zhì)收斂到塊上三角陣。思想：從A1=A出發(fā)用正交相似變換得序列{Ak}，使當k

時，Ak本質(zhì)收斂到塊上三角陣3.3QR方法56方法：設(shè)A1=Q1R1（QR分解），令A(yù)2=R1Q1，設(shè)A2=Q2R2，令A(yù)3=R2Q2，即

k=1，2，…(3.3-1){Ak}的性質(zhì)：①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設(shè)A的特征值滿足條件：|1|>|2|>…>|n|>0，vi為i對應(yīng)的特征向量，i=1，2，…，n。記X=(v1，v2，…，vn)，若有直接三角分解X-1=LU（杜利特爾分解），則(3.3-1)序列{Ak}本質(zhì)收斂于上三角陣，其主對角元素均為A的特征值。例1用QR方法求的A特征值，其中見（P322）注：若A不滿足定理條件，{Ak}不一定本質(zhì)收斂于上三角矩陣。583.2上Hessenberg矩陣的QR方法及帶原點平移的QR方法在使用QR方法之前，先A將作正交相似變換化為上Hessenberg矩陣H，然后對H作QR迭代，可大量節(jié)省運算量。①Givens變換記s=sin，c=cos，則為旋轉(zhuǎn)變換——正交陣。59推廣到n維：稱為Givens矩陣或Givens變換（旋轉(zhuǎn)變換）。易知J(i，k，)為正交陣。60②對上Hessenberg矩陣用Givens變換作QR分解令hi+1,i*=sihii+c