大规模矩阵问题的Krylov.ppt

《大规模矩阵问题的Krylov.ppt》由会员分享，可在线阅读，更多相关《大规模矩阵问题的Krylov.ppt（64页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、大规模矩阵问题的Krylov Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望Krylov子空间方法综述背景介绍投影方法Krylov子空间及其标准正交基Krylov子空间方法求解线性方程组Krylov子空间方法求解矩阵的特征值研究热点和尚未解决的问题背景大规模线性方程组的求解很多科学工程计算问题都转化为求解方程组Ax=b.如偏微分方程组的差分格式，有限元方法离散得到刚度矩阵.大规模矩阵特征值和特征向量的计算工程计算领域十分常见。如量子物理中的Kohn-Sham方程求

2、解化为哈密顿矩阵某些关键特征值对的计算.投影方法线性方程组的投影方法方程组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)KLX(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

3、.1），有从而得到迭代公式投影方法的最优性.(误差投影)设A为对称正定矩阵,x(0)为初始近似解,且K=L,则x(1)为采用投影方法得到的新近似解的充要条件是其中,且x为(1.1)的精确解.(残量投影)设A为任意方阵,x(0)为初始近似解,且 L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是其中矩阵特征值的投影方法对于特征值问题对于特征值问题Ax=x,Ax=x,其中其中A A是是nnnn的矩阵，斜的矩阵，斜交投影法是在交投影法是在mm维右子空间维右子空间K K中寻找中寻找x xi i和复数和复数 i i满满足足AxAxi i-i ix xi i L,L,其中其中L L为为mm维左子

4、空间维左子空间.当当L=KL=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使得hj+

5、1,j=0,则说明v的次数是j,且Kj是A的不变子空间.定理定理 如果记以v(1),v(2),v(m)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为,删除最后一行得到的矩阵为Hm,则在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），有从而得到迭代公式回顾不难看出，

6、当采用上述FOM算法时，需要存储所有的vi，(i=1,2,m)，当m增大时，存储量以O(mn)量级增大.而FOM计算量是O(m2n).可见其代价十分高昂.因此我们考虑重启的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算法

7、的第三步改为(iii)对i=max(1,j-k+1),j,计算hij与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)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为,删除最后一行得到的矩阵为Hm,则在Arnoldi算法

8、中，可能有较大舍入误差，改写AA是非对称阵时是非对称阵时,H,H是是HessenbergHessenberg阵，则当阵，则当A A是对称阵是对称阵时，时，H H也为对称阵，故应为三对角阵，记为也为对称阵，故应为三对角阵，记为对它采用修正对它采用修正ArnoldiArnoldi正交化过程，就得到正交化过程，就得到LanczosLanczos方法方法Krylov子空间方法解线性方程组误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法2)非对称矩阵的IOM方法和DIOM方法）对称矩阵Lanczos算法4）正定对称阵的CG法CG法是解正定对称线性方程组最有效的方法之一该方法充分利用了矩阵A的

9、稀疏性，每次迭代的主要计算量是向量乘法而实际上，CG方法也是一种Krylov子空间方法后面将给出一个数值例子Krylov子空间方法解线性方程组残量投影型方法取L=AK时的斜交投影法）GMERS方法GMERS方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法回顾.(残量投影)设A为任意方阵,x(0)为初始近似解,且 L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是其中GMERS方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法不去直接求x(1),转而去解此极小问题，这是GMERS方法的独到之处.再由之前的定理回顾定理定理 如果记以v(1),v(2),v(m

10、)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为,删除最后一行得到的矩阵为Hm,则在Arnoldi算法中，可能有较大舍入误差，改写GMERS方法在Arnoldi正交化过程之后把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法不去直接求x(1),转而去解此极小问题，这是GMERS方法的独到之处.再由之前的定理Krylov子空间方法解线性方程组残量投影型方法取L=AK时的斜交投影法）GMERS方法2)重启型GMERS方法、QGMERS、DGMERSGMERS算法具有很好的收敛性，在不考虑舍入误差的情况下，该算法能在在有限步得到精确解.但是，和FOM算法类似

11、，它需要保存大量正交基，因此我们可以采取类似重启型FOM的算法实现重启型GMERS算法同样可以采取不完全正交的方法，在正交化过程中,v(j+1)仅与最近k个v(j-k+1),v(j-1),v(j)正交就能得到QGMERS方法.为了避免存储v(1),v(2)v(m-k)，可以对进行QR分解.这种算法被称为DQGMERSKrylov子空间方法解矩阵特征值问题Arnoldi方法求解非对称矩阵特征值Lanczos方法求解对称矩阵特征值Arnoldi方法再次回顾前面的定理而所以有以下结论：1）若,则是A的不变子空间,从而Hm的所有特征值是A的特征值子集2）若,设Hm的特征值对是,即,由上述定理可得以上讨

12、论说明,可以用Hm的特征值(又称Ritz值)来近似A的特征值,用向量(又称Ritz向量)来近似相应的特征向量.事实上,Hm为A在Kyrlov子空间上的正交投影.由于Hm是上Hessenberg阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解.Arnoldi算法构造标准正交基1需要存储所有的基向量,当m很大时,存储量大2理论上为了保证收敛速度,m越大越好矛盾！Lanczos方法A是对称阵时,Hm是三对角阵，仍然采用Hm的特征值来近似A的特征值，相应的Ritz向量来近似A的特征向量.问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，QR法等

