大规模矩阵问题的Krylov.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,大规模矩阵问题的,Krylov,子空间方法综述,Krylov,子空间方法综述,背景介绍,投影方法,Krylov,子空间及其标准正交基,Krylov,子空间方法求解线性方程组,Krylov,子空间方法求解矩阵的特征值,研究热点和尚未解决的问题,背景,大规模线性方程组的求解,很多科学工程计算问题都转化为求解方程组,Ax=b.,如偏微分方程组的差分格式，有限元方法离散得到刚度矩阵,.,大规模矩阵特征值和特征向量的计算,工程计算领域十分常见。如量子物理中的,Kohn-Sham,方程求解化为哈密顿矩阵某些关键特征值对的计算,.,投影方法,线性方程组的投影方法,方程组,Ax=,b,A,是,nn,的矩阵,.,给定初始,x,(0),在,m,维空间,K(,右子空间）中寻找,x,的近似解,x,(1),满足残向量,r=b-Ax,(1),与,m,维空间,L(,左子空间）正交，即,b-Ax,(1),L,此条件称为,Petrov-Galerkin,条件,.,当空间,K=L,时，称相应的投影法为正交投影法，否则称为斜交投影法,.,K(L),X,(0),X,(1),K,L,X,(0),X,(1),解方程组的投影法的矩阵表示,设,nm,阶矩阵,V=v,(1),v,(2),v,(m,),与,W=w,(1),w,(2),w,(m,),的列分别构成,K,与,L,的一组基,.,记,z=x,(1),-x,(0),z=,Vy,有,当非奇异时（通常情况成立，见定理,1.1,），有,从而得到迭代公式,投影方法的最优性,.(,误差投影,),设,A,为对称正定矩阵,x,(0),为初始近似解,且,K=L,则,x,(1),为采用投影方法得到的新近似解的充要条件是,其中,且,x,为,(1.1),的精确解,.,.(,残量投影,),设,A,为任意方阵,x,(0),为初始近似解,且,L=AK,则,x,(1),为采用投影方法得到的新近似解的充要条件是,其中,矩阵特征值的投影方法,对于特征值问题,Ax=,x,其中,A,是,nn,的矩阵，斜交投影法是在,m,维右子空间,K,中寻找,x,i,和复数,i,满足,Ax,i,-,i,x,i,L,其中,L,为,m,维左子空间,.,当,L=K,时，称此投影方法为正交投影法,.,Kyrlov,子空间,定义为,m,维,Krylov,子空间,为,v,的次数，即使得,q(A)v,的非零首一多项式的最低次数,Krylov,子空间的标准正交基,Arnoldi,方法,基于,Gram-,Schmit,正交化方法,首先，选取一个,Euclid,范数为,1,的向量,v,(1),对,通常可取 在已知,v,(1),v,(2),v,(j,),的情况下,不妨设,v,(1),v,(2),v,(j),Av,(j,),线性无关,(,否则构造完毕,),则可求出与每个都正交的向量,而不难看出，再记,得到与,v,(1),v,(2),v,(j,),都,正交的向量 重复此过程，即可,得到一组标准正交基,.,若期间某个,j,使得,h,j+1,j,=0,则说明,v,的次数是,j,且,K,j,是,A,的不变子空间,.,定理,如果记以,v,(1),v,(2),v,(m,),为列构成的矩阵为,V,m,由,h,ij,定义的,(m+1)m,阶上,Hessenberg,矩阵为,删除最后一行得到的矩阵为,H,m,则,在,Arnoldi,算法中，可能有较大舍入误差，改写,Krylov,子空间方法解线性方程组,误差投影型方法,取,L=K,时的正交投影法,）非对称矩阵的,FOM,方法,解方程组的投影法的矩阵表示,设,nm,阶矩阵,V=v,(1),v,(2),v,(m,),与,W=w,(1),w,(2),w,(m,),的列分别构成,K,与,L,的一组基,.,记,z=x,(1),-x,(0),z=,Vy,有,当非奇异时（通常情况成立，见定理,1.1,），有,从而得到迭代公式,回顾,不难看出，当采用上述,FOM,算法时，需要存储所有的,v,i,，,(i=1,2,m),，当,m,增大时，存储量以,O(mn,),量级增大,.,而,FOM,计算量是,O(m,2,n).,可见其代价十分高昂,.,因此我们考虑重启的,FOM,算法,Krylov,子空间方法解线性方程组,误差投影型方法,取,L=K,时的正交投影法,）非对称矩阵的,FOM,方法,2),非对称矩阵的,IOM,方法和,DIOM,方法,所谓不完全正交化方法,(IOM),是指在正交化过程中,v,(j+1),仅与最近,k,个,v,(j-k+1),v,(j-1),v,(j),正交,这样做虽然破坏的正交性,但是降低了计算量,.,当然,k,选得越小，对每个,j,对应的计算量也越小,但可能要选更大的,m,才能取得满足精度要求的近似解,.,IOM,算法仅仅是把,FOM,算法的第三步改为,(iii),对,i=max(1,j-k+1),j,计算,h,ij,与,w,(j,),.,但采用,IOM,后，仍然需要存储,v,(1),v,(2),v,(m,),因为在第,(vi),步中仍然需要这些向量,.,解决这个问题可以考虑采用,H,的,LU,分解，通过自身分解的迭代更新以减少每一步的存储量,Krylov,子空间方法解线性方程组,误差投影型方法,取,L=K,时的正交投影法,）非对称矩阵的,FOM,方法,2),非对称矩阵的,IOM,方法和,DIOM,方法,）对称矩阵,Lanczos,算法,回顾,定理,如果记以,v,(1),v,(2),v,(m,),为列构成的矩阵为,V,m,由,h,ij,定义的,(m+1)m,阶上,Hessenberg,矩阵为,删除最后一行得到的矩阵为,H,m,则,在,Arnoldi,算法中，可能有较大舍入误差，改写,A,是非对称阵时,H,是,Hessenberg,阵，则当,A,是对称阵时，,H,也为对称阵，故应为三对角阵，记为,对它采用修正,Arnoldi,正交化过程，就得到,Lanczos,方法,Krylov,子空间方法解线性方程组,误差投影型方法,取,L=K,时的正交投影法,）非对称矩阵的,FOM,方法,2),非对称矩阵的,IOM,方法和,DIOM,方法,）对称矩阵,Lanczos,算法,4,）正定对称阵的,CG,法,CG,法是解正定对称线性方程组最有效的方法之一,该方法充分利用了矩阵,A,的稀疏性，每次迭代的主要计算量是向量乘法而实际上，,CG,方法也是一种,Krylov,子空间方法,后面将给出一个数值例子,Krylov,子空间方法解线性方程组,残量投影型方法,取,L=AK,时的斜交投影法,）,GMERS,方法,GMERS,方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法,回顾,.(,残量投影,),设,A,为任意方阵,x,(0),为初始近似解,且,L=AK,则,x,(1),为采用投影方法得到的新近似解的充要条件是,其中,GMERS,方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法,不去直接求,x,(1),转而去解此极小问题，这是,GMERS,方法的独到之处,.,再由之前的定理,回顾,定理,如果记以,v,(1),v,(2),v,(m,),为列构成的矩阵为,V,m,由,h,ij,定义的,(m+1)m,阶上,Hessenberg,矩阵为,删除最后一行得到的矩阵为,H,m,则,在,Arnoldi,算法中，可能有较大舍入误差，改写,GMERS,方法在,Arnoldi,正交化过程之后把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法,不去直接求,x,(1),转而去解此极小问题，这是,GMERS,方法的独到之处,.,再由之前的定理,Krylov,子空间方法解线性方程组,残量投影型方法,取,L=AK,时的斜交投影法,）,GMERS,方法,2),重启型,GMERS,方法、,QGMERS,、,DGMERS,GMERS,算法具有很好的收敛性，在不考虑舍入误差的情况下，该算法能在在有限步得到精确解,.,但是，和,FOM,算法类似，它需要保存大量正交基，因此我们可以采取类似重启型,FOM,的算法实现重启型,GMERS,算法,同样可以采取不完全正交的方法，在正交化过程中,v,(j+1),仅与最近,k,个,v,(j-k+1),v,(j-1),v,(j),正交就能得到,QGMERS,方法,.,为了避免存储,v,(1),v,(2),v,(m-k,),，可以对 进行,QR,分解,.,这种算法被称为,DQGMERS,Krylov,子空间方法解矩阵特征值问题,Arnoldi,方法求解非对称矩阵特征值,Lanczos,方法求解对称矩阵特征值,Arnoldi,方法,再次回顾前面的定理,而,所以有以下结论：,1,）若,则 是,A,的不变子空间,从而,H,m,的所有特征值是,A,的特征值子集,2,）若,设,H,m,的特征值对是,即,由上述定理可得,以上讨论说明,可以用,H,m,的特征值,(,又称,Ritz,值,),来近似,A,的特征值,用向量,(,又称,Ritz,向量,),来近似相应的特征向量,.,事实上,H,m,为,A,在,Kyrlov,子空间上的正交投影,.,由于,H,m,是上,Hessenberg,阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解,.,Arnoldi,算法构造标准正交基,1,需要存储所有的基向量,当,m,很大时,存储量大,2,理论上为了保证收敛速度,m,越大越好,矛盾！,Lanczos,方法,A,是对称阵时,H,m,是三对角阵，仍然采用,H,m,的特征值来近似,A,的特征值，相应的,Ritz,向量来近似,A,的特征向量,.,问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，,QR,法等，可以在并行机上达到,tO(logN,),的量级,.,对称,Lanczos,算法,在没有舍入误差的情形下,得到的是标准正交基相互正交的,并且至多,n,步必然终止,.,但是在出现舍入误差的情况下，计算得到的将失去正交性,.,因此,长期以来被人们认为是不稳定的,.,直到,1971,年,Paige,通过仔细的舍入误差分析,发现失去正交性恰与近似特征值的精度提高有关,.,之后,经过大量的理论分析和数值实验,人们充分认识到,Lanczos,方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一,.,目前，,Lanczos,方法是求大型稀疏对称矩阵部分极端特征值的最有效方法,改进技术,预条件法,Krylov,子空间方法的收敛性一般都和矩阵的特征值分布有关，特征值分布越集中，收敛速度越快，若特征值分布较分散，则收敛的很慢，此时可以采用预条件子来改善矩阵的特征值分布,.,选取矩阵 使得,A,的特征值比,A,集中，并且,M,-1,容易求得，于是将原问题化为问题,M,-1,Ax=M,-1,b,这是左预条件法，还可采用右预条件法和分裂预条件法,Lanczos,双正交化方法,对于非对称线性方程组，也可采用类似于,Lanczos,算法将非对称阵化为三对角阵，它的好处在于可以可以形成短迭代，而,FOM,、,IOM,、,GMERS,的计算量和存储量都随着,m,的增大而增大,在双正交化过程中，取,容易看出,A,和在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是,A,和,A,T,的方程组,.,多项式加速技术,在求解矩阵特征值问题时，可以利用,chebeshev,多项式来加快收敛,所谓多项式加速,就是采用 作为迭,代初始向量,其中 为多项式,如果,其中,p,i,为,A,的特征值对应的特征向量,则,将特征值按实部递减顺序排列,其中为,r,个“想要”的特征值,将“不想要”的特征值点集用,S,表示,.,如果多项式,能使得,S,尽可能小,则相当于让求解过程产生加速,.,而,Chebyshev,多项式由于它的特殊性质,恰好能满足这种要求,.,不过，至今,仍无有效方法确定,Chebyshev,多项式的某些参数,.,块方法、调和投影法,当矩阵的特征值分布较密集或有重特征值时，常规的,Arnoldi,方法或,Lanczos,方法就必须取,m,很大时才能收敛，此时可以采用块方法,取,可以使得无须很大的,m,就可让迭代收敛，并且此方法适合改为并行算法。,当特征值密集时，还可采用调和投影法，这也是最近研究的一个热点。当左右子空间的基满足,W,m,=,AV,m,时，称为调和投影，它除了可以通过选取位移点改善原矩阵的特征值分布外，还可以将内部特征值问题化为外部特征值问题,.,尚待解决的问题,在,Arnodli,过程和,Lanczos,双正交化过程中，会出现崩溃，怎样的矩阵容易出现崩溃。目前已有一些处理此崩溃的办法如前瞻技术，但并不十分成熟,.,在投影类方法中,Krylov,方法是否是最优的,K,是否有更好的选择,.,怎样将预条件技术运用到求解矩阵特征值问题,CG,法解方程组,Ax=b,输入：,cg(A,b,tol,)/,tol,为误差限,.,输出：初始向量,x0(,随机产生,),cputime,方程的解,x,并绘制残差随迭代步变化的曲线,.,function x=,cg(A,b,tol,),t=,cputime,;,n=,numel(b,);,x0=rand(n,1),r=b-A*x0;,s=b-A*x0;,x=x0;,j=1;,while(,norm(r,),tol,),count(j,)=j-1;,res(j,)=,norm(r,);,alpha(j,)=,dot(r,r)/dot(s,(A,*r);,b1=,dot(r,r,);,(,续,),x=,x+(alpha(j,)*s;,r=,r-(alpha(j,)*(A*s);,b2=,dot(r,r,);,beta(j+1)=b2/b1;,s=r+(beta(j+1)*s;,j=j+1;,end,t=,cputime-t,plot(count,res),xlabel,(,迭代步,),ylabel,(,残量,),title(,残量随迭代步的变化,),set(findobj(gca,Type,line,Color,0 0 1),.,Color,red,.,LineWidth,2),