第三章矩陣特征值與特征向量的計算

第三章矩陣特征值與特征向量的計算

第三章矩陣特征值與特征向量的計算許多實際的工程技術(shù)研究，如各種類型的振動問題，控制系統(tǒng)的穩(wěn)定性性問題等，最后往往歸結(jié)為求矩陣的特征值和特征向量。由理論可知，矩陣A的所有特征值可從特征方程此法缺乏實用意義。本文討論矩陣A的特征值和特征向量的數(shù)值方法第一節(jié)冪法和反冪法第二節(jié)瑞利（Rayleigh）商迭代法第三節(jié)雅可比（Jacobi）方法第四節(jié)QR算法第一節(jié)冪法和反冪法一、冪法冪法是用于求矩陣按模最大特征值及其相應(yīng)的特征向量的一種迭代方法，尤其適用于稀疏矩陣的情形。其算法的思想基于如下結(jié)論：定理1設(shè)矩陣，有n個線性無關(guān)的特征向量其對應(yīng)的特征值滿足對任意的非零向量，由A構(gòu)造向量序列對，有其中表示向量的第i個分量。證明：因為A有n個線性無關(guān)的特征向量，所以對任給的非零向量都可以用線性組合表示，即由A構(gòu)造向量序列由特征值定義知：故有同理可得：又由于，故對于足夠大的k，有其中，且當(dāng)時，，則由定理的證明過程可知，矩陣的特征值的計算步驟：1）先任取一非零向量Vo，一般取2）按（1）式計算3）當(dāng)k足夠大時，取由可知，在實際計算時，為避免出現(xiàn)過大或過小的數(shù)，常將迭代向量先規(guī)范化，然后再計算。其迭代過程轉(zhuǎn)化為：對于第k步證：由上式可知由時算法步驟：1）輸入矩陣A(aij)，初始向量，誤差限，最大迭代次數(shù)N.2）置3）求整數(shù)，使4）計算置5）判斷，輸出停機；否則，轉(zhuǎn)66）若，置，，轉(zhuǎn)3；否則輸出失敗信息，停機function[lambda,V]=powerl(A,X,epsilon,max1)%inputAistheNxNnonsingularmatrix%XisanNx1matrix%epsilonPisthetolerance%max1isthemaximumnumberofiterations%outputlambdaisthedominanteigenvalue%Visthedominanteigenvector%initializeparameterslambda=0;cnt=0;err=1;state=1;while((cnt<=max1)&(state==1))Y=A*X;

%normalizeY

[m,j]=max(abs(Y));c1=m;

dc=abs(lambda-c1);Y=(1/c1)*Y;%updateXandlambdaandcheckforconvergence

dv=norm(X-Y);err=max(dc,dv);X=Y;lambda=c1;state=0;

if(err>epsilon)state=1;endcnt=cnt+1;endV=X;c1=Y(j);例用規(guī)范化的冪法求矩陣的按模最大特征值和對應(yīng)特征向量初始向量為：k01234567127444.4237744.9233344.9957244.9995944.9995344.9995319514.8432214.9762314.9986514.9998814.9998314.999831-184-29.64262-29.95048-29.99722-29.99974-29.99968-29.999681111111110.346720.334130.333370.333340.333330.333330.333331-0.67153-0.66727-0.66670-0.66667-0.66667-0.66667-0.66667127444.4237744.9233344.9957244.9995944.9995344.99953討論：1）當(dāng)為矩陣的m重根時，按模求矩陣的最大特征值和向量仍成立。2）當(dāng)時則序列為擺動序列，當(dāng)k充分大時，有二、冪法的加速冪法的收斂速度是線性的，且依賴于比值原點移位法：矩陣A與的特征值有以下關(guān)系，若是A的特征值，則就是的特征值。且相應(yīng)的特征向量不變。由迭代公式：適當(dāng)?shù)剡x取，使得且：這樣用冪法計算的最大模特征值及相應(yīng)特征向量的收斂速度比對A用冪法計算要快。此法稱為原點移位法。缺點是取較困難。三、反冪法反冪法用于計算按模最小的特征值及對應(yīng)的特征向量。設(shè)A為n階非奇異實數(shù)方陣，則A-1存在。若A的特征值滿足對應(yīng)的特征向量為。因為則其滿足對應(yīng)的特征向量仍是xi(i=1,2,···,n)。因此，對A作冪法迭代可得A的模最大特征值及其相應(yīng)的特征矢量；而用A-1代替A做冪法迭代，則可得A的模最小特征值及其相應(yīng)的特征矢量。用A-1代替A的作法就稱為反冪法或逆冪法（反迭代法或逆迭代法）。則可得反冪法的迭代向量為：由于計算時比較麻煩，往往不能破壞矩陣A的一些性質(zhì)（如稀疏性），因此常用求解方程組代替迭代。由于矩陣在求解過程中保持不變，因此可使用三角分解法求解，則反冪法計算的主要步驟為：1、對A進行三角分解A=LU;2、3、解方程組用帶原點移位的反冪法來修正特征值，并求相應(yīng)的特征向量是非常有效的。（近似值已知的條件下）A-1移位反冪法算法流程1、輸入A=(aij)，近似值，初始向量，誤差限，最大容許迭代次數(shù)N。2、置3、作三角分解4、計算置5、若，則置，輸出，停機；否則轉(zhuǎn)76、若，則置，轉(zhuǎn)4；否則輸出失敗信息，停機。function[lambda,V]=invpow(A,X,alpha,epsilon,max1)%inputAistheNxNnonsingularmatrix%XisanNx1matrix%alphaisgivenshift%epsilonPisthetolerance%max1isthemaximumnumberofiterations%outputlambdaisthedominanteigenvalue%Visthedominanteigenvector%initializeparameters[n,n]=size(A);A=A-alpha*eye(n);lambda=0;cnt=0;err=1;state=1;while((cnt<=max1)&(state==1))

