第三章矩阵特征值与特征向量的计算

第三章矩阵特征值与特征向量的计算

第三章矩阵特征值与特征向量的计算许多实际的工程技术研究，如各种类型的振动问题，控制系统的稳定性性问题等，最后往往归结为求矩阵的特征值和特征向量。由理论可知，矩阵A的所有特征值可从特征方程此法缺乏实用意义。本文讨论矩阵A的特征值和特征向量的数值方法第一节幂法和反幂法第二节瑞利（Rayleigh）商迭代法第三节雅可比（Jacobi）方法第四节QR算法第一节幂法和反幂法一、幂法幂法是用于求矩阵按模最大特征值及其相应的特征向量的一种迭代方法，尤其适用于稀疏矩阵的情形。其算法的思想基于如下结论：定理1设矩阵，有n个线性无关的特征向量其对应的特征值满足对任意的非零向量，由A构造向量序列对，有其中表示向量的第i个分量。证明：因为A有n个线性无关的特征向量，所以对任给的非零向量都可以用线性组合表示，即由A构造向量序列由特征值定义知：故有同理可得：又由于，故对于足够大的k，有其中，且当时，，则由定理的证明过程可知，矩阵的特征值的计算步骤：1）先任取一非零向量Vo，一般取2）按（1）式计算3）当k足够大时，取由可知，在实际计算时，为避免出现过大或过小的数，常将迭代向量先规范化，然后再计算。其迭代过程转化为：对于第k步证：由上式可知由时算法步骤：1）输入矩阵A(aij)，初始向量，误差限，最大迭代次数N.2）置3）求整数，使4）计算置5）判断，输出停机；否则，转66）若，置，，转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);例用规范化的幂法求矩阵的按模最大特征值和对应特征向量初始向量为：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）当为矩阵的m重根时，按模求矩阵的最大特征值和向量仍成立。2）当时则序列为摆动序列，当k充分大时，有二、幂法的加速幂法的收敛速度是线性的，且依赖于比值原点移位法：矩阵A与的特征值有以下关系，若是A的特征值，则就是的特征值。且相应的特征向量不变。由迭代公式：适当地选取，使得且：这样用幂法计算的最大模特征值及相应特征向量的收敛速度比对A用幂法计算要快。此法称为原点移位法。缺点是取较困难。三、反幂法反幂法用于计算按模最小的特征值及对应的特征向量。设A为n阶非奇异实数方阵，则A-1存在。若A的特征值满足对应的特征向量为。因为则其满足对应的特征向量仍是xi(i=1,2,···,n)。因此，对A作幂法迭代可得A的模最大特征值及其相应的特征矢量；而用A-1代替A做幂法迭代，则可得A的模最小特征值及其相应的特征矢量。用A-1代替A的作法就称为反幂法或逆幂法（反迭代法或逆迭代法）。则可得反幂法的迭代向量为：由于计算时比较麻烦，往往不能破坏矩阵A的一些性质（如稀疏性），因此常用求解方程组代替迭代。由于矩阵在求解过程中保持不变，因此可使用三角分解法求解，则反幂法计算的主要步骤为：1、对A进行三角分解A=LU;2、3、解方程组用带原点移位的反幂法来修正特征值，并求相应的特征向量是非常有效的。（近似值已知的条件下）A-1移位反幂法算法流程1、输入A=(aij)，近似值，初始向量，误差限，最大容许迭代次数N。2、置3、作三角分解4、计算置5、若，则置，输出，停机；否则转76、若，则置，转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的特征值，并求相应的特征向量。解：对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第二节瑞利（Rayleigh）商迭代法定理：设A为n阶实对称方阵，其特征值满足：且相应的特征向量正交，即

则按规范化的向量Uk有：证明：显然其迭代收敛速度比幂法加快了一倍。其迭代矩阵为：第三节Jacobi方法

雅可比方法是求实对称矩阵全部特征值及对应特征向量的一个变换方法。也是一种迭代，它的基本思想是通过一系列的正交相似变换，化实对称矩阵为一个近似对角阵，此对角元素就是该矩阵的特征值，这些正交阵的乘积矩阵的各列就是相应的特征向量。一、预备知识：（1）矩阵A与B均为实数方阵，若存在非奇异阵P，使得则矩阵A与B相似，A、B有相同的特征值。（2）若B是上或下三角矩阵或对角阵，则B的主对角元素即是B的特征值。（3）若Q为实矩阵，且，则Q为正交矩阵。正交矩阵的乘积仍为正交矩阵。（4）实对称矩阵的特征值为实数，且存在标准正交的特征向量系。5）对任何实对称矩阵A，总存在正交矩阵Q，使得其中为A的特征值，的列是相应的特征向量。6）矩阵为旋转矩阵。iijj其几何意义为：在n维空间中，以i,j轴形成的平面上将i,j轴旋转一个角度对于二维坐标系：二、经典Jacobi方法以二阶矩阵为例：旋转矩阵：其中：为使A的相似矩阵B成为对角阵，只需适当选取θ使P为正交矩阵可得：其中：则旋转矩阵P确定。A的特征值为：其特征向量的计算为：设：其中为P的行向量。由：即：对应与特征值的特征向量为：三、推广：对于n阶矩阵与二阶矩阵的对角化类似，对于n阶实对称矩阵，记，对A作一系列旋转变换，即：其中仍是对称阵，为旋转矩阵。它与单位矩阵的区别在于(i,i),(i,j),(j,i),(j,j)四个位置的元素不同。具有如下的性质：（1）为正交矩阵。（2）PA只改变A的第i行与第j行元素，APT

