大型矩阵迭代分析研究 计算机科学与技术专业

目录摘要31前言52迭代法的思想62.1迭代法的基本概念62.2迭代法的收敛条件63常见迭代法73.1Jacobi迭代法73.1.1Jacobi迭代法的计算公式73.1.2Jacobi迭代法的收敛条件83.2Gauss迭代法93.2.1Gauss迭代法的计算公式93.2.2Gauss迭代法的收敛条件103.3SOR迭代法103.3.1SOR迭代法的计算公式103.3.2SOR迭代法的收敛条件113.4CG迭代法113.5PCG迭代法123.6GMRES迭代法134数值算例分析13算例4.113算例4.2145总结16参考文献17致谢18摘要这篇论文介绍了大型矩阵运算时要用到的迭代法和其思想，并简单介绍了常用的几种迭代方法：如基于矩阵分解原理的Jacobi迭代法，Gauss迭代法和SOR迭代法；属于共轭方向法的CG法和PCG法；最后介绍了适用于解决非线性方程组的GMRES迭代法及他们其中一些迭代法的收敛条件.并通过对不同矩阵的具体运算比较了他们的迭代次数，简单分析了几种迭代方法的优劣.关键字:迭代法思想，Jacobi迭代法，Gauss迭代法，SOR迭代法,CG迭代法,PCG迭代法，GMRES迭代法，迭代次数，收敛体条件.AbstractThisessayintroducestheiterativemethodwhichisneededduringlargematrixoperation.Meanwhile,someofcommonmethodsofiterationarealsosuggestedinthisessay,forinstance:Jacobianinteration,GaussinterationandsuccessiveoverrelaxationmethodwhichwereallbasedontheMatrixdecompositiontechnique;conjugategradientalgorithmandPreprocessingconjugategradientalgorithmwhichbothbelongtoconjugatedirectionmethod;GMRESalgorithmwhichisusedtosolvenon-linearequationandsomeoftheirownconditionofconvergence.Inaddition,bycomparingtheiriterationsteps,Iroughlyanalysedtheirmeritsanddemerits.Keywords:iterativemethod，Jacobianinteration,Gaussinteration,successiveoverrelaxationmethod,conjugategradientalgorithm,Preprocessingconjugategradientalgorithm,GMRESalgorithm,iterationstepsandconditionofconvergence.1前言矩阵自19世纪被英国数学家凯利提出来，历经了200多年的发展，在现代生活中有及其广泛的应用，下面我将列举几项矩阵在现实生活中的应用：矩阵可用于预测水质量状况以及用于处理水污染等问题;矩阵还可以用于模拟飞机飞行所需要的飞行环境;在物理学领域的磁电场的探查测量方面也有很广泛的应用;许多图片在电脑上的存储方式就是以矩阵的格式存储的,以便对图片进行调整和美化的时候可以利用数学公式对图片进行处理等。除此以外，矩阵还在许多领域有着极为广泛的应用，这里就不一一举例了。而这些矩阵通常都是大型矩阵，如何高效的处理大型矩阵和与其相关的线性方程组成了一个值得深思的问题。而迭代法在处理零元素较多的大型稀疏矩阵时具有运算和存储两方面的优势。比较典型的有基于矩阵分解原理的Jacobi迭代法，这是由伟大的普鲁士数学家雅可比在19世纪提出来的，在当时Jacobi迭代法具有许多优点，其计算公式简单，每一次的迭代中只需计算一次矩阵和向量的乘法。这也是最早期的迭代法。同雅可比迭代法一样，高斯-塞德尔和超松弛迭代法也是基于矩阵分解原理，这些都是传统的迭代法，出现时间较久。此外，还有其他一些优异的矩阵迭代算法，如以极小化的方法来讨论方程组的解的著名的共轭梯度（CG）算法，CG法出现在20世纪中期，它有许多优异的性质：比如相较于传统迭代法，CG法收敛速度较快；此外CG法只有极少量的非向量运算。然而由于实际运算中舍入误差的积累以及正交性的逐渐丧失，CG法长时间并未得到广泛应用。直到20世纪70年代，预处理共轭梯度算法（PCG算法）出现了，它通过对CG算法中未加工过的矩阵进行某种形式的加工改造，使它具有更适宜迭代的特点。近些来，PCG法有了进一步的发展。它现在已经有了MICCG，块预处理，高阶LU分解和行和一致预处理等预处理方法，预处理技术正在不断完善中；此外共轭梯度算法也发展为可应用于不定阵，非奇异阵及复矩阵等多种矩阵形式的较为系统的方法。经过了近200年的发展，迭代法变得越来越系统，且新的优异的迭代法也不断涌出，本文主要对其中几种比较典型的迭代方法进行介绍和探究，本文的内容主要分为：迭代法基本概念及收敛性的介绍，几种迭代法步骤及特性的介绍(Jacobi,Gauss,SOR,CG,PCG,GMRES)，具体的算例分析和总结。2迭代法的思想2.1迭代法的基本概念迭代法的由来与线性方程组的求解有较深的关系，下面我们来举例说明有线性方程组Ax=b,（2.1）其中A为大型稀疏矩阵（具有较多的零元素），此时直接求解很困难，而用迭代法是比较合理有效的.例如：求解线性方程组8x1−3x2记为Ax=b，其中有A=8−32411−16312,x该方程组的精确解为x∗=3，2，1x1=183x则有x=B0x+fB0=038−2任意选取初始值，例如x(0)=0，0，0T，代入（2.3）式可以得到新的值x(1)=x1(1)，xx然而向量序列x(k)2.2迭代法的收敛条件不难看出若有limk→∞x(由上可知，若要研究x(ε(对于线性方程组x=Bx+f,设有唯一解为x∗x∗=Bx∗+与x(k+1)=B0x(k)+f相减得，ε(k+1)=ε(k)=Bε要考察x(k)的收敛性，就要研究B在什么条件下limk→∞ε3常见的几种迭代法3.1Jacobi迭代法3.1.1Jacobi迭代法的计算公式上面迭代法的思想是从线性方程组的角度来看，我们不妨从矩阵的角度来看.矩阵分解技术在矩阵中具有广泛应用，它通过对矩阵的分解来构造迭代格式.如将矩阵A分解为M+N故此Ax=b即可写成Mx+Nx=b即Mx=b-Nx当M可逆时有x=由此可构造迭代格式x其中−M−1N为迭代矩阵，且M下面我们将讨论几种基于矩阵分解原理的常见迭代法.将线性方程组中的系数A=(aA=a11_0:=D-L-U.设aii≠0(i=1,2…,n),选取M为A的对角元素部分，即选取M=D（对角矩阵），A=Dx其中B=I-D−1A=D−1（L+U）≡BJ，f=D−1b.由迭代法有公式Dx(k+1)=（L+U）x(k)或aiixi(k+1)因此，对应于求解线性方程组Ax=b的雅可比迭代法的计算公式为x雅可比迭代法计算形式比较简单，每一次只需将算出来的值再代入迭代算式里，循环反复即可得到结果，但也是因为这个原因，迭代次数往往较多，故在实际生活中已经很少被应用.3.1.2雅可比迭代法的收敛性下面我们来讨论雅可比迭代法的收敛性.由上面介绍的迭代法的收敛条件可知迭代法的收敛性于它自身的谱半径有关，有定理如下定理3.1：设Ax=b，其中A=D-L-U为非奇异矩阵，且对角矩阵D也非奇异，则解线性方程组雅可比迭代法收敛的充要条件是ρBJ<1,其中BJ=D然而实际运用中，大型矩阵的谱半径求解运算量十分大，因而我们需要其他的一些判定条件.再给出判定条件之前，先简单介绍一下对角占优矩阵和可约与不可约矩阵的概念.（对角占优矩阵）设A=若A的元素满足a称A为严格对角占优矩阵.若A的元素满足a且上式至少有一个不等式成立，则称A为弱对角占优矩阵[1](可约与不可约矩阵)设A=(aij)n×n(nPT其中A11为r阶方程，A22为n-r阶方程（1≤r<n）,则称A为可约矩阵，否则，若不存在这样的置换矩阵P使上式成立，则称A介绍了这几个矩阵的概念之后，下面给出关于雅可比矩阵收敛的几个定理：定理3.2：设Ax=b，如果A为严格对角占优矩阵，则解Ax=b的雅可比迭代法收敛；如果A为弱对角占优矩阵，且A为不可约矩阵，则解Ax=b的雅可比迭代法收敛[1]定理3.3：设矩阵A对称，且对角元aii>0(i=1,2,…,n),则解线性方程组Ax=b的雅可比迭代法收敛的充分必要条件是A及2D-A均为正定矩阵，其中D=diag(3.2Gauss迭代法3.2.1Gauss迭代法的计算公式类似的，Gauss迭代法也是根据矩阵分解技术，将A的下三角部分作为分裂矩阵M，即M=D-L（下三角矩阵），A=M-N，这样就可以得到Ax=b的Gauss迭代法x其中B=I-（D−L)−1A=(D−L)−1U≡BGs，f=由迭代法可构造Gauss迭代公式(D−L)x(k+1)=U或写为Dx(k+1)=Lx(k+1)+即a于是，高斯-塞德尔迭代法的计算公式为x同Jacobi迭代法一样，Gauss迭代法也是比较传统的迭代方法，出现了近200年，早已不适合实际运用.3.2.2高斯-塞德尔迭代的收敛性与雅可比迭代法一样，解线性方程组的高斯-塞德尔迭代法收敛的充要条件是谱半径ρBGs<同样的由于实际中谱半径的计算问题，关于高斯-塞德尔迭代法收敛性也有其他一些定理定理3.4：对于Ax=b，如果A为严格对角占优矩阵，则解Ax=b的高斯-塞德尔迭代法收敛；如果A为弱对角占优矩阵，且A为不可约矩阵，则解Ax=b的高斯-塞德尔迭代法收敛[5]定理3.5：设矩阵A对称，且对角元aii>0(i=1,2,…,n),则解线性方程组Ax=b的高斯-塞德尔迭代法收敛的充分必要条件是A正定[3.3SOR迭代法3.3.1SOR迭代法的计算公式出现于20世纪70年代的超松弛迭代法（简称SOR迭代法）同样是基于矩阵分解原理，此时迭代法还有矩阵分解技术已经经历了100多年的发展，在SOR迭代法中的分裂矩阵引入了一个参数ω，其分裂矩阵为：M=1其中ω>0于是，可同上利用迭代法，构造的SOR迭代矩阵为Lω从而得到SOR迭代法.解线性方程组的SOR方法为x其中Lω=(D−ωL)−1由迭代法可构造公式(D−ωL)x(k+1)=((1-ω)D+ωU)x或Dx(k+1)=Dx(k+1)+ω(b+Lx(k+1)+由此，得到解求解Ax=b的SOR迭代法的计算公式为x3.3.2SOR迭代法的收敛条件SOR迭代法收敛的充分必要条件与Jacobi，Gauss迭代法一样，也是谱半径ρLω<1，根据算式可以看出SOR迭代法的谱半径此时有关于松弛因子ω的SOR法的收敛条件有定理如下：定理3.6：设解线性方程组Ax=b的SOR迭代法收敛，则0<ω<2[6]定理3.7：对于线性方程组Ax=b，若满足：A是一个对称正定矩阵，且A=D-L-U0<ω<2那么解Ax=b的SOR迭代法是收敛的[6]定理3.8：对于线性方程组Ax=b，如果：A是一个严格对角占优矩阵（或弱对角占优不可约矩阵）0<ω那么解Ax=b的SOR迭代法是收敛的[7]SOR迭代法的收敛速度与松弛因子ω有关，选取不同的ω，迭代次数可能相差非常大.而ω的选取规律至今没有很好的理论.3.4CG方法出现于20世纪的共轭梯度法（简称CG方法），既是一种用于解决大型线性方程组的算法，又是处理大型非线性方程组的比较有效的算法.因为仅需利用一阶导数信息，它既避免了最速下降法收敛性较差的缺点，又避免了牛顿法计算存储量大的缺点.相较于传统迭代法，它具有存储量小，步收敛性，稳定性高的优点，而且避免的其它参数的设置和计算.它是一种对应于求一个二次函数极值的变分方法，其仍然选择一组搜索方向p(0),p考虑线性方程组Ax=b的求解问题，其中A是给定的n维向量.定义一个二次函数φx可以得知：φ的梯度为φx=Ax-其中φ的Hessian矩阵就是矩阵A，因此，φ有唯一的极小点x∗=A−1b，从而有Ax根据上式可以看出求解方程组就对应于求解二次泛函数φ的极小点其算法具体为任取x(0)∈Rn,计算对k=0，1，…,计算αxr(k+1)=p（3）当r(k)=0，或（p(k)，Ap(k)）=0，计算停止，此时则x(k)=x∗.由于A正定，故当（p(k)，A由于{r(k)}互相正交，故而在r(0),r(1),…,r(n)中至少有一个零向量.若r(k)=0，则x(k)=x3.5PCG方法共轭梯度法的收敛速度与特征值的分布有一定关系，如果系数矩阵的特征值较均匀地分布在一个长区间上时，此时共轭梯度法的收敛性就会变得较差.而如果我们能令系数矩阵的特征值分布在一个较小的区间内，则此时算法的收敛性就比较好.因此，带有预条件的预处理共轭梯度算法诞生了.具体的预优共轭梯度算法如下(1)输入A，M，b和r0=b−Ax0,ρ0=r(2)ω=Apxk=xzk=Mβk=ρ若有ρk<ρ0ε,则输出xk这就是预处理的共轭梯度算法，也简称为PCG算法.由上式可以得知PCG算法的迭代次数与预矩阵M息息相关.要想迭代次数尽可能的少，预矩阵M的选取很关键，满足下列条件时，预矩阵M即为一个优异的预矩阵：M对称正定；M的稀疏性与A应大致相同；M−1具有某种特殊结构，如块对角，或三角矩阵的乘积等[2]3.6GMRES方法下面将介绍的一种算法是较为现代的迭代算法——GMRES算法，该算法属于空间逼近法.GMRES算法拥有非常优异的高效性和稳定性，常常被用于求解大型和对称不定矩阵和非对称线性方程组.具体算法如下：输入A，b，初始向量x0，最大迭代次数mr0=b−Ax0,ℎqj+1ℎj+1,j若满足ℎj+1,j<ε或j=m进行步骤qj+1对最小二乘问题进行求解min{Hjy若满足Axj−b2<ε这种算法不会发生中断，而且具有节约存储量和运算量少等优点。而m取多大比较好，至今还没有理论上的结果。理论仅有：对于任意的m，GMRES算法具有收敛性，若满足在系数矩阵A有正定的对称部分（即12（AT4数值算例分析算例4.1考虑下列二维Poisson方程边值问题−其中τ={(x,y)|0≤x≤1,0≤y≤1},利用分片线性连续有限元方法离散上述方程可得到相应的迭代方程组A其中h为网格步长，Aℎ为对称正定矩阵.由有限元方法的性质可知，Aℎ为稀疏矩阵，为了确保解的精度.需要h足够小，即为A表4.1不同网格步长下Jacobi,Gauss,CG法的迭代步数不难看出雅可比迭代法的收敛次数是最大的，故而现在实际生活中几乎不怎么运用，同样高斯-塞德尔迭代法的迭代次数也太大，对于电脑来说过于占用存储空间，所以也很少应用。而共轭梯度算法对于这两种算法来说，收敛速度要明显快很多.其次，我们来探讨超松弛迭代法以及松弛因子ω对迭代结果的影响.下面是ω取不同值时的迭代次数(这里的四个矩阵与上面的四个矩阵相同).表4.2不同网格步长下选取不同松弛因子的SOR法的迭代步数我们可以很容易发现对于这个步长为12其中松弛因子为1.9时，所选的256阶矩阵对应的迭代次数比64阶矩阵对应的迭代次数还要少，可见松弛因子对于SOR迭代法的影响还是非常大的.此外，实际过程中最佳松弛因子的选择只能通过实际试验得到.在这个算例的最后，我们来探讨预处理的共轭梯度算法.上面我们已经讨论过了共轭梯度算法，发现它的迭代速度要比雅可比迭代法，高斯-塞德尔迭代法快很多.尽管如此，由于实际计算时舍入误差的积累以及rk在这里我取了三个预条件，分别为：B为A的对角阵，B为A的三对角阵以及B为A的上三角阵，同样的依旧是上面那四个矩阵，下面我们来看对应的迭代次数.h预条件h预条件1111B为A的对角阵2653104205B为A的三对角阵2959117234B为A的上三角阵29不收敛不收敛不收敛表4.3不同网格步长下选取不同预条件时PCG法的迭代步数由表格可知B为A的对角阵是一个较优的预条件，此时的迭代次数要比没有预条件的共轭梯度算法少很多，而B为A的上三角阵则在不适合作为预条件.B为A的三对角阵虽然也可以作为预条件但不如B为A的对角阵作为预条件那么好.而在实际生活中，像高压断路器电场，水坝应力分析，三维电阻力探测等领域，预处理的共轭梯度算法及其改进算法都是很好的选择，当然在选择了好的预条件的前提下.算例4.2考虑下列二维Helmholfz方程−其中x利用分片线性连续有限元方法离散上述方程，令步长h为125此