%SolvesystemAY=X

Y=A\X;%normalizeY[m,j]=max(abs(Y));c1=Y(j);dc=abs(lambda-c1);Y=(1/c1)*Y;%updateXandlambdaandcheckforconvergence

dv=norm(X-Y);err=max(dc,dv);X=Y;lambda=c1;state=0;

if(err>epsilon)state=1;endcnt=cnt+1;endlambda=alpha+1/lambda;V=X;例：用反冪法求矩陣，接近2.93的特征值，并求相應(yīng)的特征向量。解：對A-2.93I作三角分解得k012010.988684100.93-0.997372511.86492.07244367.959059512.69231314.278436-7.4019254-12.803851-14.2676296.883790612.83758114.266268k0123.07526883.00878783.00009547.959059512.69231314.278436第二節(jié)瑞利（Rayleigh）商迭代法定理：設(shè)A為n階實對稱方陣，其特征值滿足：且相應(yīng)的特征向量正交，即

則按規(guī)范化的向量Uk有：證明：顯然其迭代收斂速度比冪法加快了一倍。其迭代矩陣為：第三節(jié)Jacobi方法

雅可比方法是求實對稱矩陣全部特征值及對應(yīng)特征向量的一個變換方法。也是一種迭代，它的基本思想是通過一系列的正交相似變換，化實對稱矩陣為一個近似對角陣，此對角元素就是該矩陣的特征值，這些正交陣的乘積矩陣的各列就是相應(yīng)的特征向量。一、預(yù)備知識：（1）矩陣A與B均為實數(shù)方陣，若存在非奇異陣P，使得則矩陣A與B相似，A、B有相同的特征值。（2）若B是上或下三角矩陣或?qū)顷?，則B的主對角元素即是B的特征值。（3）若Q為實矩陣，且，則Q為正交矩陣。正交矩陣的乘積仍為正交矩陣。（4）實對稱矩陣的特征值為實數(shù)，且存在標(biāo)準(zhǔn)正交的特征向量系。5）對任何實對稱矩陣A，總存在正交矩陣Q，使得其中為A的特征值，的列是相應(yīng)的特征向量。6）矩陣為旋轉(zhuǎn)矩陣。iijj其幾何意義為：在n維空間中，以i,j軸形成的平面上將i,j軸旋轉(zhuǎn)一個角度對于二維坐標(biāo)系：二、經(jīng)典Jacobi方法以二階矩陣為例：旋轉(zhuǎn)矩陣：其中：為使A的相似矩陣B成為對角陣，只需適當(dāng)選取θ使P為正交矩陣可得：其中：則旋轉(zhuǎn)矩陣P確定。A的特征值為：其特征向量的計算為：設(shè)：其中為P的行向量。由：即：對應(yīng)與特征值的特征向量為：三、推廣：對于n階矩陣與二階矩陣的對角化類似，對于n階實對稱矩陣，記，對A作一系列旋轉(zhuǎn)變換，即：其中仍是對稱陣，為旋轉(zhuǎn)矩陣。它與單位矩陣的區(qū)別在于(i,i),(i,j),(j,i),(j,j)四個位置的元素不同。具有如下的性質(zhì)：（1）為正交矩陣。（2）PA只改變A的第i行與第j行元素，APT

只改變A的第i列與第j列元素，PAPT

只改變A的第i、j行,第i、j列元素。雅可比方法就是用一系列的旋轉(zhuǎn)相似變換逐漸將A化為對角陣的過程。注：旋轉(zhuǎn)矩陣中的i,j中的位置是由Ak-1中要化為零的元素決定的。為使，應(yīng)滿足：對于Ak-1中的元素其變化為1、非對角元素平方和：此外：則2、主對角元素平方和又3）變換前后全部元素的平方和不變或結(jié)論：每作一次旋轉(zhuǎn)相似變換，會使對角線元素平方和增加，矩陣的非主對角線上元素的平方和就減少一次，繼續(xù)變換直到足夠小，則主對角元既是矩陣的近似特征值，而各次旋轉(zhuǎn)矩陣的乘積之行向量即為矩陣的近似向量。由到的計算步驟為：1、由選主元法，確定i,j，使2、由計算轉(zhuǎn)角，同時確定旋轉(zhuǎn)變化陣3、計算中的元素其余元素不變。4、當(dāng)足夠小時停止運算，否則轉(zhuǎn)向（1）繼續(xù)運算。定理：設(shè)為對稱陣，是由上述經(jīng)典Jacobi法確定的矩陣，則當(dāng)時：證明：記其中為的非對角元素構(gòu)成的對稱陣。則取為Ak-1