只改变A的第i列与第j列元素，PAPT

只改变A的第i、j行,第i、j列元素。雅可比方法就是用一系列的旋转相似变换逐渐将A化为对角阵的过程。注：旋转矩阵中的i,j中的位置是由Ak-1中要化为零的元素决定的。为使，应满足：对于Ak-1中的元素其变化为1、非对角元素平方和：此外：则2、主对角元素平方和又3）变换前后全部元素的平方和不变或结论：每作一次旋转相似变换，会使对角线元素平方和增加，矩阵的非主对角线上元素的平方和就减少一次，继续变换直到足够小，则主对角元既是矩阵的近似特征值，而各次旋转矩阵的乘积之行向量即为矩阵的近似向量。由到的计算步骤为：1、由选主元法，确定i,j，使2、由计算转角，同时确定旋转变化阵3、计算中的元素其余元素不变。4、当足够小时停止运算，否则转向（1）继续运算。定理：设为对称阵，是由上述经典Jacobi法确定的矩阵，则当时：证明：记其中为的非对角元素构成的对称阵。则取为Ak-1

中按模最大的非对角元素，即所以例：用Jacobi方法求矩阵：的特征值（精度达到）。解：取i=1,j=2,因为，故则：按上述方法依次变换可得：矩阵A的特征值（精确值）为：特点：1、Jacobi方法具有精度高、收敛快、算法稳定的特点，且求得的特征向量具有较好的正交性。2、缺点是计算量大，随n的增加而收敛减速。适用于求较低阶的对称满秩矩阵的特征值和特征向量。1、输入初始阵A，确定容差Epsilon2、由选主元法，确定i,j，使3、由即确定旋转变化阵4、计算中的元素其余元素不变。5、计算非对角阵的F范数E(A)6、若E(A)<Epsilion，则主对角线值即为特征值，的各列即为相应的特征向量；否则，重复上述过程。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));四、过关Jacobi法经典Jacobi法在每次相似变换前，都需从(n2-n)/2中选取按模最大的元素作为变换对象，因此比较费机时。实用中常采用过关法：1、取某一较小的数作为关，如取2、按某种顺序检查比大的元素作为相似变换消元对象，直到满足所有非对角元中元素均小于Epson。3、再取比更小的作为关并继续运算。4、依次类推直到满足精度为止。第四节QR算法应用于求任意方阵全部特征值及其特征向量的最有效方法。与Jacobi方法不同的是QR方法是先进行QR分解，然后通过逆序相乘来实现正交相似变换。QR分解可用平面旋转或Householder变换来实现。一、Householder变换定义：设向量且满足，则称为Householder矩阵（简称H矩阵）或镜面反射阵构造H矩阵主要是确定向量W。对任意非零向量，可令，这时其中镜面阵的性质：1、是正交阵，即HTH=I2、是对称阵，即H-1=H=HT3、H仅有两个不等的特征值，其中1是n-1重根，-1是单重根。W为其相应的特征向量。显然有即向量aa’关于平面S对称（超平面span{w}）4、H阵对向量变换时，具有镜面反射作用：若令则有对任意非零向量，可选择H阵使其中，。即可以将向量中除第一个分量外的其它分量全化为零。此可取则：向量变换为同理可得：为避免计算可能产生的误差，取所以对向量构造H阵的计算步骤为：1、计算2、计算向量3、计算4、写出矩阵例：设有向量，构造H矩阵，使，其中解：因为于是：可得二、QR分解适当选取Householder矩阵，经n-1步的变换可将一般的矩阵做QR分解使A=QR，其中Q是正交阵，R为上三角阵。设而为逐次选定的H矩阵，对A逐次左乘使得：变换阵H的构造：设对A已进行了k-1次变换，其变换阵为：则第k步的变换主要使其中则取其中而这样构成的使得k次变换为这样经n-1次变换后，矩阵化为记则有对于满秩矩阵实行QR分解时，当n较大的，计算量较大。一般采用如下算法：1、作相似变换，将A化为上海森堡（upperHessenberg）矩阵-拟上三角矩阵

，次对角线下方全为零的特殊矩阵。注：对于实对称阵经变换后，H