大规模矩阵问题的Krylov.ppt

《大规模矩阵问题的Krylov.ppt》由会员分享，可在线阅读，更多相关《大规模矩阵问题的Krylov.ppt（64页珍藏版）》请在三一文库上搜索。

1、大规模矩阵问题的Krylov 子空间方法综述,Krylov子空间方法综述,背景介绍 投影方法 Krylov子空间及其标准正交基 Krylov子空间方法求解线性方程组 Krylov子空间方法求解矩阵的特征值 研究热点和尚未解决的问题,背景,大规模线性方程组的求解 很多科学工程计算问题都转化为求解方程组Ax=b.如偏微分方程组的差分格式，有限元方法离散得到刚度矩阵. 大规模矩阵特征值和特征向量的计算 工程计算领域十分常见。如量子物理中的Kohn-Sham方程求解化为哈密顿矩阵某些关键特征值对的计算.,投影方法,线性方程组的投影方法 方程组Ax=b,A是nn的矩阵. 给定初始x(0),在m维空间K(

2、右子空间）中寻找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为对称正定矩阵,

3、x(0)为初始近似解,且K=L,则x(1)为采用投影方法得到的新近似解的充要条件是 其中, 且x为(1.1)的精确解.,. (残量投影)设A为任意方阵, x(0)为初始近似解,且 L=AK, 则x(1)为采用投影方法得到的新近似解的充要条件是 其中,矩阵特征值的投影方法 对于特征值问题Ax=x,其中A是nn的矩阵，斜交投影法是在m维右子空间K中寻找xi和复数i满足Axi- ixi L,其中L为m维左子空间.当L=K时，称此投影方法为正交投影法.,Kyrlov子空间,定义为m维Krylov子空间 为v的次数，即使得q(A)v的非零首一多项式的最低次数,Krylov子空间的标准正交基,Arnold

4、i方法 基于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+1,j=0, 则说明v的次数是j,且Kj是A的不变子空间.,定理 如果记以v(1),v(2),v(m)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为 ,删除最后一行得到的矩阵为Hm

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

6、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, 计算hij与w(j).,但采用IOM后，仍然需要存储v(1), v(2), v(m),因为在第

7、(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算法中，可能有较大舍入误差，改写,A是非对称阵时,H是Hessenberg阵，则当A是对称阵时，H也为对称阵，故应为三对角阵，

8、记为 对它采用修正Arnoldi正交化过程，就得到Lanczos方法,Krylov子空间方法解线性方程组,误差投影型方法 取L=K时的正交投影法 ）非对称矩阵的FOM方法 2 ) 非对称矩阵的IOM方法和DIOM方法 ）对称矩阵Lanczos算法 4 ）正定对称阵的CG法,CG法是解正定对称线性方程组最有效的方法之一 该方法充分利用了矩阵A的稀疏性，每次迭代的主要计算量是向量乘法而实际上，CG方法也是一种Krylov子空间方法 后面将给出一个数值例子,Krylov子空间方法解线性方程组,残量投影型方法 取L=AK时的斜交投影法 ）GMERS方法,GMERS方法把方程组的求解化为一个极小问题的求

9、解，回顾之前介绍投影类方法,回顾,. (残量投影)设A为任意方阵, x(0)为初始近似解,且 L=AK, 则x(1)为采用投影方法得到的新近似解的充要条件是 其中,GMERS方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法 不去直接求x(1),转而去解此极小问题，这是GMERS方法的独到之处.再由之前的定理,回顾,定理 如果记以v(1),v(2),v(m)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为 ,删除最后一行得到的矩阵为Hm,则 在Arnoldi算法中，可能有较大舍入误差，改写,GMERS方法在Arnoldi正交化过程之后把方程组的求解化

10、为一个极小问题的求解，回顾之前介绍投影类方法 不去直接求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)正交