中按模最大的非對角元素，即所以例：用Jacobi方法求矩陣：的特征值（精度達到）。解：取i=1,j=2,因為，故則：按上述方法依次變換可得：矩陣A的特征值（精確值）為：特點：1、Jacobi方法具有精度高、收斂快、算法穩(wěn)定的特點，且求得的特征向量具有較好的正交性。2、缺點是計算量大，隨n的增加而收斂減速。適用于求較低階的對稱滿秩矩陣的特征值和特征向量。1、輸入初始陣A，確定容差Epsilon2、由選主元法，確定i,j，使3、由即確定旋轉(zhuǎn)變化陣4、計算中的元素其余元素不變。5、計算非對角陣的F范數(shù)E(A)6、若E(A)<Epsilion，則主對角線值即為特征值，的各列即為相應(yīng)的特征向量；否則，重復(fù)上述過程。function[V,D]=jacobil(A,epsilon)%inputAistheNxNmatrix%epsilonPisthetolerance%Visthematrixofeigenvector%Disthematrixofeigenvalue%initializeparametersD=A;[n,n]=size(A);V=eye(n);state=1;%conculaterowpandcolumnqoftheoff-diagonalelementofgreatestmagnitude[m1,p]=max(abs(D-diag(diag(D))));[m2,q]=max(m1);p=p(q);while(state==1)t=(D(p,p)-D(q,q))/(2*D(p,q));

if(t<0)b=-(sqrt(t^2+1)-abs(t));elseb=sqrt(t^2+1)-abs(t);end

c=1/sqrt(1+b^2);s=b*c;R=[c,s;-s,c];D([p,q],:)=R'*D([p,q],:);D(:,[p,q])=D(:,[p,q])*R;V(:,[p,q])=V(:,[p,q])*R;[m1,p]=max(abs(D-diag(diag(D))));[m2,q]=max(m1);p=p(q);

delt=norm(D-diag(diag(D)),'fro')

%if(abs(D(p,q))<epsilon*sqrt(sum(diag(D).^2)/n))

if(delt<epsilon)state=0;endendD=diag(diag(D));四、過關(guān)Jacobi法經(jīng)典Jacobi法在每次相似變換前，都需從(n2-n)/2中選取按模最大的元素作為變換對象，因此比較費機時。實用中常采用過關(guān)法：1、取某一較小的數(shù)作為關(guān)，如取2、按某種順序檢查比大的元素作為相似變換消元對象，直到滿足所有非對角元中元素均小于Epson。3、再取比更小的作為關(guān)并繼續(xù)運算。4、依次類推直到滿足精度為止。第四節(jié)QR算法應(yīng)用于求任意方陣全部特征值及其特征向量的最有效方法。與Jacobi方法不同的是QR方法是先進行QR分解，然后通過逆序相乘來實現(xiàn)正交相似變換。QR分解可用平面旋轉(zhuǎn)或Householder變換來實現(xiàn)。一、Householder變換定義：設(shè)向量且滿足，則稱為Householder矩陣（簡稱H矩陣）或鏡面反射陣構(gòu)造H矩陣主要是確定向量W。對任意非零向量，可令，這時其中鏡面陣的性質(zhì)：1、是正交陣，即HTH=I2、是對稱陣，即H-1=H=HT3、H僅有兩個不等的特征值，其中1是n-1重根，-1是單重根。W為其相應(yīng)的特征向量。顯然有即向量aa’關(guān)于平面S對稱（超平面span{w}）4、H陣對向量變換時，具有鏡面反射作用：若令則有對任意非零向量，可選擇H陣使其中，。即可以將向量中除第一個分量外的其它分量全化為零。此可取則：向量變換為同理可得：為避免計算可能產(chǎn)生的誤差，取所以對向量構(gòu)造H陣的計算步驟為：1、計算2、計算向量3、計算4、寫出矩陣?yán)涸O(shè)有向量，構(gòu)造H矩陣，使，其中解：因為于是：可得二、QR分解適當(dāng)選取Householder矩陣，經(jīng)n-1步的變換可將一般的矩陣做QR分解使A=QR，其中Q是正交陣，R為上三角陣。設(shè)而為逐次選定的H矩陣，對A逐次左乘使得：變換陣H的構(gòu)造：設(shè)對A已進行了k-1次變換，其變換陣為：則第k步的變換主要使其中則取其中而這樣構(gòu)成的使得k次變換為這樣經(jīng)n-1次變換后，矩陣化為記則有對于滿秩矩陣實行QR分解時，當(dāng)n較大的，計算量較大。一般采用如下算法：1、作相似變換，將A化為上海森堡（upperHessenberg）矩陣-擬上三角矩陣

，次對角線下方全為零的特殊矩陣。注：對于實對稱陣經(jīng)變換后，H