数值分析第8章矩阵特征值问题计算.ppt

第8章 矩阵特征问题的计算,8.1 引言8.2 幂法及反幂法8.3 豪斯霍尔德方法8.4 QR方法,8.1 引 言,工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题.,下面先复习一些矩阵的特征值和特征向量的基础知识.,定义1 已知n阶矩阵A=(aij)，则,称为A的特征多项式.,一般有n个根(实的或复的，复根按重数计算)称为A的特征值.用(A)表示A的所有特征值的集合.,A的特征方程,设为A的特征值，相应的齐次方程组,注：当A为实矩阵时，()=0为实系数n次代数方程，其复根是共轭成对出现.,的非零解x称为矩阵A的对应于的特征向量.,例1 求A的特征值及特征向量，其中,解 矩阵A的特征方程为,求得矩阵A的特征值为：,对应于各特征值矩阵A的特征向量分别为：,定理1 设为ARnn的特征值,且Ax=x(x0)，则有,-p为A-pI的特征值，即(A-pI)x=(-p)x；,c为的cA特征值(c0为常数)；,下面叙述有关特征值的一些结论：,k为Ak的特征值，即Akx=kx；,设A为非奇异矩阵，那么0,且-1为A-1的特征值，即A-1x=-1x.,定理2 设i(i=1,2,n)为n阶矩阵A=(aij)的特征值，则有,称为A的迹；,定理3 设ARnn，则有,定理4 设A 为分块上三角矩阵，即,其中每个对角块Aii均为方阵，则,定理5 设A与B为相似矩阵（即存在非奇异矩阵P使B=P-1AP），则,定理5说明，一个矩阵A经过相似变换，其特征值不变.,一个亏损矩阵是一个没有足够特征向量的矩阵，亏损矩阵在理论上和计算上都存在困难.,A与B有相同的特征值；,如果y是B的特征向量，则Py是A的特征向量.,定义2 如果实矩阵A有一个重数为k的特征值,且对应于的A的线性无关的特征向量个数 k，则A称为亏损矩阵.,定理6 ARnn可对角化，即存在非奇异矩阵P使,的充分必要条件是A具有n个线性无关的特征向量.,如果ARnn有 m个(mn)不同的特征值1,2,m，则对应的特征向量 x1,x2,xm 线性无关.,定理7(对称矩阵的正交约化)设ARnn为对称矩阵，则,存在一个正交矩阵P使的,且1,2,n为A的特征值，而P(u1,u2,un)列向量uj为A的对应于j 的单位特征向量.,A的特征值均为实数；,A有n个线性无关的特征向量；,定义3 设n阶矩阵A=(aij)，令,下面讨论矩阵特征值界的估计.,；,集合 称为复平面上以aii为圆心，以ri为半径的n阶矩阵A的n个Gerschgorin圆盘.,定理8(Gerschgorin圆盘定理),特别地，如果A的一个圆盘Di是与其它圆盘分离(即孤立圆盘)，则Di中精确地包含A的一个特征值.,设n阶矩阵A(aij)，则A的每一个特征值必属于下面某个圆盘之中,如果A有m个圆盘组成一个连通的并集S，且S与余下n-m个圆盘是分离的，则S内恰包含A的m个特征值.,或者说 A的特征值都在n个圆盘的并集中.,证明 只就给出证明.设为A的特征值，即,Ax=x，其中x=(x1,x2,xn)T0.,或,记，考虑Ax=x的第k个方程，即,于是,即,这说明，A的每一个特征值必位于A的一个圆盘中，并且相应的特征值一定位于第k个圆盘中(其中k是对应特征向量x绝对值最大的分量的下标).,利用相似矩阵性质，有时可以获得A的特征值进一步的估计，即适当选取非奇异对角阵,并做相似变换.适当选取 可使某些圆盘半径及连通性发生变化.,例2 估计矩阵A的特征值范围，其中,解 矩阵A的3个圆盘为,由定理8，可知A的3个特征值位于3个圆盘的并集中，由于D1是孤立圆盘，所以D1内恰好包含A的一个特征值1(为实特征值)，即,A的其它两个特征值2,3包含在D2,D3的并集中.,现在取对角阵,做相似变换,矩阵A1的3个圆盘为,显然，3个圆盘都是孤立圆盘，所以，每一个圆盘都包含A的一个特征值(为实特征值)，且有估计,当A为实矩阵，如果限制用正交相似变换，由于A有复的特征值，A不能用正交相似变换约化为上三角阵.用正交相似变换能约化到什么程度呢？,定理9(Schur定理)设ARnn，则存在酉矩阵U使,其中rii(i=1,2,n)为 A的特征值.,下面给出理论上有关通过酉相似变换及正交变换可以约化一般矩阵A到什么程度的问题.,其中Rii(i=1,2,m)为一阶或二阶方阵，且每个一阶Rii是A的实特征值，每个二阶对角块Rii的两个特征值是 A的两个共轭复特征值.,定理10(实Schur分解)设ARnn，则存在正交矩阵Q使,定义4 设ARnn为对称矩阵，对于任一非零向量x，称,我们转向实Schur型的实际计算.,为对应于向量x的瑞利(Rayleigh)商.,定理11 设ARnn为对称矩阵(其特征值次序记为12n)，则,1.(对任何非零xRn)；,2.；,3.,证明 只证1，关于2,3自己作练习.,由于A为实对称矩阵，可将 1,2,n 对应的特征向量 x1,x2,xn 正交规范化，则有(xi,xj)=ij，设x0为Rn中任一向量，则有,于是,从而1成立.结论1说明瑞利商必位于n和1之间.,关于计算矩阵A的特征值问题，当n2,3时，我们还可按行列式展开的办法求()=0的根.但当n较大时，如果按展开行列式的办法，首先求出()的系数，再求()的根，工作量就非常大，用这种办法求矩阵的特征值是不切实际的，由此需要研究求A的特征值及特征向量的数值解法.,本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法（迭代法），另一类是正交相似变换的方法（变换法）.,幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法，所不同的是幂法是计算矩阵的主特征值(矩阵按模最大的特征值称为主特征值，其模就是该矩阵的谱半径)和相应特征向量的一种向量迭代法，而反幂法则是计算非奇异(可逆)矩阵按模最小的特征值和相应特征向量的一种向量迭代法.下面分别介绍幂法与反幂法.,8.2 幂法及反幂法,现讨论求1及x1的方法.,设实矩阵A=(aij)有一个完全的特征向量组，即A有n个线性无关的特征向量，设矩阵A的特征值为1,2,n,相应的特征向量为x1,x2,xn.已知A的主特征值1是实根，且满足条件,8.2.1 幂法（又称乘幂法）,幂法的基本思想是:任取非零的初始向量v0,由矩阵A构造一向量序列vk,称为迭代向量，由假设，v0可唯一表示为,于是,其中,由假设 故 从而,为1的特征向量.,所以当k充分大时，有,即为矩阵A的对应特征值1 的一个近似特征向量.,用(vk)i 表示vk的第i个分量，则当k充分大时，有,即为A的主特征值1的近似值.,由于,这种由已知非零向量v0及矩阵A的乘幂Ak构造向量序列vk以计算A的主特征值1(2.7)及相应特征向量(2.5)的方法就称为幂法.,迭代公式实质上是由矩阵A的乘幂 Ak与非零向量v0相乘来构造向量序列vk=Akv0，从而计算主特征值1及其对应的特征向量，这就是幂法的思想.,的收敛速度由比值,来确定，r越小收敛越快，但当r1时收敛可能很慢.,定理12 设ARnn有n个线性无关的特征向量，主特征值1满足条件,|1|2|n|，,则对任何非零向量v0(a10)，幂法的算式成立.,又设A有n个线性无关的特征向量，1对应的r个线性无关的特征向量为x1,x2,xr，则由(2.2)式有,如果A的主特征值为实的重根,即1=2=r,且,|r|r+1|n|，,为A的特征向量，这说明当A的主特征值是实的重根时，定理5的结论还是正确的.,应用幂法计算A的主特征值1及其对应的特征向量时，如果|1|1(或|1|1)，迭代向量 vk的各个不等于零的分量将随 k 而趋向于无穷(或趋向于零)，这样在计算机实现时就可能“溢出”.为克服这个缺点，就需要将迭代向量加以规范化.,设有一向量v0，将其规范化得向量为,其中max(v)表示v的绝对值最大的分量.即如果有,则max(v)=vq，且q为所有绝对值最大的分量中的最小下标.,在定理12的条件下幂法可这样进行：任取一初始向量v00(a10)，构造规范化向量序列为,由(2.3)式,收敛速度由比值r=|2/1|确定.总结上述结论，有,同理，可得到,定理13 设ARnn有n个线性无关的特征向量，主特征值1满足|1|2|n|，则对任意非零初始向量v0=u0(a10)，有幂法计算公式为,例1 用幂法计算矩阵,的主特征值与其对应的特征向量.,解取 v0=u0=(0,0,1)T,则,直到k=8 时的计算结果见下表,从而,见书p303-例3.,8.2.2 幂法的加速方法,1、原点平移法,由前面讨论知道，应用幂法计算A的主特征值的收敛速度主要由比值 r=|2/1|来决定，但当r 接近于1时，收敛可能很慢.这时，一个补救办法是采用加速收敛的方法.,其中p为参数，设A的特征值为i，则对矩阵B的特征值为i-p，而且A,B的特征向量相同.,引进矩阵 B=A-pI.,如果要计算A的主特征值1,只要选择合适的数p，使1-p为矩阵B=A-pI 的主特征值，且,那么，对矩阵B=A-pI应用幂法求其主特征值1-p,收敛速度将会加快.这种通过求B=A-pI的主特征值和特征向量，而得到A的主特征值和特征向量的方法叫原点平移法.对于A的特征值的某种分布，它是十分有效的.,例4 设AR44有特征值,比值r=|2/1|0.9.做变换,B=A-12I(p=12),则B的特征值为,应用幂法计算B的主特征值1的收敛速度的比值为,虽然常常能够选择有利的p值,使幂法得到加速,但设计一个自动选择适当参数p的过程是困难的.,下面考虑当A的特征值是实数时，怎样选择p使采用幂法计算1得到加速.,且使收敛速度的比值,设A的特征值都是实数，且满足,则对实数p，使矩阵A-pI的主特征值为1-p或n-p时，当我们计算1及x1时，首先应选取p使,显然，当2-p=-(n-p)时，即 P=(2+n)/2P*时为最小值，这时收敛速度的比值为,当希望计算n时，应选取 p=(1+n-1)/2P*使得应用幂法计算n得到加速.,当A的特征值都是实数，满足,且2,n能初步估计出来，我们就能确定P*的近似值.,例2 用原点平移加速法求例1中矩阵A的主特征值与其对应的特征向量.,对B应用幂法，仍取 v0=(0,0,1)T,则,解 取p=-2.5,做平移变换B=A-pI，则,迭代5步的计算结果见下表,可得到B的主特征值为 113.5000,主特征向量为 v1(0.5,1.0,0.7500)T,因此，A的主特征值为 1=1+p 11.0000,主特征向量仍为 x1=(0.5,1,0.7500)T.,原点位移的加速方法，是一个矩阵变换方法.这种变换容易计算，又不破坏矩阵A的稀疏性，但p的选择依赖对A的特征值分布的大致了解.,见书p306-例5.,设ARnn为对称矩阵，称,为向量x的瑞利商，其中(x,x)=xTx为内积.由定理11知道，实对称矩阵A的特征值1及n可用瑞利商的极限值表示.下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征值的加速上来.,2、瑞利商(Rayleigh)加速,定理14 设ARnn为对称矩阵，特征值满足,对应的特征向量xi满足(xi,xj)=ij(单位正交向量)，应用幂法公式(2.9)计算A的主特征值1，则规范化向量uk的瑞利商给出1的较好的近似值为,由此可见，R(uk)比k更快的收敛于1.,证明 由(2.8)式及,得,幂法的瑞利商加速迭代公式可以写为,其中A为n阶实对称矩阵.,对给定的误差限，当|kk-1|时，取近似值,8.2.3 反幂法,反幂法是用于求非奇异矩阵A的按模最小的特征值和对应特征向量的方法.而结合原点平移法的反幂法则可以求矩阵A的任何一个具有先了解的特征值和对应的特征向量。,设矩阵A非奇异,其特征值i(i=1,2,n),满足,其相应的特征向量x1,x2,xn线性无关，则 A-1 的特征值为1/i，对应的特征向量仍为xi(i=1,2,n).,此时,A-1的特征值满足,因此,对A-1应用幂法,可求出其主特征值1/n k 和特征向量 xn uk.,从而求得A的按模最小特征值 n 1/k 和对应的特征向量 xn uk，这种求A-1的方法就称为反幂法.,为了避免求A-1,可通过解线性方程组Avk=uk-1得到vk，采用LU分解法，即先对A进行LU分解A=LU,此时反幂法的迭代公式为,反幂法的迭代公式为,对给定的误差，当|kk-1|时，得,显然，反幂法的收敛速度取决于比值，比值越小，收敛越快.,定理15 设ARnn为非奇异矩阵，且有n个线性无关的特征向量，其对应的特征值满足|1|2|n-2|n|0，则对任意非零初始向量u0(an0)，由反幂法计算公式构造的向量序列vk，uk满足,在反幂法中也可以用原点平移法加速迭代过程，或求其它特征值与其对应的特征向量.,如果矩阵(A-pI)-1存在，显然其特征值为,对应的特征向量仍然是x1,x2,xn，现对矩阵(A-pI)-1应用幂法，得到反幂法的迭代公式,如果p是A的特征值j的一个近似值，且设j与其它特征值是分离的，即,就是说1/(j-p)是矩阵(A-pI)-1的主特征值，可用反幂法(2.12)计算特征值及特征向量.,设ARnn有 n个线性无关的特征向量 x1,x2,xn，则,其中,同理可得：,定理16 设ARnn有n个线性无关的特征向量，矩阵A的特征值及对应的特征向量分别记为i 及xi(i=1,2,n)，而p为j的近似值，(A-pI)-1存在，且,则对任意非零初始向量u0(aj0)，由反幂法计算公式(2.12)构造的向量序列vk，uk满足,且收敛速度为,由该定理知，对A-pI(其中pj)应用反幂法，可用来计算特征向量xj，只要选择p是j的一个较好的近似且特征值分离情况较好，一般r很小，常常只要迭代一二次就可完成特征向量的计算.,反幂法迭代公式中的vk是通过解方程组,求得的,为了节省工作量,可以先将A-pI进行三角分解,于是求vk相对于解两个三角形方程组,实验表明,按下述方法选择u0是较好的:选u0使,用回代求解(2.13)即得v1,然后再按公式(2.12)进行迭代.,反幂法计算公式见书p311.前面已提到.,见书p311-例6.,8.3 豪斯霍尔德方法,(1)用初等反射矩阵作正交相似变换约化一般实矩阵A为上海森伯格矩阵.,8.3.1 引言,本节讨论两个问题,(2)用初等反射矩阵作正交相似变换约化对称矩阵A为对称三对角矩阵.,于是，求原矩阵特征值问题，就转化为求上海森伯格矩阵或对称三对角矩阵的特征值问题.,8.3.2 用正交相似变换 约化一般实矩阵为上海森伯格矩阵,设ARnn，下面来说明，可选择初等反射矩阵U1,U2,Un-2使A经正交相似变换约化为一个上海森伯格阵.,(1)设,其中c1=(a21,an1)TRn-1,不妨设c10，否则这一步不需要约化.于是,可选择初等反射阵 使，其中,令,则,其中,(2)第k步约化：重复上述过程，设对A已完成第1步,第k-1步正交相似变换，即有,或,且,其中 为k阶上海森伯格阵，,设ck0,于是可选择初等反射阵Rk使 其中，Rk计算公式为,令,则,其中 为k+1阶上海森伯格阵，第k步约化只需计算 及(当A为对称矩阵时，只需要计算).,(3)重复上述过程，则有,定理17(豪斯霍尔德约化矩阵为上海森伯格阵)设ARnn则存在初等反射矩阵U1,U2,Un-2 使,为上海森伯格矩阵.,总结上述结论，有,算法1(豪斯霍尔德约化矩阵为上海森伯格阵)设ARnn，本算法计算U0TAU0=H(上海森伯格型)，其中U0=U1U2Un-2为初等反射阵的乘积.,1.U0I,2.对于k=1,2,n-2,(1)计算初等反射阵Rk使,本算法约需要5n3/3次乘法运算，要明显形成U0还需要附加2n3/3次乘法.,(2)约化计算,(3)U0U0Uk,例7 用豪斯霍尔德方法将矩阵,约化为上海森伯格阵.,解 选取初等反射阵R1使,其中c1=(2,4)T.,(1)计算R1:,则有,(2)约化计算:,则得到上海森伯格阵为,8.3.3 用正交相似变换 约化对称矩阵为对称三对角矩阵,定理18(豪斯霍尔德约化对称矩阵为对称三对角矩阵)设ARnn为对称矩阵，则存在初等反射矩阵U1,U2,Un-2使,为对称三对角矩阵.,证明 由定理17,存在初等反射矩阵U1,U2,Un-2 使 为上海森伯格阵，且An-1亦是对称阵，因此，An-1为对称三对角阵.,由上面讨论可知，当A为对称阵时，由AkAk+1=Ak Uk Ak一步约化计算中只需计算Rk及Rk A22(k)Rk.又由于A的对称性，故只需计算Rk A22(k)Rk的对角线以下元素.注意到,引进记号,则有,对对称阵A用初等反射阵正交相似约化为对角三对角阵大约需要2n3/3次乘法.,用正交矩阵进行相似约化有一些特点，如构造的，Uk容易求逆，且Uk的元素数量级不大，这个算法是十分稳定的.,算法2见书p318.,8.4 QR 方 法,8.4.1 QR算法,Rutishauser(1958)利用矩阵的三角分解提出了计算矩阵特征值的LR算法，Francis(1961,1962)利用矩阵的QR分解建立了计算矩阵特征值的QR算法.,QR方法是一种变换方法，是计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一.,(1)上海森伯格矩阵的全部特征值问题；,(2)计算对称三对角矩阵的全部特征值问题，,目前QR方法主要用来计算：,且QR方法具有收敛快，算法稳定等特点.,对于一般矩阵ARnn(或对称矩阵)，则首先用豪斯霍尔德方法将A化为上海森伯格阵B(或对称三对角阵)，然后再用QR方法计算B的全部特征值.,设ARnn，且A对进行QR分解，即,A=QR，,其中R为上三角阵,Q为正交阵,于是可得到一新矩阵,B=RQ=QTAQ.,显然，B是由A经过正交相似变换得到，因此B与A的特征值相同.再对B进行QR分解，又可得一新矩阵,重复这一过程可得到矩阵序列：,设A=A1,将A1进行QR分解A1=Q1R1,作矩阵A2=R1Q1=Q1TR1Q1,QR算法，就是利用矩阵的QR分解，按上述递推法则构造矩阵序列Ak的过程.只要A为非奇异矩阵，则由QR算法就完全确定Ak.,定理19(基本QR方法)设A=A1Rnn,构造QR算法:,证明(1),(2)显然，现证(3).用归纳法，显然，当k=1时有,设Ak-1有分解式,于是,由第5章定理30或定理31知，将Ak进行QR分解,即将Ak用正交变换(左变换)化为上三角矩阵,这就是说Ak+1可由Ak按下述方法求得：,其中QkT=Pn-1P2P1，故,(1)左变换Pn-1P2P1Ak=Rk(上三角阵);,(2)右变换RkP1TP2TPn-1T=Ak+1.,定理20(QR方法的收敛性)设A=(aij)Rnn,(1)如果A的特征值满足:,(2)A有标准型A=XDX-1,其中D=diag(1,2,n),且设X-1有三角分解X-1=LU(L为单位下三角阵，U为上三角阵)，则由QR算法产生的Ak本质上收敛于上三角矩阵，即,若记Ak=(aij(k)，则,(1),(2)当ij时，,当ij时aij(k)的极限不一定存在.,证明可参阅1.,定理21 如果对称矩阵A满足定理20的条件,则由QR算法产生的Ak收敛于对角阵D=diag(1,2,n).,证明 由定理20即知.,设ARnn，且A有完备的特征向量集合，如果A的等模特征值中只有实重特征值或多重共轭复特征值，则由QR算法产生的Ak本质收敛于分块上三角矩阵(对角块为一阶和二阶子块)，且对角块中每一个22子块给出A的一对共轭复特征值，每一个一阶对角子块给出A的实特征值，即,关于QR算法收敛性的进一步结果为：,其中m+2l=n,Bi为22子块,它给出A一对共轭复特征值.,8.4.2 带原点位移的QR方法,p322-自学.,8.4.3 用单步QR方法计算上海森伯格阵特征值,p325-自学.,