11、就能得到QGMERS方法. 为了避免存储v(1),v(2)v(m-k)，可以对 进行QR分解.这种算法被称为DQGMERS,Krylov子空间方法解矩阵特征值问题,Arnoldi方法求解非对称矩阵特征值 Lanczos方法求解对称矩阵特征值,Arnoldi方法 再次回顾前面的定理 而 所以有以下结论：,1）若 ,则 是A的不变子空间,从而Hm的所有特征值是A的特征值子集 2）若 ,设Hm的特征值对是 ,即 ,由上述定理可得,以上讨论说明,可以用Hm的特征值 (又称Ritz值)来近似A的特征值,用向量 (又称Ritz向量)来近似相应的特征向量.事实上, Hm为A在Kyrlov子空间上的正交投影.

12、 由于Hm是上Hessenberg阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解.,Arnoldi算法构造标准正交基 1需要存储所有的基向量,当m很大时,存储量大 2理论上为了保证收敛速度,m越大越好 矛盾！,Lanczos 方法 A是对称阵时,Hm是三对角阵，仍然采用Hm的特征值来近似A的特征值，相应的Ritz向量来近似A的特征向量. 问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，QR法等，可以在并行机上达到tO(logN)的量级.,对称Lanczos算法,在没有舍入误差的情形下,得到的是标准正交基相互正交的,并且至多n步必然终

13、止.但是在出现舍入误差的情况下，计算得到的将失去正交性.因此,长期以来被人们认为是不稳定的.直到1971年,Paige通过仔细的舍入误差分析,发现失去正交性恰与近似特征值的精度提高有关.之后,经过大量的理论分析和数值实验,人们充分认识到Lanczos方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一. 目前，Lanczos方法是求大型稀疏对称矩阵部分极端特征值的最有效方法,改进技术,预条件法 Krylov子空间方法的收敛性一般都和矩阵的特征值分布有关，特征值分布越集中，收敛速度越快，若特征值分布较分散，则收敛的很慢，此时可以采用预条件子来改善矩阵的特征值分布. 选取矩阵 使得A的特征值比A集

14、中，并且M-1容易求得，于是将原问题化为问题M-1Ax=M-1b,这是左预条件法，还可采用右预条件法和分裂预条件法,Lanczos双正交化方法 对于非对称线性方程组，也可采用类似于Lanczos算法将非对称阵化为三对角阵，它的好处在于可以可以形成短迭代，而FOM、IOM、GMERS的计算量和存储量都随着m的增大而增大 在双正交化过程中，取 容易看出A和在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是A和AT的方程组.,多项式加速技术 在求解矩阵特征值问题时，可以利用chebeshev多项式来加快收敛 所谓多项式加速,就是采用 作为迭 代初始向量,其中 为多项式,如果 , 其中p

15、i为A的特征值对应的特征向量 则 将特征值按实部递减顺序排列,其中为r个“想要”的特征值,将“不想要”的特征值点集用S表示.如果多项式 ,能使得S尽可能小,则相当于让求解过程产生加速.而Chebyshev多项式由于它的特殊性质,恰好能满足这种要求.不过，至今,仍无有效方法确定Chebyshev多项式的某些参数.,块方法、调和投影法 当矩阵的特征值分布较密集或有重特征值时，常规的Arnoldi方法或Lanczos方法就必须取m很大时才能收敛，此时可以采用块方法 取 可以使得无须很大的m就可让迭代收敛，并且此方法适合改为并行算法。 当特征值密集时，还可采用调和投影法，这也是最近研究的一个热点。当左

16、右子空间的基满足Wm=AVm时，称为调和投影，它除了可以通过选取位移点改善原矩阵的特征值分布外，还可以将内部特征值问题化为外部特征值问题.,尚待解决的问题,在Arnodli过程和Lanczos双正交化过程中，会出现崩溃，怎样的矩阵容易出现崩溃。目前已有一些处理此崩溃的办法如前瞻技术，但并不十分成熟. 在投影类方法中Krylov方法是否是最优的,K是否有更好的选择. 怎样将预条件技术运用到求解矩阵特征值问题,CG法解方程组Ax=b,输入：cg(A,b,tol) /tol为误差限. 输出：初始向量x0(随机产生),cputime, 方程的解x,并绘制残差随迭代步变化的曲线.,function x=

17、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),