矩阵特征值问题的解法.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,类似地,当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-11 step 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矩阵.,约化矩阵为h