特征值与特征向量算法的研究.doc

特征值与特征向量算法的研究摘 要：目前在很多科学领域中进行研究时，问题常会转化成特征值与特征向量的求解。本文就求解特征值、特征向量的几个重要的算法作出了研究。如：幂法，反幂法，QR算法，Jacobi迭代法等。讨论了各算法的原理及各算法在MATLAB中的运行情况，从而比较出在面对不同性质的矩阵时每个算法都各有千秋。幂法计算简单，特别适用于高阶稀疏矩阵，但其收敛速度较慢，要想加快幂法的收敛速度可采用反幂法及位移技术。QR方法被人们称为数值数学，最值得注意的算法之一，它是目前求任意矩阵全部特征值和特征向量最有效的方法。Jacobi方法是古典方法，它收敛快，精度高，便于并行计算且算法稳定。但比较适用于求低阶的对称矩阵的全部特征值和特征向量。关键词：特征值 特征向量 幂法 QR算法 雅可比算法Abstract: At present While carrying on research in a lot of scientific fields，the questions often change into how to solve the eigenvalue and eigenvector. The degree paper do research in some important arithmetic on eigenvalue and eigenvector, such as power method, inverse power method, QR arithmetic and Jacobi arithmetic etc. In this paper, we discuss the theory of arithmetic, also including how to use them in the MATLAB. Then we can come to the conclusion that the power method is easy to run, it is fit to sparse matrix, but the speed is too slow! If you want to speed the rate of convergence, you can use inverse power method. QR is diffused as numerical mathematics, one of the noteworthiness arithmetic; it is the best arithmetic which can solve all eigenvalue and eigenvector of any matrix. Jacobin arithmetic is a classicality, the rate of convergence is fast, and the precision are high too. It is easy to parallel calculate, and the result is steady but it is fit to calculate all eigenvalue and eigenvector of symmetric matrix.Keywords:Eigenvalue eigenvector power method QR arithmetic Jacobin arithmetic目录摘要(1)1绪论11问题提出与研究的目的和意义（3）12国内外研究现状（3）13论文结构与研究方法（3）14论文使用的软件环境（4）2 MATLAB语言及其在数值计算方面应用的简介 （4）21幂法（4）22反幂法（6）23移位反幂法（8）24 QR算法 （10）25雅可比（Jacobi）迭代法（12）3记单侧旋转法的对称矩阵特征值的求法（16）4几种算法的比较（16）5 MATLAB计算仿真结果 （17）在MATLAB中用幂法求其特征值与特征向量（17）6尚待深入研究的问题（17）参考文献（18）致谢（18）一、绪 论1.1问题提出与研究的目的和意义代数特征值问题是数值代数的一个重要研究领域。对它的研究具有重要的理论意义和应用价值。许多科学技术和工程问题，如结构力学中的固有频率分析问题以及控制系统的稳定性问题,最终都转化为特征值问题。代数特征值问题的解法长期在各个领域都起着不可忽视的作用，因为它充分地显示出所谓经典数学与实用数值分析之间的差异。本文就求解特征值、特征向量的几个重要的算法作出了研究。如：幂法，反幂法，QR算法，Jacobi迭代法等。同时，讨论了各算法的原理及各算法的MATLAB实现，是本文的研究目的。1.2国内外研究综述矩阵逆特征值问题(亦称矩阵反问题)涉及的领域有数学物理、地球物理、量子化学、光学、力学、结构设计、模态识别、自动控制等等.例如，弹簧一一一质量系统的参数识别，极点配置，但是，由于逆问题本身的复杂性以及理论模型与实际问题的差异，矩阵逆特征值问题的研究进展并不是很快，只是关于对称三对角线矩阵的逆特征问题的研究相对地比较成熟，有很多提法及解的适定性条件.而对于五对角线矩阵逆特征值问题的研究较少，有关结论也不多.目前，求解特征值问题的方法可分为两大类：一类称为变换方法，另一类称为向量迭代法。变换方法是直接对原矩阵进行处理，通过一系列变换，使之变成一易于求解特征值的形式，如Jacobi方法，Givens方法，QR方法以及QL方法等。变换方法由于要存储矩阵元素，因而，它只适用于求解阶数较底的矩阵，它一般和向量迭代法和起来使用。向量迭代法是通过一系列矩阵向量乘积而求得特征值和特征向量。由于向量迭代法可采用压缩存贮技术，因而它适合于求解大型矩阵尤其是大型稀疏矩阵的特征值问题。目前常用的计算大型对称稀疏矩阵特征值问题的迭代法有子空间迭代法和Lanczos方法.子空间迭代法可同时求得几个端部特征值及特征向量.子空间迭代法有一个很大的缺点,就是收敛速度较慢,从而运算量较大,并且随着迭代次数的增加,舍入误差的影响也会增大. Lanczos方法利用三项递推公式,存储量较小,计算速度较快,理论分析也较为成熟.但Lanczos方法并不是对各种类型的矩阵都非常有效,对有些问题, Lanczos方法的收敛速度是很慢的.1975年，E.R.Davidson提出了计算大型对称矩阵极端特征值的一种新方法Davidson方法。Davidson方法对于计算分子物理学中的特征值问题是非常有效的。因为在分子物理学中的矩阵一般都是对角占优矩阵，对于对角占优矩阵，Davidson方法能够使得不要求的特征向量的分量从逼近要求的特征向量的子空间的基中逐步减小，从而大大加快了特征值的收敛速度。继E.R.Davidson之后，又有许多人对Davidson方法进行了研究和改进，使之能够计算非对角占优矩阵和非对称矩阵的极端特征值。这些改进使得Davidson方法变得更为灵活和有效。不精确Newton方法也是目前计算大型对称稀疏矩阵特征值的一种有效方法。1.3论文结构与研究方法特征值和特征向量不仅是一个重要的分支，而且业已成为现代各科技领域处理大量有限维空间形式与数量关系的强有力的工具。特别是计算机的广泛应用为矩阵论的应用开辟了广阔的前景。例如：系统工程、优化方法和稳定性的理论等，都与矩阵论有密切的关系。特别在工程实践中（例如，机械的振动、气轮机叶片的振动、电磁振荡、物理学中某些临界值的确定等）特征值与特征向量的计算有着非常重要的应用，其相应的算法对不同实际问题有着略微的差异1.4论文使用的软件环境MATLAB C语言二、MATLAB语言及其在数值计算方面应用的简介2.1幂法2.1.1幂法的定义幂法是计算矩阵主特征值（指模最大的特征值）及对应的特征向量的一种向量迭代12.1.2幂法的基本思想任取初始非零向量，由矩阵A构造迭代格式=A（K=1，2，） (2.1.1)于是得迭代序列=A=A= (2.1.2).=A=由假设矩阵A有n个线性无关的特征向量（K=1，2，n）于是给定的初始向量可以用这组特征向量线性表示，即 (2.1.3)并设。把代入，得 (2.1.4) (2.1.5) 同理可得 (2.1.6)设主特征值满足。把上式改写成如下形式： 同理有 因为| (i=2,3,n) 所以 (i=2,3,n)从而可得 (2.1.7)和 (2.1.8)则说明迭代向量 为特征值所对应的特征值的近似向量.而式(8)表示向量和向量近似地线性相关,其中比例系数是一个常数,该数即为所求的主特征值.具体计算时,可以取 (2.1.9) 其中表示向量的第i个分量.2.1.3算法基本思想：设，记其中符号表示向量中模最大的分量。对于任意的非零初始向量，按迭代格式,；生成向量序列及数列2。2.1.4程序在Matlab中，为利用幂法计算矩阵的主特征值及特征向量，首先编写求数列中模最大的分量的程序2function t=maxnorm(a)n=length(a);t=0;for i=1:n if abs(a(i)/max(abs(a)=1 t=a(i); endend 调用上述程序，用幂法计算矩阵的特征值及对应的特征向量的程序如下：function mt,my=maxtr(a,eps)n=length(a);x0=diag(ones(n);k=1x=a*x0while norm(x-x0)eps k=k+1 q=x; y=x/maxnorm(x) x=a*y; x0=q;endmt=maxnorm(x)my=yRayleigh迭代法计算矩阵的特征值幂法是一种求一般矩阵主特征值的向量迭代法，若矩阵是对称的，则可以用Rayleigh商迭代法function raymaxtr(a,eps)n=length(a);x0=diag(ones(n);k=1x=a*x0;y=xry=y*xwhile norm(x-x0)eps k=k+1 q=x y=x/norm(x,2) x=a*y x0=q ry=y*xend2.2反幂法2.2.1反幂法的定义及基本思想如果矩阵A为非奇异矩阵，则存在且A的特征值均不为零。设A有n个线性无关的特征向量，它们所对应的特征值为并按模的大小排列为由可得 (2.1.10)矩阵的特征值(i=1,2,n),并有于是的主特征值为，而且所对应于的特征向量任是，因此对矩阵应用幂法求主特征值，就是对A求绝对值最小的特征值。用代替A做幂法计算，称为反幂法1。现仅考虑的主特征值为单重的情况,即有任给初始向量,可坐如下迭代: (k=1,2,) (2.1.11) 式(11)中需要计算A的逆矩阵.为了避免求逆矩阵把式(2.1.11)改写成然后把迭代向量作归一化的计算,计算公式为 (2.1.12)则 为矩阵A按模最小的特征值,特征向量为,其收敛率为.2.2.2程序（利用迭代格式：）反幂法求矩阵的按模最小的特征值和对应的特征向量function mx,mt,my=invmaxtr(a,eps)n=length(a);x0=diag(ones(n);x=inv(a)*x0;k=0;while norm(x-x0)eps k=k+1 q=x; y=x/maxnorm(x) x=inv(a)*y; x0=q;endmt=1/maxnorm(x)my=y2.3移位反幂法移位特征值需要一个特征值的较好的近似值，然后用迭代可以得到一个精确的解。首先使用其他方法，如QM法和Givens法得到一个初始近似值。这里只考虑特征值是唯一的情况3。2.3.1定理（移位特征值）设，V为A的特征对。如果是任意量，则-，V是A-I的特征对3。2.3.2定理（逆特征值）设，V为A的特征对。如果0，则V是1/，V是矩阵的特征对3。2.3.3定理 设，V为A的特征对。如果，则1/（-），V是矩阵（A-I）的特征对3。2.3.4移位反幂法的基本思想设n*n矩阵A有不同的特征值。考虑特征值可选择常量，使得=1/（-I）是的主特征值。而且，如果选择适当的，则通过下列递归公式可得到序列和： （2.3.1）和 （2.3.2）其中且 （2.3.3）这两个序列收敛到矩阵的主特征对，。最后，通过 （2.3.4）计算得到A对应的特征值3。在Matlab中用移位反幂法计算n*n矩阵A的特征值和对应的特征向量。设n个特征值满足，且是一个实数，对每个i=1，2，j-1，j+1，n满足。2.3.5程序Functionlambda,v=invpow(A,X,alpha,epsilon,max1)%Tnput A is an nxn matrix% -X is the nxl starting vector% -alpha is the given shift % -epsilon is the tolerance % -maxl is the maximum number of iterations %output-lambda is the dominant eigenvector% -v is the dominant eigenvector%Initialize the maxtrix A -alphaI and parametersn,n=size(A);A=A-alpha*eye(n);Lambda=0;Cnt=0;err=1;state=1;while (cntepsilon) state=1l end cnt=cnt+1;endlambda=alpha+1/c1;V=X;2.4 QR算法2.4.1 QR方法的定义QR方法是求任意矩阵的全部特征值的一种有效方法，它是JACOBI方法的推广3。2.4.2基本思想利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换，使其变化为一个近似的上三角矩阵来求全部特征值3。这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。82.4.3构造原理 实对称矩阵可用正交相似变换将其化为对角形矩阵，但对非对称矩阵，一般用正交相似变换化不成对角矩阵，但SCHUR分解定理给我们一个有关这方面的结果。2.4.3.1定理（实SCHUR分解定理）设矩阵ARn*n，则存在一个正交矩阵QRn*n，使QTAQ= 其中每个Bii是1*1或2*2的小矩阵，若Bii为1*1的，其元素就是A的实特征值，否则Bii的特征值是A一对共轭复特征值。 2431定理指出了求矩阵A的全部特征值也可用正交相似变换的方法来做，正交相似变换的结果虽然不是对角矩阵，而是分块三角形矩阵，但它同样能很方便地求出全部特征值，有关一般矩阵的正交相似变换，我们不加证明地给出一个结论。42.4.3.2定理设非奇异矩阵ARn*n，且有n个不同的特征值，记A=A(1)。如果对整数k，有矩阵A(k)的QR分解为A(k)=QkRk，则令A(k+1)=QTkA(k)Qk，当k时有A(k)本质上收敛于分块上三角形矩阵，这里“本质上收敛”指A(k)的主对角线上的元素或子块有确定的极限，其它元素或子块不管是否有极限。此定理给出了求解一般矩阵全部特征值的方法。由2431定理， 令A(k+1)=(Q1Q2.Qk)TA(Q1Q2.Qk)，则Qk也是正交矩阵，A(k+1)=QTkA(k)Qk说明A(k+1)也是原矩阵A的正交相似变换，从而A(k+1)与A有相同的特征值，n任意，此外，由A(k)=QkRk，则有QTkA(k)= QTkA(k)Rk=Rk，故有A(k+1)=QkRk，这说明A(k+1)可直接交换Qk与Rk的乘积顺序得到，于是可得如下QR算法。8对A(k)作QR分解A(k)=QkRk。 逆序相乘A(k)的分解矩阵，A(k)=RkQk。判别A(k+1)是否为主对角线为1*1或2*2的子块形式的分块上三角形矩阵，若是对角线上各子块的特征值为所求特征值，终止，否则k+1k，转。2.4.4QR分析从QR 算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值，那麽涉及的计算量是很大的，因此应该先对A作预处理。应用中常先对 做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素比三角矩阵多了一条次对角线，其形式为： 上Hessenberg矩阵 下Hessenberg矩阵实际上，Hessenberg矩阵虽然不是三角矩阵，但它很接近三角矩阵，由于其每列只比三角矩阵多一个非零元，故选用旋转变换做QR分解更简单些，因为对上Hessenberg矩阵H的第1列到第n-1列，依次做旋转变换使H的主对角线下的元素都变为零，则H化为上三角矩阵R了，6用矩阵描述就是记，则J为正交矩阵，解出H，可得H=J-1R，因为J-1也是正交矩阵，于是得H的QR分解容易验证按这个方法做对H做QR分解，然后使用QR算法则构造的迭代序列都是上Hessenberg矩阵中进行，于是整个QR算法都在上Hessenberg矩阵中进行，这当然使QR算法的计算量大量减少。2.4.5程序应用QR算法求出矩阵全部特征值后，可以再应用反幂法进一步求特征向量function qrmaxtr(a,e)k=1n=length(a);a0=zeros(n)while norm(diag(a-a0)e k=k+1 a0=a; q,r=qr(a) a=r*q;endt=diag(a)2.5雅可比（Jacobi）迭代法2.5.1可比迭代格式雅可比迭代计算元线性方程组 （2.5.1）写成矩阵形式为。若将式（ 2.5.1）中每个方程的留在方程左边，其余各项移到方程右边；方程两边除以则得到下列同解方程组：记，构造迭代形式或 （2.5.2）迭代计算式（2.5.2）称为简单迭代或雅可比迭代5。任取初始向量 ，由式（2.5.2）可得到迭代向量序列2.5.2雅可比迭代矩阵设由，得到等价方程：记不难看出，正是迭代式（4.3）的迭代矩阵，是常数项向量。于是式（2.5.2）可写成矩阵形式：9 （2.5.3）其中：2.5.3雅可比迭代收敛条件对于方程组，构造雅可比迭代格式其中，。当迭代矩阵的谱半径时，迭代收敛，这是收敛的充分必要条件。迭代矩阵的某范数时，迭代收敛。要注意的是范数小于1只是判断迭代矩阵收敛的充分条件，当迭代矩阵的一种范数|B|1，并不能确定迭代矩阵是收敛还是发散。6例如，则，但它的特征值是0.9和0.8。是收敛矩阵。当方程组的系数矩阵具有某些性质时，可直接判定由它生成的雅可比迭代矩阵是收敛的。定理4.3 若方程组的系数矩阵，满足下列条件之一，则其雅可比迭代法是收敛的。（1）为行对角占优阵，即（2）为列对角占优阵，即2.5.4程序functionV,D=jacobil(A,epsilon)%Input -A is an nxn matrix% -epsilon is the tolerance%Output -V is the num matrix of eigenvectors% -D is the diagonal nxn matrix of eigenvalues%Initialize V,D,and parametersD=A;n,n=size(A);V=eye(n);state=1;%Calculate row p and column q of the off-diagonal element of greatest magnitude in Am1 p=max(abs(D-diag(diag(D);m2 q=max(m1);p=p(q);while(state=1) t=D(p,q)/(D(q,q)-D(p,p); c=1/sqrt(t2+1); s=c*t; 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); if(abs(D(p,q)epsilon*sqrt(sum(diag(D).2)/n) state=0; endendD=diag(diag(D);三、记单侧旋转法的对称矩阵特征值的求法3.1分析求解对称方阵特征值及特征向量的常用方法是雅可比算法。但雅可比算法在每次消去计算前都要对非主对角元素选择最大元，导致非实际计算开销很大。在消去计算时，必须对方阵同时进行行、列旋转变换，这里称之为双侧旋转变换。在双侧旋转变换中，方阵的行、列方向都有数据相关关系，使得整个计算中的相关关系复杂，在求解高阶矩阵时难于并行化7。针对雅可比算法的缺点，将双侧旋转改为单侧旋转，提出了一种求对称矩阵特征值的快速算法。这一算法对矩阵仅实施列变换，使得数据相关关系仅在同列之间，因此方便数据划分，适合并行计算。3.2算法的基本思想设A为一对称方阵，是A的特征值，x是A的特征向量，即有：Axx。因A=AT，则ATAxAxx，说明了2是ATA的特征值，x是ATA的特征向量，即ATA的特征值是A的特征值的平方，并且它们的特征向量相同。如果对A进行列变换，得到方阵Q使其各列两两正交，即AV=Q，这里V为表示正交化过程的正交方阵。因QTQ(AV)TAV=VTATAV=ding(1,2,，n)，为n阶对角方阵。可见这里1,2,，n是矩阵ATA的特征值，V的各列是ATA的特征向量。由此可得：i i2（i=1,2,n）。设Q的第i列为q，则qiTqi=i，|qi|2（qiTqi）1/2因此在将A进行列变换得到各列两两正交的方阵Q后，其各列的谱范数|qi|2即为A的特征值的绝对值。设V的第i列为vi，由于Av=Q，则有Avi=qi。又因vi为A的特征向量，即有Aviivi，则iviqi。因此i的符号可由向量qi及vi的对应分量是否同号判别，实际上在具体算法中只要判断它们的第一个分量是否同号即可。若相同，则i取正值，否则取负值。3.3与经典雅可比的比较用雅可比算法求对称矩阵特征值，会遇到额外计算开销大，数据相关关系复杂、不易并行等缺陷，针对以上缺陷该文提出了基于单侧旋转方法求对称矩阵特征值的算法。这一算法不但收敛速度快而且对矩阵仅实施列变换，数据相关关系仅在同列之间，因此方便数据划分，容易并行。四、几种算法的比较本论文介绍的幂法，反幂法，QR方法及Jacobi方法是求矩阵特征值和特征向量的常用数值方法，它们都是造构造迭代产生的矩阵序列来达到目的的。幂法计算简单，特别适用于高阶稀疏矩阵，但其收敛速度不能令人满意，要想加快幂法的收敛速度可采用反幂法及位移技术。 QR方法是60年代发展起来的，被人们称为数值数学，最值得注意的算法之一，它是目前求任意矩阵全部特征值和特征向量最有效的方法。Jacobi方法是古典方法，它收敛快，精度高，便于并行计算且算法稳定。用Jacobi方法求出的特征向量在较好的正交性，不过它的计算量较大，当阶数n增大时收敛速度减慢，因此Jacobi方法适用于求低阶的对称矩阵的全部特征值和特征向量。五、MATLAB计算仿真结果在MATLAB中用幂法求其特征值与特征向量设矩阵 在MATLAB中输入命令：a=2 -1 0;-1 2 -1;0 -1 2