13、，可以在并行机上达到tO(logN)的量级.对称对称LanczosLanczos算法算法,在没有舍入误差的情形下在没有舍入误差的情形下,得到得到的是标准正交基相互正交的的是标准正交基相互正交的,并且至多并且至多n n步必然终止步必然终止.但但是在出现舍入误差的情况下，计算得到的将失去正交是在出现舍入误差的情况下，计算得到的将失去正交性性.因此因此,长期以来被人们认为是不稳定的长期以来被人们认为是不稳定的.直到直到19711971年年,Paige,Paige通过仔细的舍入误差分析通过仔细的舍入误差分析,发现失去正交性恰与发现失去正交性恰与近似特征值的精度提高有关近似特征值的精度提高有关.之后之后

14、,经过大量的理论分经过大量的理论分析和数值实验析和数值实验,人们充分认识到人们充分认识到LanczosLanczos方法是求解大方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一型稀疏对称矩阵特征值问题的最有效方法之一.目前，目前，LanczosLanczos方法是求大型稀疏对称矩阵部分极端特征值方法是求大型稀疏对称矩阵部分极端特征值的最有效方法的最有效方法改进技术预条件法预条件法KrylovKrylov子空间方法的收敛性一般都和矩阵的特征值子空间方法的收敛性一般都和矩阵的特征值分布有关，特征值分布越集中，收敛速度越快，若分布有关，特征值分布越集中，收敛速度越快，若特征值分布较分散，则收敛的

15、很慢，此时可以采用特征值分布较分散，则收敛的很慢，此时可以采用预条件子来改善矩阵的特征值分布预条件子来改善矩阵的特征值分布.选取矩阵选取矩阵使得使得A A的特征值比的特征值比A A集中，并集中，并且且MM-1-1容易求得，于是将原问题化为问题容易求得，于是将原问题化为问题MM-1-1Ax=MAx=M-1-1b,b,这是左预条件法，还可采用右预条件法和分裂预条这是左预条件法，还可采用右预条件法和分裂预条件法件法LanczosLanczos双正交化方法双正交化方法对于非对称线性方程组，也可采用类似于对于非对称线性方程组，也可采用类似于LanczosLanczos算算法将非对称阵化为三对角阵，它的好

16、处在于可以可以形成法将非对称阵化为三对角阵，它的好处在于可以可以形成短迭代，而短迭代，而FOMFOM、IOMIOM、GMERSGMERS的计算量和存储量都随的计算量和存储量都随着着mm的增大而增大的增大而增大 在双正交化过程中，取在双正交化过程中，取容易看出容易看出A A和在其中扮演对偶的角色，此方法特别适用和在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是于需要求解两个系数矩阵分别是A A和和A AT T的方程组的方程组.多项式加速技术在求解矩阵特征值问题时，可以利用在求解矩阵特征值问题时，可以利用chebeshevchebeshev多项式来加快收敛多项式来加快收敛所谓多项式

17、加速所谓多项式加速,就是采用就是采用作为迭作为迭 代初始向量代初始向量,其中其中为多项式为多项式,如果如果,其其中中p pi i为为A A的特征值对应的特征向量的特征值对应的特征向量则则 将特征值按实部递减顺序排列将特征值按实部递减顺序排列,其中为其中为r r个个“想想要要”的特征值的特征值,将将“不想要不想要”的特征值点集用的特征值点集用S S表表示示.如果多项式如果多项式,能使得能使得S S尽可能小尽可能小,则相当于让则相当于让求解过程产生加速求解过程产生加速.而而ChebyshevChebyshev多项式由于它的多项式由于它的特殊性质特殊性质,恰好能满足这种要求恰好能满足这种要求.不过，

18、至今不过，至今,仍无仍无有效方法确定有效方法确定ChebyshevChebyshev多项式的某些参数多项式的某些参数.块方法、调和投影法块方法、调和投影法当矩阵的特征值分布较密集或有重特征值时，当矩阵的特征值分布较密集或有重特征值时，常规的常规的ArnoldiArnoldi方法或方法或LanczosLanczos方法就必须取方法就必须取mm很大时很大时才能收敛，此时可以采用块方法才能收敛，此时可以采用块方法取取可以使得无须很大的可以使得无须很大的mm就可让迭代收敛，并且此方法就可让迭代收敛，并且此方法适合改为并行算法。适合改为并行算法。当特征值密集时，还可采用调和投影法，这也是当特征值密集时，

19、还可采用调和投影法，这也是最近研究的一个热点。当左右子空间的基满足最近研究的一个热点。当左右子空间的基满足WWmm=AV=AVmm时，称为调和投影，它除了可以通过选取位时，称为调和投影，它除了可以通过选取位移点改善原矩阵的特征值分布外，还可以将内部特移点改善原矩阵的特征值分布外，还可以将内部特征值问题化为外部特征值问题征值问题化为外部特征值问题.尚待解决的问题在Arnodli过程和Lanczos双正交化过程中，会出现崩溃，怎样的矩阵容易出现崩溃。目前已有一些处理此崩溃的办法如前瞻技术，但并不十分成熟.在投影类方法中Krylov方法是否是最优的,K是否有更好的选择.怎样将预条件技术运用到求解矩阵

20、特征值问题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-tplot(count,res),xlabel(迭代步迭代步),ylabel(残量残量),title(残量随迭代步的变化残量随迭代步的变化),set(findobj(gca,Type,line,Color,0 0 1),.Color,red,.LineWidth,2)