矩阵特征值问题的解法.ppt

1,Numerical eigenvalue of matrix矩阵特征值问题的解法,2,给出.若有 使得：,则称 为矩阵 的特征值,为相应的特征向量。特征值 为特征方程的根。,3,与矩阵想干的一些重要结果:eigenvalueofmarix.doc,4,特征值的估计与扰动问题,特征值的估计,称之为Gerschgorin圆盘（盖尔圆）.,Gerschgorin 圆盘定理,5,第二圆盘定理,设 为 阶实方阵，如果 的 个Gerschgorin圆盘与其他圆盘不相连，则恰好有 的 个特征值落在该 个圆盘的并集之中，即：,特别地：孤立圆盘仅含有一个特征值.,为 的一个重新排列,则 中含有 的 个特征值.,6,试讨论A的特征值的分布.,解 由A确定的3个圆盘分别为,所以 315-222-63-2,例 1 设矩阵,R1=-41,R2=2,R3=+42,x,y,0,-2,-4,-6,2,3,4,5,实际上，1=4.20308,2=-0.442931,3=-3.76010,7,关于实对称矩阵的极大极小定理,为矩阵 关于向量 的Rayleigh（雷利）商.,为 阶实对称矩阵，则其特征值皆为实数，记做：,并且存在规范正交特征向量系，满足：,8,由于，对于任意，可以取，使得：.,证明:假设 为 的规范正交特征向量组，则对任何向量，有,9,于是,因而,特别地，若取，这时,从而.同理可证,10,按模最大特征值和特征向量的乘幂法,设A是n阶矩阵，其n个特征值按模从大到小排序为,又假设关于1，2，n的特征向量v1，v2，vn线性无关.,11,任意取定初始向量x0,.,建立迭代公式：,12,因为,故当k时，xk1ka1v1.,因此，xk可看成是关于特征值1的近似特征向量,有一严重缺点，当|1|1（或|1|1时）Vk中不为零的分量将随K的增大而无限增大，计算机就可能出现上溢（或随K的增大而很快出现下溢）,13,在实际计算时，须按规范法计算，每步先对向量xk进行“规范化”。迭代格式改为,14,对任意给定的初始向量x0,类似地,15,当10时,当10时,16,按模最大特征值1及其相应的特征向量v1的乘幂法的计算公式：,17,定理,设 ARnn为非亏损矩阵，其主特征根 1为实根，且|1|2|n|。则从任意非零向量(满足)出发,迭代 收敛到主特征向量，收敛到1。,每个不同的特征根只对应一个Jordan块,若有 1=2，则此法不收敛。,任取初始向量时，因为不知道，所以不能保证 1 0，故所求得之 不一定是，而是使得 的第一个，同时得到的特征根是m。,18,算法:乘幂法 Step 1 Set k=1;Step 2 Find index such that|V0 index|=|V0|;Step 3 Set V0=V0/V0 index;/*规格化 V0*/Step 4 While(k Nmax)do steps 5-11Step 5 V=A V0;/*由Uk1 计算 Vk*/Step 6=V index;Step 7 Find index such that|V index|=|V|;Step 8 If V index=0 then Output(“A has the eigenvalue 0”;V0);STOP./*矩阵是奇异的，用户尝试新的 V0*/Step 9 err=|V0 V/V index|;V0=V/V index;/*计算 Uk*/Step 10 If(err TOL)then Output(;V0);STOP./*成功*/Step 11 Set k+;Step 12 Output(Maximum number of iterations exceeded);STOP./*失败*/,19,求矩阵A的按模最大的特征值,解 取x(0)=(1,0)T,计算x(k)=Ax(k-1),结果如下,例1,可取0.41263,x1(0.017451,0.014190)T.,20,如用规范化乘幂法解例1,仍取u(0)=v(0)=(1,0)T,则有,故可取 1 0.412627,x1(1,0.813138)T.,用乘幂法求A的按模最大的特征值和相应特征向量.,例3 设,解 取初值u(0)=v(0)=(1,1,1)T,计算得,21,可取 1 6.000837,x1(1,0.714316,-0.249895)T.,实际上,A的3个特征值分别为1=6,2=3,3=2.,22,Remark1:具体计算时，U(0)的选取很难保证一定有10。但是，由于舍入误差的影响，只要迭代次数足够多，如，就会有，因而最后结论是成立的。对于 的情形，由于对任意l均有上面的结论，故只要取另外的l使 即可。,Remark2:以上讨论只是说明了乘幂法的基本原理。当 太小或太大时，将会使U(k)分量的绝对值过小或过大，以致运算无法继续进行。因此，实际计算时，常常是每进行m步迭代进行一次规范化，如用,其中，max(U(m)表示向量U(m)的绝对值最大的分量。,23,代替U(k)继续迭代。由于特征向量允许差一个非零常数因子，因而从V(k)往后继续迭代与从U(k)往后继续迭代的收敛速度是相同的，但规范化的做法有效防止了溢出现象。至于m的选取，可以自由掌握，如取m1，5等等。,Remark3:若主特征值是重特征值，如,则有,从而,24,由此可得乘幂法的算法。但是应该注意到，在重特征值的情形下，从不同的非零初始向量出发迭代，可能得到主特征值的几个线性无关的特征向量。,Remark4:由上述推导可知，乘幂法收敛的快慢取决于比值 的大小，该比值越小收敛越快。由此便提出了乘幂法的加速收敛方法，如Rayleigh商加速法、原点平移法等。,25,且|1=|2|3 n,这时,迭代式可写成,则对充分大的k有,Remark5:对于1-2，或1与2共轭等情形，也可类似进行计算。,v(2i)12i(a1x1+a2x2),v(2i+1)12i+1(a1x1-a2x2),26,于是有,x1v(k+1)+1v(k),x2v(k+1)-1v(k),对于规范化的幂法，由于,u(k+2)=v(k+2)/k+2=Au(k+1)/k+2,=Av(k+1)/k+1k+2=A2u(k)/k+1k+2,于是有,x1k+1u(k+1)+1u(k),x2k+1u(k+1)-1u(k),27,的按模最大特征值和相应的特征向量。,例4 用乘幂法求矩阵,解 取初始向量u(0)=v(0)=(1,1,2)T,计算可得,28,29,加速技术,由于,所以,乘幂法收敛速度取决于比值|2/1|,当|2/1|1时,收敛是很慢的.,1.Aitken 加速方法,由(5.2)式可知,x2=13u(13)-1u(12)=(0,0.631924,0.631924)T.,x1=13u(13)+1u(12)=(4.315961,8.631924,8.631924)T,实际上,A的特征值为1=4,2=-4,3=1.,30,可见,序列k线性收敛于1.,会达到加速收敛的目的.,构造Aitken序列,如把Aitken加速方法用于例3,则有,31,2.原点位移法,作矩阵B=A-pE,则B的特征值为mi=i-p(i=1,2,n),而且对应的特征向量相同.,32,则对B应用乘幂法可达到加速收敛的目的。,解 由于A的特征值为1=6,2=3,3=2,故取p=2.5,则B的特征值为m1=3.5,m2=0.5,m3=-0.5,则,如果选取p,使m1仍然是B的按模最大特征值，且满足,取初始向量u(0)=v(0)=(1,1,1)T,由规范化计算公式:,例,用原点位移法求例3中矩阵A的按模最大的特征值和特征向量.,33,计算可得,34,这是因为|2/1|=1/2,而|m2/m1|=1/7,故对B应用乘幂法远比对A应用乘幂法收敛的快.,取,16+2.5=6.000102,x1u(6)=(1,0.714287,0.249995)T,35,反幂法 inverse power method,Q:在每一步我们怎样计算?,A:先对A进行三角分解，再解线性系统.,若知道某一特征根 i 的大致位置 p，即对任意 j i 有|i p|j p|，并且如果(A pI)1存在，则可以用反幂法求(A pI)1的主特征根 1/(i p)，收敛将非常快。,思路,36,也可将上式改写成,式(5.3)称为反幂法.显然有,每一步求v(k)需要求解线性方程组,可采用LU分解法求解.,37,反幂法还可结合原点位移法应用.设已求得矩阵A的特征值i的某个近似值,对B应用反幂法可求出精度更高的i和xi.,设已求得例3中矩阵A的特征值的近似值16.003,和相应的特征向量x1(1,0.714405,-0.249579)T,试用带原点位移的反幂法求1和x1的更精确的值.,作原点位移,令B=A-E,则B的特征值为,例6,解 取p=6.003,作矩阵B=A-6.003E,则,38,取初始向量u(0)=(1,0.714405,-0.249579)T,对B用反幂法计算可得:,可见收敛速度非常快,这是因为B的3个特征值为1=-4.003,2=-3.003,3=-0.003,|3/2|0.000999很小.,39,3 QR方法,QR方法在特征值计算问题的发展上具有里程碑意义。在1955年的时候人们还觉得特征值的计算是十分困扰的问题，到1965年它的计算基于QR方法的程序已经完全成熟。直到今天QR方法仍然是特征值计算的有效方法之一。,40,方法的基本思想,首先作 的 分解：,（对角元非负）,取：然后作 的 分解.,一般地：,于是得矩阵序列：,41,可以证明：,(1),(2),为一阶或二阶方阵.,于是 的特征值即为 的特征值.,42,记AA且有AQ 1R1.将等号右边两个矩阵因子的次序交换，得ARQ，且，(3)即AA.不难证明:即Ak+1AkA，矩阵序列Ak有相同的特征值.记,容易得到 是Ak的一个QR分解,43,相似约化为上Hessenberg矩阵,对一般n阶矩阵，QR算法的每一个迭代步需要O(n)次乘法运算.如果矩阵阶数稍大，这个算法几乎没有实际的应用价值.通常采用的方法是先将矩阵相似约化为上Hessenberg形式的矩阵，在此基础上应用QR迭代.这时，一个QR迭代步的乘法运算次数只需O(n）次.,44,所谓上Hessenberg矩阵是指一个n阶矩阵A，如果当ij+1时，aij=0，称A为上Hessenberg矩阵.例如：一个5阶的上Hessenberg矩阵具有如下的形式：,下面介绍QR方法时，都假设矩阵A是一个上Hessenberg矩阵.,45,Hessenberg矩阵的 方法,先把 经相似变换约化为Hessenberg矩阵，即：,且有很多零元.,设 为Hessenberg矩阵,作 分解:,46,问题在于 是否仍为Hessenberg矩阵？,可以证明：,若 为Hessenberg矩阵，,则：仍为Hessenberg矩阵。,47,如果A是一个满秩的上Hessenberg矩阵，可以证明，经过一个QR迭代步得到的A2QA1Q仍然是上Hessenberg矩阵.因为上Hessenberg矩阵次对角线以下的元素全为0，因此，只要证明，当k时，由迭 代格式产生的矩阵Ak的次对角元趋向于零就可以了.,48,QR算法的收敛性,定理设n阶矩阵A的n个特征值满足|n|0，其相应的n个线性无关特征向量为x1，x2，xn.记X(x1,x2,xn),Y=X.如果Y存在LU分解，那么，由迭代式产生的矩阵Ak基本收敛于上三角矩阵R.这里，基本收敛的含义指Ak的元素中除对角线以下的元素趋于零外，可以不收敛于R的元素.,49,QR算法的迭代过程,一个QR迭代步的计算 对l=1，2，n-1，构造n-1个平面旋转矩阵Pl,l+1,使A1的次对角元全部零化，实现A1的QR分解的计算，这里,50,用Pl,l+1右乘(24)，所得结果也放回矩阵A相应的元素中.,51,QR算法的迭代控制,当迭代步数k充分大时，由迭代格式产生的Ak的次对角元趋于0.在 实 际计算中，控制迭代次数常用的一种办法是，预先给定一个小的正数，在一个迭代步的计 算结束后，对l=n-1,n-2,，1，依次判别次对角元的绝对值是否满足 或更严格的准则是 或不太严格的准则是 如果上面三个不等式中有一个成立，把 看做实际上为零.,52,带原点位移的算法,由算法收敛性可以看出，算法的收敛速度 依赖于矩阵相邻特征值的比 值.为了加快算法的收敛速度，在迭代过程中，对矩阵Ak确定一个原点位移量sk，称Ak-skI为带原点位移量的矩阵，再对Ak-skI应用算法.这时，迭代格式改为 称为带原点位移的QR算法,53,计算特征值问题的QR方法，实际上总是分成2个阶段：,54,Jacobi方法是求实对称矩阵全部特征值和特征向量的一种矩阵变换方法。,4 Jacobi 方法,实对称矩阵A具有下列性质：,（1）A的特征值均为实数；,（2）存在正交矩阵R，使RTAR=diag(1,2,n),而,55,R的第i个列向量恰为i的特征向量;,直接求正交矩阵R是困难的.Jacobi提出用一系列所谓平面旋转矩阵逐次将A约化为对角矩阵.,平面解析几何中的平面坐标旋转变换,表示平面上坐标轴旋转角的变换.,（3）若记A1=RTAR,则A1仍为对称矩阵.,平面旋转矩阵,在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为,56,Rpq()具有下列性质:,一般地,在n维向量空间Rn中,沿着xpyq平面旋转角的变换矩阵为,称Rpq()为平面旋转矩阵.,57,设实对称矩阵A=(aij)nn,记B=RpqT()ARpq()=(bij)nn则它们元素之间有如下关系:,（1）Rpq()为正交矩阵,即Rpq-1()=RpqT();,（2）如果A为对称矩阵,则RpqT()ARpq()也为对称矩阵,且与A有相同的特征值.,（3）RpqT()A仅改变A的第p行与第q行元素,ARpq()仅改变A的第p列与第q列元素.,58,所以有,从而,有(5.5)、(5.6)式可得,如果apq0,适当选取角,使,59,只需角满足,从而,如果取|apq|=,若记,于是,则上式可记为,60,由式(5.7),令t=tan,则t满足方程,t2+2t-1=0,经典Jacobi算法是对A(0)=A施行一系列平面旋转变换:,为保证|/4,取绝对值较小的根,有,于是,Jacobi 方法,A(1)=R1TA(0)R1,A(2)=R2TA(1)R2,A(k)=RkTA(k-1)Rk,每一步变换选择A(k-1)=(aij(k-1)nn 的非对角线元素中绝对值最大者apq(k-1)(称为主元素)作为歼灭对象,构造平面旋,61,是给定的精度要求,则A的特征值可取为iaii(k),i=1,2,n.,转矩阵Rk=Rpq(),经变换得到A(k)=(aij(k)nn,且apq(k)=0,这时由(5.8)式有,从而,由此递推得到,当k充分大时,或者(A(k),或者,另外,由于 A(k)=RkTA(k-1)Rk=RkTRk-1TR1TAR1R2Rk=RTAR,62,的全部特征值.,解 记 A(0)=A,取p=1,q=2,apq(0)=a12(0)=2,于是有,因此,R=R1R2Rk 的列向量xj(j=1,2,n)为A的近似特征向量.,例7 用Jacobi 方法计算对称矩阵,从而有,63,所以,再取p=2,q=3,apq(1)=a23(1)=2.020190,类似地可得,以下依次有,64,65,从而A的特征值可取为 12.125825,28.388761,34.485401,为了减少搜索非对角线绝对值最大元素时间,对经典的Jacobi方法可作进一步改进.,1.循环Jacobi方法:按(1,2),(1,3),(1,n),(2,3),(2,4),(2,n),(n-1,n)的顺序,对每个(p,q)的非零元素apq作Jacobi变换,使其零化,逐次重复扫描下去,直至(A)为止.,2.过关Jacobi方法:取单调下降收敛于零的正数序列k,先以1为关卡值,依照1中顺序,将绝对值超过1的非对角元素零化,待所有非对角元素绝对值均不超过1时,再换下一个关卡值2,直到关卡值小于给定的精度.,66,Jacobi方法具有方法简单紧凑,精度高,收敛较快等优点,是计算对称矩阵全部特征值和相应特征向量的有效方法,但计算量较大,一般适用于阶数不高的矩阵.,67,矩阵的约化,68,Householder 矩阵,称 为Householder矩阵.,性质：1）（Hermite矩阵）,2）正交性：,69,正交矩阵 作用于 上，仍有：,即不改变向量的长度.,证明：若，则只需取 即可.,若，并确定w使,据此应有：,70,即：w应与向量 平行.,因为，所以,又因为，所以可取,这时 即为所求的Householder矩阵.,可以设计，使得 变为所需要的形状.,求（找），使得：,71,这里：,还可构造，使得：,72,即要求，且 的后面 个元素为零.,73,设：作,使得：,令,74,则,显然，这样构造的 仍然是Householder矩阵.,约化矩阵为Hessenberg形式,相似变换：,其中，当,75,矩阵的 分解,分解：，为正交矩阵，为上三角阵.,作法：用 表示 的列向量，令,取Householder矩阵 使得,其中，则：,76,然后,构造Householder矩阵,其中 为 阶Householder矩阵，使得：,77,则，具有如下形式：,构造出一串Householder 矩阵，使,记，显然为正交矩阵.,78,求实对称矩阵特征值的二分法,设 已知存在，使得：,不妨设：,79,Sturm 序列,定义：多项式序列,80,有：,性质1，仅有实根.,这是因为矩阵 的任意k阶主子矩阵都是对称矩阵，,81,性质3 若 是 的根，则,82,事实上,因为,而,故有：,83,从而 每个的零点均为单重实零点.,证明：（归纳法）,由于,又：,故：,84,同理：,故：,于是由,得出,类似可得,于是又由 便得,对于，均有：,85,性质5 在 的邻域，所有多项式 都为正，,多项式序列 称为Sturm序列：,1.相邻多项式无公共根,2.若,3.相邻两多项式的根彼此交错,86,例如：,注：若，则取它的符号与前一个,同号.,例如：矩阵,87,此时,表6.4,88,89,证明：由性质5知，时 的同号数,现让 从 起从左向右变动。,当 不经过所有 的零点 时，显然不变；,90,事实上，对于充分小的,在 内 均不变号，,根据性质3和5以及关于 的符号规定，,它们的符号变化为表6.5，6.6之一.,总之，通过 的零点时，不改变.,91,92,当 通过 的零点 时，根据同样的论证，,计算实对称三对角矩阵特征值的二分法,93,求特征值 的二分法如下：,（1）取 的中点,（2）计算 的同号数,（3）若 则,否则,（4）若 则令,94,特征向量的计算,用反幂法求 的特征向量：,（优势特征值）,此时：,即：,这样可求得 相应于 的特征向量.,