Ch8 矩阵特征值计算.ppt

第8章矩阵特征值计算,8.1特征值性质和估计8.2幂法及反幂法8.3正交变换与矩阵分解8.4QR方法,1,8.1特征值性质和估计,工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题.,2,8.1.1特征值问题及性质,求A的特征值问题等价于求A的特征方程,设矩阵ARnn，特征值问题是求C和非零向量xRn，使,的根.,Ax=x,其中是矩阵A的特征值，x是矩阵A属于特征值的特征向量.,3,定理1设为ARnn的特征值,且Ax=x,x0，则有,(2)-为A-I的特征值,且(A-I)x=(-)x;,(1)c为的cA特征值(c0为常数),且(cA)x=(c)x;,(3)k为Ak的特征值,且Akx=kx；,(4)若A为非奇异矩阵,则-1为A-1的特征值,且A-1x=-1x.,特征值的基本性质.,4,定理2(1)设矩阵ARnn可对角化,即存在非奇异矩阵P使,的充分必要条件是A具有n个线性无关的特征向量.,(2)如果矩阵ARnn有m个(mn)不同的特征值1,2,m，则对应的特征向量x1,x2,xm线性无关.,5,定理3设ARnn为对称矩阵,则,(3)存在一个正交矩阵P使的,且1,2,n为A的特征值,而P(u1,u2,un)的列向量uj为A的对应于j的单位特征向量.,(1)A的特征值均为实数;,(2)A有n个线性无关的特征向量;,6,定义设ARnn为对称矩阵,对于任一非零向量x0,称,为对应于向量x的瑞利(Rayleigh)商.,定理4设ARnn为对称矩阵(其特征值依次记为12n),则,(1).(对任何非零xRn);,(2).;,7,证明只证(1).,由于A为实对称矩阵,可将1,2,n对应的特征向量x1,x2,xn正交规范化,即有(xi,xj)=ij,设x0为Rn中任一向量,则,于是,瑞利商必位于n和1之间,(1)成立.,8,8.1.2特征值估计与扰动,定义1设n阶矩阵A=(aij)，令,(1);,(2)集合称为复平面上以aii为圆心，以ri为半径的n阶矩阵A的n个格什戈林(Gerschgorin)圆盘.,9,定理5(Gerschgorin圆盘定理),特别地，如果A的一个圆盘Di是与其它圆盘分离(即孤立圆盘)，则Di中精确地包含A的一个特征值.,(1)设n阶矩阵A(aij)，则A的每一个特征值必属于下面某个圆盘之中,(2)如果A有m个圆盘组成一个连通的并集S，且S与余下n-m个圆盘是分离的，则S内恰包含A的m个特征值.,或者说A的特征值都在n个圆盘的并集中.,10,证明只就(1)给出证明.,Ax=x，其中x=(x1,x2,xn)T0.,或,记,考虑Ax=x的第k个方程,即,于是,即,设为A的特征值,11,A的每一个特征值必位于A的一个圆盘中,并且相应的特征值一定位于第k个圆盘中(其中k是对应特征向量x绝对值最大的分量的下标).,利用相似矩阵性质,有时可以获得A的特征值进一步的估计,即适当选取非奇异对角阵,并做相似变换.适当选取可使某些圆盘半径及连通性发生变化.,12,例1估计矩阵的特征值范围.,解矩阵A的3个圆盘为,由圆盘定理可知,A的3个特征值位于3个圆盘的并集中,由于D1是孤立圆盘,所以D1内恰好包含A的一个特征值1(为实特征值),即,A的其它两个特征值2,3包含在D2,D3的并集中.,13,取对角阵,做相似变换,矩阵A1的3个圆盘为,3个圆盘都是孤立圆盘,所以,每一个圆盘都包含A的一个特征值(为实特征值),且有估计,14,定理6(Bauer-Fike定理)设是A+IRnn的一个特征值,且P-1AP=D=diag(1,2,n),则有,其中|p为矩阵的p范数，p=1,2,.,15,证明(A)时显然成立,故只考虑(A).这时D-I非奇异，设x是A+I对应于的特征向量，由(A+I-I)x=0左乘P-1可得,而对角矩阵(D-I)-1的范数为,所以有,这时总有(A)中的一个取到m值.,16,cond(P)是特征值扰动的放大系数,但将A对角化的相似变化矩阵不是唯一的,所以取cond(P)的下确界,称为特征值问题的条件数.只要v(A)不很大,矩阵微小扰动只带来特征值的微小扰动.但是v(A)难以计算,有时只对一个P,用cond(P)代替v(A).,特征值问题的条件数和解线性方程组时的矩阵是两个不同的概念.例:二阶矩阵A=diag(1,10-10),有v(A)=1,但解线性方程组的矩阵条件数cond(A)=1010.,17,计算矩阵A的特征值.当n=2,3时,可按行列式展开的办法求特征方程p()=0的根.但当n较大时,如果按展开行列式的办法,首先求出p()的系数,再求p()的根,工作量就非常大,用这种办法求矩阵特征值是不切实际的,由此需要研究A的特征值及特征向量的数值方法.,介绍(计算机上)两类方法:1.幂法及反幂法(迭代法)，2.正交相似变换的方法(变化法).,18,8.2幂法及反幂法,幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法.幂法是计算矩阵的按模最大的特征值和相应特征向量的一种向量迭代法,特别适用于大型稀疏矩阵.反幂法是计算非奇异(可逆)矩阵按模最小的特征值和相应特征向量的一种向量迭代法,特别是计算海森伯格阵或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一.,19,8.2.1幂法(又称乘幂法),现讨论求1及x1的方法.,设实矩阵A=(aij)有一个完全的特征向量组,即A有n个线性无关的特征向量,设矩阵A的特征值为1,2,n,相应的特征向量为x1,x2,xn.已知A的主特征值1是实根,且满足条件,20,幂法的基本思想是:任取非零的初始向量v0,由矩阵A构造一向量序列vk,称为迭代向量，由假设,v0可唯一表示为,21,于是,其中,由假设故从而,为1的特征向量.,22,所以当k充分大时,有,则vk为矩阵A的对应特征值1的一个近似特征向量.,用(vk)i表示vk的第i个分量,则当k充分大时,有,即为A的主特征值1的近似值.,由于,这种由已知非零向量v0及矩阵A的乘幂Ak构造向量序列vk以计算A的主特征值1及相应特征向量的方法就称为(乘)幂法.,23,幂法的思想:由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列vk=Akv0，从而计算主特征值1及其对应的特征向量.,的收敛速度由比值,来确定，r越小收敛越快，但当r1时收敛可能很慢.,24,定理7设ARnn有n个线性无关的特征向量,主特征值1满足,|1|2|n|,则对任何非零向量v0,又设A有n个线性无关的特征向量,1对应的r个线性无关的特征向量为x1,x2,xr,则,如果A的主特征值为实的r重根,即1=2=r,且,|r|r+1|n|，,25,为1的特征向量,这说明当A的主特征值是实的重根时,定理7的结论还是正确的.,26,应用幂法计算A的主特征值1及其对应的特征向量时,如果|1|1(或|1|2|n|,则对任意非零初始向量v0=u0(a10),有幂法计算公式为(uk,vk),则有,32,例2用幂法计算矩阵,的主特征值和相应的特征向量.,解取v0=u0=(1,1,1)T,则,33,计算结果如下表,34,这个结果是用8位浮点数字进行运算得到的,uk的分量值是舍入值.于是得到,及相应的特征向量(0.7482,0.6497,1)T.1和相应的特征向量的真值(8位数字)为,35,例用幂法计算矩阵,的主特征值与其对应的特征向量.,解取v0=u0=(0,0,1)T,则,36,直到k=8时的计算结果见下表,从而,37,8.2.2加速方法,1、原点平移法,由前面讨论知道，应用幂法计算A的主特征值的收敛速度主要由比值r=|2/1|来决定，但当r接近于1时，收敛可能很慢.这时，一个补救办法是采用加速收敛的方法.,其中p为参数，设A的特征值为i，则对矩阵B的特征值为i-p，而且A,B的特征向量相同.,引进矩阵B=A-pI.,38,如果要计算A的主特征值1,只要选择合适的数p，使1-p为矩阵B=A-pI的主特征值，且,那么，对矩阵B=A-pI应用幂法求其主特征值1-p,收敛速度将会加快.这种通过求B=A-pI的主特征值和特征向量，而得到A的主特征值和特征向量的方法叫原点平移法.对于A的特征值的某种分布，它是十分有效的.,39,例3设AR44有特征值,比值r=|2/1|0.9.做变换,B=A-12I(p=12),则B的特征值为,应用幂法计算B的主特征值1的收敛速度的比值为,虽然常常能够选择有利的p值,使幂法得到加速,但设计一个自动选择适当参数p的过程是困难的.,40,下面考虑当A的特征值是实数时，怎样选择p使采用幂法计算1得到加速.,且使收敛速度的比值,设A的特征值都是实数，且满足,则对实数p，使矩阵A-pI的主特征值为1-p或n-p时，当我们计算1及x1时，首先应选取p使,41,显然，当2-p=-(n-p)时，即P=(2+n)/2P*时为最小值，这时收敛速度的比值为,当希望计算n时，应选取p=(1+n-1)/2P*使得应用幂法计算n得到加速.,当A的特征值都是实数，满足(2.10)式,且2,n能初步估计出来，我们就能确定P*的近似值.,42,例4用原点平移加速法求例2中矩阵A的主特征值与其对应的特征向量.,对B应用幂法，仍取v0=(1,1,1)T,则,解取p=0.75,做平移变换B=A-pI，则,43,计算结果如下表,由此得B的主特征值为11.7865914，A的主特征值1为1=1+0.752.5365914.,44,原点位移的加速方法，是一个矩阵变换方法.这种变换容易计算，又不破坏矩阵A的稀疏性，但p的选择依赖对A的特征值分布的大致了解.,及相应的特征向量为x1(0.7482,0.6497,1)T.,与例2结果比较，上述结果比例2迭代15次还好.若迭代15次11.7865258(相应的12.5365258)就是真值.,45,设ARnn为对称矩阵，称,为向量x的瑞利商，其中(x,x)=xTx为内积.由定理4知道，实对称矩阵A的特征值1及n可用瑞利商的极限值表示.下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征值的加速上来.,2、瑞利商(Rayleigh)加速,46,定理9设ARnn为对称矩阵，特征值满足,对应的特征向量xi满足(xi,xj)=ij(单位正交向量),应用幂法公式(2.9)计算A的主特征值1，则规范化向量uk的瑞利商给出1的较好的近似值为,由此可见，R(uk)比k更快的收敛于1.,47,证明由(2.8)式及,得,48,幂法的瑞利商加速迭代公式可以写为,其中A为n阶实对称矩阵.,对给定的误差限，当|kk-1|0，则对任意非零初始向量u0(an0)，由反幂法计算公式构造的向量序列vk，uk满足,(1),(2),收敛速度的比值r=|n/n-1|.,54,在反幂法中也可以用原点平移法加速迭代过程，或求其它特征值与其对应的特征向量.,如果矩阵(A-pI)-1存在，显然其特征值为,对应的特征向量仍然是x1,x2,xn，现对矩阵(A-pI)-1应用幂法，得到反幂法的迭代公式,55,如果p是A的特征值j的一个近似值，且设j与其它特征值是分离的，即,就是说1/(j-p)是矩阵(A-pI)-1的主特征值，可用反幂法(2.12)计算特征值j及特征向量.,设ARnn有n个线性无关的特征向量x1,x2,xn,则,56,其中,同理可得下面的定理.,57,定理11设ARnn有n个线性无关的特征向量，矩阵A的特征值及对应的特征向量分别记为i及xi(i=1,2,n)，而p为j的近似值，(A-pI)-1存在，且,(1),(2),则对任意非零初始向量u0(aj0)，由反幂法计算公式(2.12)构造的向量序列vk，uk满足,且收敛速度为,58,由该定理知，对A-pI(其中pj)应用反幂法，可用来计算特征向量xj，只要选择p是j的一个较好的近似且特征值分离情况较好，一般r很小，常常只要迭代一二次就可完成特征向量的计算.,反幂法迭代公式中的vk是通过解方程组,求得的,为了节省工作量,可以先将A-pI进行三角分解,于是求vk相对于解两个三角形方程组,59,实验表明,按下述方法选择u0是较好的:选u0使,用回代求解三角形方程组(2.13)即得v1，然后再按公式(2.12)进行迭代.,60,反幂法计算公式：,1.分解计算P(A-pI)=LU,且保留L,U及P信息,2.反幂迭代法,(1)解Uv1=(1,1,1)T求v1,(2)解k=2,3,解Lyk=Puk-1求yk,解Uvk=yk求vk,k=max(vk),计算uk=vk/k,61,例5用反幂法求,的对应于计算特征值=1.2679(精确特征值为)的特征向量(用5位浮点数进行计算).,解用部分选主元的三角分解将A-pI(其中p=1.2679)分解为,P(A-pI)=LU,其中,62,由Uv1=(1,1,1)T,得,由LUv2=Pu1,得,3对应的特征向量是,由此看出u2是x3的相当好的近似.,特征值,3的真值为,63,8.3正交变换与矩阵分解,8.3.1豪斯霍尔德变换,定义设向量wRn，且wTw=1，称矩阵,为初等反射矩阵，这个矩阵也称为豪斯霍尔德变换.如果记w=(w1,w2,wn)，则有,64,定理12设有初等反射矩阵，其中H=I-2wwT，其中wTw=1，则,(1)H是对称矩阵，即HT=H.,(2)H是正交矩阵，即H-1=H.,(3)设A为对称矩阵，那么A1=H-1AH=HAH即亦是对称矩阵.,证明只证H的正交性，其中显然.,HTH=H2=(I-2wwT)(I-2wwT)=I-4wwT+4w(wTw)wT=I.,设向量u0，则显然,是一个初等反射矩阵.,65,初等反射矩阵的几何意义:考虑以为w法向量且过原点O的超平面S:wTx=0.任取向量vRn,则v=x+y,其中xS,yST.于是,Hx=(I-2wwT)x=x-2wwTx=x.,从而,Hv=x-y=v.,其中v为v关于平面S的镜面反射.,初等反射矩阵在计算上的意义是它能用来约化矩阵.设向量Hx0,可选择一初等反射矩阵H,使Hx=y.,Hy=(I-2wwT)y=y-2wwTy=-y.,66,定理13设x,y为两个不相等的n维向量,|x|2=|y|2,则存在一个初等反射矩阵H，使Hx=y.,证明令,则得到一个初等反射矩阵,而且,因为,所以,w是使Hx=y成立的唯一长度等于1的向量(不计符号).,67,定理14(约化定理)设x=(x1,x2,xn)T0,则存在初等反射矩阵H，使Hx=-e1，其中,证明记y=-e1,设xy,则有|x|2=|y|2,由定理13,存在变换H=I-2wwT,使Hx=y=-e1,其中,68,其中记u=(x1+,x2,xn)T，.,如果和x1异号,那么计算x1+时有效数字可能损失,我们取和x1有相同的符号,即取,在计算时,可能上溢或下溢,为了避免溢出,将x规范化,记u=x+e1=(u1,u2,un)T，于是,69,则有H使Hx=e1，其中,例6设x=(3,5,1,1)T，则|x|2=6.取=6,可直接验证Hx=-e1=(-6,0,0,0)T.,70,8.3.2吉文斯变换,设向量x,yR2，则变换,或,是平面上向量的一个旋转变换，其中,为正交变换.,R2中变换,其中x=(x1,x2,xn)T,y=(y1,y2,yn)T，而,71,称为Rn中平面xi,xj的旋转变换，也称为吉文斯变换.P=P(i,j,)=P(i,j)称为平面旋转矩阵.,72,显然，P=P(i,j,)具有性质：,(1)P与单位阵I只是在(i,i),(i,j),(j,i),(j,j)位置元素不一样，其它相同.,(2)P为正交矩阵(P-1=PT).,(3)P(i,j)A(左乘)只需计算第i行与第j行元素，即对A=(aij)mn有,其中c=cos,s=sin.,(4)AP(i,j)(右乘)只需计算第i列与第j列元素,利用平面旋转变换，可得向量x中的指定元素变为零.,73,定理15(约化定理)设x=(x1,xi,xj,xn)T,其中不全为零，则可选择平面旋转阵P=P(i,j,)，使,其中,于是，由c,s的取法得,证明取c=cos=xi/xi,s=sin=xj/xi，由P(i,j,)x=x=(x1,xi,xj,xn)T,利用矩阵乘法,显然有,74,8.3.3矩阵的QR分解与舒尔分解,定理16设ARnn非奇异，则存在正交矩阵P，使PA=R，其中R为上三角矩阵.,证明1.用吉文斯变换给出构造P的方法.,(1)第1步约化，由设有j(j=1,2,n)使aj10，则可选择吉文斯变换P(1,j)，将aj1处的元素化为零.即若aj10(j=2,3,n)，则存在P(1,j)使得,可简记为P1A=A(2)，其中P1=P(1,n)P(1,2).,75,(2)第k步约化，设上述过程已完成第1步至第k-1步，于是有,由设有j(njk)使ajk(k)0(j=k+1,n)，则可选择吉文斯变换P(k,j)(j=k+1,n)使,PkA(k)=P(k,n)P(k,k+1)A(k)=PkPk-1P1A=A(k+1),其中Pk=P(k,n)P(k,k+1).,76,(3)继续上述约化过程，最后则有,Pn-1P2P1A=A(n)=R(上三角形矩阵).,令P=Pn-1P2P1,它是一个正交矩阵,有PA=R.,2.用豪斯霍尔德变换构造出正交矩阵P.记A(0)=A，它的第一列记为，不妨设a1(0)0，找到矩阵,于是,其中,77,一般地，设,其中D(j-1)为j-1阶方阵，其对角线以下元素均为0，矩阵A(j-1)为n-j+1阶方阵，设其第一列为a1(j-1)，可选择n-j+1的豪斯霍尔德矩阵变换,使,根据Hj构造nn阶的变换矩阵Hj为,78,于是有,它和A(j-1)有类似的形式，只是D(j)为j阶方阵，其对角线以下元素是0，这样经过n-1步运算得到,Hn-1H2H1A=A(n-1)=R(上三角形矩阵).,其中P=Hn-1H2H1是一个正交矩阵,从而有PA=R.,79,定理17(QR分解定理)设ARnn为非奇异矩阵,则存在正交矩阵Q与上三角矩阵R，使,A=QR，,且当R的对角元素为正时，分解是唯一的.,证明从定理16可知,只要令Q=PT就有A=QR.下面证明分解的唯一性，设有两种分解,A=Q1R1=Q2R2，,其中Q1,Q2为正交矩阵,R1,R2为对角元素均为正的上三角矩阵,则,由假设及对称正定矩阵ATA的楚列斯基(平方根法)分解的唯一性，则得R1=R2，从而可得Q1=Q2.证毕.,80,定理16保证了A可分解为A=QR.若A非奇异，则R也非奇异.如果不规定R的对角元素为正，则分解不是唯一的.一般按吉文斯变换或豪斯霍尔德变换方法作出的分解A=QR，R的对角元不一定是正的，设上三角矩阵R=(rij)，只要令,则为正交矩阵，为对角元是|rii|的上三角矩阵，这样便是符合定理17的唯一QR分解.,81,例7用豪斯霍尔德变换作矩阵A的QR分解：,解按(3.2)式找豪斯霍尔德矩阵H1R33，使,则有,82,再找H2R22,使H2(1.44949,3.44949)T=(*,0)T,得,这是一个上三角矩阵，但对角元素皆为负数.,83,只要令D=-I，则有R=-H2H1A是对角元素为正的上三角矩阵，取,则得A=QR.,除了QR分解,矩阵的Schur分解也是重要的工具,它解决矩阵ARnn可约化到什么程度的问题，对复矩阵ACnn，则存在酉矩阵U,使UHAU为一个上三角矩阵R，其对角线元素就是A的特征值,A=URUH称为A的舒尔分解，对于实矩阵A，其特征值可能有复数，A不能用正交相似变换约化为上三角矩阵，但它可约化为以下形式.,84,其中对角块Rii(i=1,2,m)为一阶或二阶方阵,且每个一阶Rii是A的实特征值，每个二阶对角块Rii的两个特征值是A的两个共轭复特征值.,定理18(实Schur分解)设ARnn，则存在正交矩阵Q使,A=QRQT称为A的实舒尔分解,可以通过逐次正交变换使A趋于实舒尔型矩阵，以求A的特征值.,85,8.3.4用正交相似变换约化一般实矩阵为上海森伯格矩阵,设ARnn，下面来说明,可选择初等反射矩阵U1,U2,Un-2使A经正交相似变换约化为一个上海森伯格阵.,(1)设,86,其中c1=(a21,an1)TRn-1，不妨设c10,否则这一步不需要约化.于是,可选择初等反射阵使，其中,令,87,则,其中,88,(2)第k步约化：重复上述过程，设对A已完成第1步,第k-1步正交相似变换，即有,或,且,89,其中为k阶上海森伯格阵，,设ck0,是可选择初等反射阵Rk使其中，Rk计算公式为,90,令,则,91,其中为k+1阶上海森伯格阵，第k步约化只需计算及(当A为对称矩阵时，只需要计算).,(3)重复上述过程，则有,92,定理19(豪斯霍尔德约化矩阵为上海森伯格阵)设ARnn则存在初等反射矩阵U1,U2,Un-2使,为上海森伯格矩阵.,总结上述结论，有,本算法约需要5n3/3次乘法运算.如果要把U0也算出来还需要增加2n3/3次乘法.,93,例8用豪斯霍尔德方法将矩阵,约化为上海森伯格阵.,解选取初等反射阵R1使,其中c1=(2,4)T.,(1)计算R1:,94,则有,(2)约化计算:令,则得到上海森伯格阵为,95,如果A是对称的，则H=U0TAU0也对称，这时H是一个对称三对角矩阵.,定理20(豪斯霍尔德约化对称矩阵为对称三对角矩阵)设ARnn为对称矩阵，则存在初等反射矩阵U1,U2,Un-2使,为对称三对角矩阵.,96,证明由定理19,存在初等反射矩阵U1,U2,Un-2使为上海森伯格阵，且An-1亦是对称矩阵，因此，An-1为对称三对角矩阵.,由上面讨论可知，当A为对称阵时，由AkAk+1=AkUkAk一步约化计算中只需计算Rk及RkA22(k)Rk.又由于A的对称性，故只需计算RkA22(k)Rk的对角线以下元素.注意到,97,引进记号,则有,对对称阵A用初等反射阵正交相似约化为对角三对角阵大约需要2n3/3次乘法.,用正交矩阵进行相似约化有一些特点，如构造的Uk容易求逆，且Uk的元素数量级不大，这个算法是十分稳定的.,98,8.4QR方法,8.4.1QR算法,QR方法是一种