《矩阵特征值与特征向量的计算yjs.ppt》由会员分享，可在线阅读，更多相关《矩阵特征值与特征向量的计算yjs.ppt（60页珍藏版）》请在三一办公上搜索。
1、1,第8章 矩阵特征值和特征向量的计算,很多工程计算中，会遇到特征值和特征向量的计算，如：机械、结构或电磁振动中的固有值问题；物理学中的各种临界值等。这些特征值的计算往往意义重大。,求解线性方程组的迭代法，重要一点是判断迭代法的收敛性；判断方法之一就是看迭代矩阵的特征值的模是否都小于1。,2,PA()是的高次的多项式，它的求根是很困难的。设法通过数值方法是求它的根。,通常对某个特征值，可以用些针对性的方法来求其近似值。,若要求所有的特征值，则可以对A做一系列的相似变换，“收敛”到对角阵或上（下）三角阵，,从而求得所有特征值的近似。,n阶方阵A的特征值是特征方程 PA()=det(A-E)=0
2、的根.,A的特征向量是齐次线性方程组(A-E)x=0 的非零解.,特征根和特征向量的定义（复习）,3,定理1：AR nn，1,n为A的特征值，则,（2）A的行列式值等于全体特征值之积，即,（1）A的迹数等于特征值之和，即,特征根和特征向量的基本结论。,定理2 设为AR nn的特征值且Ax=x，其中x不为0,则,（1）c为cA的特征值（c为常数且不为0);,（2）-p为A-pI 的特征值，即(A-pI)x=(-p)x;,（3）k为Ak的特征值；,（4）设A为非奇异阵，那么 且 为 特征值，即,4,定义 设矩阵A,BR nn，若有可逆阵P，使 则称A与B相似。,定理 若矩阵A,BR nn且相似，则
3、（1）A与B的特征值完全相同；（2）若x是B的特征向量，则Px便为A的特征向量。,相似矩阵及定义其性质,5,8.1 幂法和反幂法,8.1.1 幂法,幂法是用来求矩阵A按模最大的特征值和相应的特征向量的方法.也称为主特征值和主特征向量。,设A是单构矩阵,即A有n个线性无关的特征向量.,A的n个特征值为|1 2 n,对应的特征向量为1,2,n 线性无关.我们要求1 和1.,幂法的基本思想是取初始非零向量x0Rn,作迭代 xk+1=Axk=Ak+1x0,k=0,1,2,产生迭代序列xk.,由于1,2,n 线性无关,从而有 x0=11+22+nn(8.3),6,xk=Akx0=11k1+22k2+nn
4、kn,设|12n,这时,上式可写成,若10,则对充分大的k有,因而有,从而特征向量1 xk.,乘幂法的收敛速度取决于|2/1|的大小.,故有,8.1.1 幂法,7,因此，常把每一步计算的迭代向量xk规范化。,对非零向量x,用max(x)表示x的按绝对值最大的分量,称向量y=x/max(x)为向量x的规范化向量.,例如,设向量x=(2,1,-5,-1)T,则max(x)=-5,y=(-0.4,-0.2,1,0.2)T.可见规范化向量y总满足y=1.,幂法的规范化计算公式为:,任取初始向量x0=y0 0,计算,可得,实际计算时，考虑到当11时，xk的非零向量趋于无穷；当11时,xk趋于零；导致计算
5、机会出现上溢或下溢。,8.1.1 幂法,8,所以,其收敛速度由比值|2/1|来确定.,又由于,所以,因此,当k充分大时可取:1 mk,1 xk.,8.1.1 幂法,9,算法8.1 幂法,程序见p174。,(1)输入矩阵A，非零初始向量y0，最大迭代次数N，精度，置k:=0，u=0；,(2)计算 mk=max(yk)；,(3)计算,(4)若|mk-u|,则输出mk,xk,停算；,(5)若k=N，则停算，输出计算失败信息；否则，置k=k+1,u=mk,转步2；,10,用乘幂法求A的按模最大的特征值和相应特征向量.,例 8.1 设,解 取初值x0=y0=(1,1,1)T,计算得,可取 1 6.000
6、837,1(1,0.714316,-0.249895)T.,实际上,A的3个特征值分别为1=6,2=3,3=2.,11,8.1.2 加速技术,由于,所以,乘幂法收敛速度取决于比值|2/1|,当|2/1|1时,收敛是很慢的.,1.Aitken 加速方法,由上式可知,可见,序列mk线性收敛于1.,构造Aitken序列,会达到加速收敛的目的.,12,2.原点位移法,作矩阵B=A-aI,则B的特征值为qi=i-a(i=1,2,n),而且对应的特征向量相同.,如果选取a,使q1仍然是B的按模最大特征值，且满足,则对B 应用乘幂法可达到加速收敛的目的。,则对A aI计算1 a及对应的特征向量比对A计算收敛
7、得快，,此即为原点位移法。,13,2.原点位移法,程序见P176，例8.2,计算1 a及特征向量的迭代公式,，k=0,1,2,原点平移法是一个矩阵变换过程，变换简单且不破坏原矩阵的稀疏性。但由于预先不知道特征值的分布，所以应用起来有一定困难，,14,反幂法是求矩阵按模最小的特征值和相应特征向量的方法.,8.1.3 反幂法,设A是n阶非奇异矩阵,其特征值为,|1|2|n-1|n|0,对应的特征向量为1,2,n,则有A-1的特征值为,对应的特征向量为n,n-1,1.,要想求n和n只需对A-1应用乘幂法,任取初始向量x0=y00,作,15,也可将上式改写成,式(8.8)称为反幂法.显然有,每一步求x
8、k需要求解线性方程组,可采用LU分解法求解.,8.1.3 反幂法,16,(8.9),8.1.3 反幂法,17,程序见P178，例8.3,则每次迭代只需解二个三角形方程组，因此实用的公式为,8.1.3 反幂法,18,Jacobi方法是求实对称矩阵全部特征值和特征向量的一种矩阵变换方法。,8.2 Jacobi 方法,实对称矩阵A具有下列性质：,(1)A的特征值均为实数;其对应的特征向量线性无关且两两正交。,(2)存在正交矩阵Q，使QTAQ=diag(1,2,n),而且,Q的第i个列向量恰为i的特征向量;,（3）若记A1=QTAQ,则A1仍为对称矩阵.,直接找Q不大可能。我们可以构造一系列特殊形式的
9、正交阵Q1,.,Qn对A作正交变换，使得对角元素比重逐次增加，非对角元变小。,当非对角元已经小得无足轻重时，可以近似认为对角元就是A的所有特征值。Jacobi方法就是这样一类方法。,19,平面解析几何中的平面坐标旋转变换,表示平面上坐标轴旋转角的变换.,8.2.1 平面旋转矩阵(旋转正交相似变换）,在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为,一般地,在n维向量空间Rn中,沿着xi yj平面旋转角的变换矩阵为,20,称Rij()为平面旋转矩阵或Givens变换矩阵.,Rij()具有下列性质:,i j,8.2.1 平面旋转矩阵(旋转正交相似变换）,21,8.2.1 平面旋
10、转矩阵(旋转正交相似变换）,22,(2)Rij()为正交矩阵,即Rij-1()=RijT();,(3)如果A为对称矩阵,则RijT()ARij()也为对称矩阵,且与A有相同的特征值.,(4)RijT()A仅改变A的第i行与第j行元素,ARij()仅改变A的第i列与第j列元素.,8.2.1 平面旋转矩阵(旋转正交相似变换）,23,设实对称矩阵A=(apq)nn,记B=RijT()ARij()=(bpq)nn则它们元素之间有如下关系:,所以有,8.2.1 平面旋转矩阵(旋转正交相似变换）,（8.14）,24,从而,由上面两式可得,如果aij0,适当选取角,使,只需角满足,8.2.1 平面旋转矩阵(
11、旋转正交相似变换）,25,由式(8.15),令t=tan,则t满足方程,t2+2dt-1=0,为保证|/4,取绝对值较小的根,有,于是,且,8.2.1 平面旋转矩阵(旋转正交相似变换）,（8.20）,26,非对角线元素的平方和，,由,8.2.1 平面旋转矩阵(旋转正交相似变换）,得,27,我们有以下的收敛性定理保证上述计算过程。,8.2.2 Jacobi方法,28,Jacobi收敛性定理,证,则有,反复利用上式，即得,29,设k充分大时，有,这里需要说明一点:并不是对矩阵A的每一对非对角线非零元素进行一次这样的变换就能得到对角阵.因为在用变换消去 的时候,只有第 i 行,第 j 行,第 i 列
12、,第 j 列元素在变化,如果 或 为零,经变换后又往往不是零了.,因此,Qk=RT1RT2RTk 的列向量xj(j=1,2,n)为A的近似特征向量.,8.2.2 Jacobi方法,30,的全部特征值.,解 记 A0=A,取i=1,j=2,aij(0)=a12(0)=2,于是有,例 用Jacobi 方法计算对称矩阵,从而有,31,所以,再取i=2,j=3,aij(1)=a23(1)=2.020190,类似地可得,以下依次有,例题,32,从而A的特征值可取为 12.125825,28.388761,34.485401,特征向量为R1TR2TRkT,例题,33,为了减少搜索非对角线绝对值最大元素时间
13、,对经典的Jacobi方法可作进一步改进.,1.循环Jacobi方法:按(1,2),(1,3),(1,n),(2,3),(2,4),(2,n),(n-1,n)的顺序,对每个(i,j)的非零元素aij作Jacobi变换,使其零化,逐次重复扫描下去,直至S(A)为止.,2.过关Jacobi方法:取单调下降收敛于零的正数序列k,先以1为关卡值,依照1中顺序,将绝对值超过1的非对角元素零化,待所有非对角元素绝对值均不超过1时,再换下一个关卡值2,直到关卡值小于给定的精度.,8.2.2 Jacobi方法,具体算法和程序见p183，p184。,34,用Jacobi方法求得的结果精度一般都比较高，特别是求得
14、的特征向量正交性很好。所以Jacobi方法是求实对称矩阵全部特征值和特征向量的一个较好的方法。,它的弱点是计算量大，对原矩阵是稀疏矩阵，旋转变换后不能保持其稀疏的性质。,一般适用于阶数不高的矩阵.,8.2.2 Jacobi方法,35,8.3 QR方法,60年代出现的QR算法是目前计算中小型矩阵的全部特征值与特征向量的最有效方法。（实矩阵、非奇异。）理论依据：任一非奇异实矩阵都可分解成一个正交矩阵Q和一个上三角矩阵R的乘积，而且当R的对角元符号取定时，分解是唯一的。,同理可得：Ak相似于A(k=2,3,),故他们有相同特征根。,36,QR方法收敛性,37,QR方法收敛性,38,QR方法运算量很大
15、，为了减少运算量，常在使用QR方法之前把矩阵A简化为拟上三角矩阵。或称之为海森伯格矩阵(次对角元以下的元素全为零)。,8.3.2 化一般矩阵为拟上三角矩阵,形状为,39,为镜面反射矩阵，或Householder变换矩阵。,Houholder矩阵H=H(v)有如下性质：,(2),(3)记S为以v为法向量的平面,则几何上x与y=Hx关于平面S对称。因为,镜面反射变换,40,x,据前面定义和性质，有下面的定理。,证,镜面反射变换,由此可得,41,镜面反射变换,程序见P187,定理得证。,42,从而有,镜面反射变换,再利用c的值反过来计算,43,与平面旋转变换不同的是，镜面反变换可成批的消去向量的非零
16、元.,将任意矩阵A简化为海森伯格矩阵的步骤如下：,镜面反射变换,44,镜面反射变换,45,镜面反射变换,46,镜面反射变换,例：用Householder变换将矩阵A化为上Hessenberg阵,解：求Housholder矩阵,其中,47,镜面反射变换,48,镜面反射变换,用Household方法对矩阵A作正交相似变换,使A相似与上Hessenberg阵，程序见P190,49,、用 Givens变换对上Hessenberg阵作QR分解,具体步骤：,假设,（否则进行下一步）,取旋转矩阵R(1,2),50,、用 Givens变换对上Hessenberg阵作QR分解,51,、用 Givens变换对上H
17、essenberg阵作QR分解,52,、用 Givens变换对上Hessenberg阵作QR分解,53,、用 Givens变换对上Hessenberg阵作QR分解,54,、用 Givens变换对上Hessenberg阵作QR分解,因此，最多做n-1旋转变换，即得,55,、用 Givens变换对上Hessenberg阵作QR分解,因为,均为正交矩阵，故,其中Q=R21TR32TRnn-1T仍为正交矩阵。可算出完成这一过程的计算量O(4n2),比一般矩阵的QR分解的计算量O(n3)少一个数量级。,需要说明的是：通常用QR方法计算全部特征值，然后用反幂法求其相应的特征向量。,56,解：首先将A化为上Hessenberg矩阵，取,57,H即与A相似的上Hessenberg阵。,将H进行QR分解。,记B1=H,58,59,重复上述过程，迭代11次，得,算法和程序见P192，P193,60,练习题,第194页 习题88.2,8.3-8.8,
