大型矩陣迭代

1、蘇州大學(xué)本科生畢業(yè)設(shè)計（論文）目錄 摘要 31 前言 52 迭代法的思想 62.1 迭代法的基本概念 62.2 迭代法的收斂條件 63 常見迭代法 73.1 Jacobi迭代法 73.1.1 Jacobi迭代法的計算公式 73.1.2 Jacobi迭代法的收斂條件 83.2 Gauss迭代法 93.2.1 Gauss迭代法的計算公式 93.2.2 Gauss迭代法的收斂條件 103.3 SOR迭代法 103.3.1 SOR迭代法的計算公式 103.3.2 SOR迭代法的收斂條件 113.4 CG迭代法 113.5 PCG迭代法 123.6 GMRES迭代法 134 數(shù)值算例分析 13算例4.1

2、 13算例4.2 145 總結(jié) 16參考文獻 17致謝 18摘要這篇論文介紹了大型矩陣運算時要用到的迭代法和其思想，并簡單介紹了常用的幾種迭代方法：如基于矩陣分解原理的Jacobi迭代法，Gauss迭代法和SOR迭代法；屬于共軛方向法的CG法和PCG法；最后介紹了適用于解決非線性方程組的GMRES迭代法及他們其中一些迭代法的收斂條件.并通過對不同矩陣的具體運算比較了他們的迭代次數(shù)，簡單分析了幾種迭代方法的優(yōu)劣.關(guān)鍵字:迭代法思想，Jacobi迭代法，Gauss迭代法，SOR迭代法,CG迭代法,PCG迭代法，GMRES迭代法，迭代次數(shù)，收斂體條件.AbstractThis essay intro

3、duces the iterative method which is needed during large matrix operation. Meanwhile, some of common methods of iteration are also suggested in this essay, for instance: Jacobian interation, Gauss interation and successive over relaxation method which were all based on the Matrix decomposition techni

4、que; conjugate gradient algorithm and Preprocessing conjugate gradient algorithm which both belong to conjugate direction method; GMRES algorithm which is used to solve non-linear equation and some of their own condition of convergence. In addition, by comparing their iteration steps, I roughly anal

5、ysed their merits and demerits.Keywords: iterative method，Jacobian interation, Gauss interation, successive over relaxation method, conjugate gradient algorithm, Preprocessing conjugate gradient algorithm, GMRES algorithm, iteration steps and condition of convergence.1 前言矩陣自19世紀被英國數(shù)學(xué)家凱利提出來，歷經(jīng)了200多年的

6、發(fā)展，在現(xiàn)代生活中有及其廣泛的應(yīng)用，下面我將列舉幾項矩陣在現(xiàn)實生活中的應(yīng)用：矩陣可用于預(yù)測水質(zhì)量狀況以及用于處理水污染等問題;矩陣還可以用于模擬飛機飛行所需要的飛行環(huán)境;在物理學(xué)領(lǐng)域的磁電場的探查測量方面也有很廣泛的應(yīng)用;許多圖片在電腦上的存儲方式就是以矩陣的格式存儲的,以便對圖片進行調(diào)整和美化的時候可以利用數(shù)學(xué)公式對圖片進行處理等。除此以外，矩陣還在許多領(lǐng)域有著極為廣泛的應(yīng)用，這里就不一一舉例了。而這些矩陣通常都是大型矩陣，如何高效的處理大型矩陣和與其相關(guān)的線性方程組成了一個值得深思的問題。而迭代法在處理零元素較多的大型稀疏矩陣時具有運算和存儲兩方面的優(yōu)勢。比較典型的有基于矩陣分解原理的Ja

7、cobi迭代法，這是由偉大的普魯士數(shù)學(xué)家雅可比在19世紀提出來的，在當(dāng)時Jacobi迭代法具有許多優(yōu)點，其計算公式簡單，每一次的迭代中只需計算一次矩陣和向量的乘法。這也是最早期的迭代法。同雅可比迭代法一樣，高斯-塞德爾和超松弛迭代法也是基于矩陣分解原理，這些都是傳統(tǒng)的迭代法，出現(xiàn)時間較久。此外，還有其他一些優(yōu)異的矩陣迭代算法，如以極小化的方法來討論方程組的解的著名的共軛梯度（CG）算法，CG法出現(xiàn)在20世紀中期，它有許多優(yōu)異的性質(zhì)：比如相較于傳統(tǒng)迭代法，CG法收斂速度較快；此外CG法只有極少量的非向量運算。然而由于實際運算中舍入誤差的積累以及正交性的逐漸喪失，CG法長時間并未得到廣泛應(yīng)用。直到

8、20世紀70年代，預(yù)處理共軛梯度算法（PCG算法）出現(xiàn)了，它通過對CG算法中未加工過的矩陣進行某種形式的加工改造，使它具有更適宜迭代的特點。近些來，PCG法有了進一步的發(fā)展。它現(xiàn)在已經(jīng)有了MICCG，塊預(yù)處理，高階LU分解和行和一致預(yù)處理等預(yù)處理方法，預(yù)處理技術(shù)正在不斷完善中；此外共軛梯度算法也發(fā)展為可應(yīng)用于不定陣，非奇異陣及復(fù)矩陣等多種矩陣形式的較為系統(tǒng)的方法。經(jīng)過了近200年的發(fā)展，迭代法變得越來越系統(tǒng)，且新的優(yōu)異的迭代法也不斷涌出，本文主要對其中幾種比較典型的迭代方法進行介紹和探究，本文的內(nèi)容主要分為：迭代法基本概念及收斂性的介紹，幾種迭代法步驟及特性的介紹(Jacobi,Gauss,S

9、OR,CG,PCG,GMRES)，具體的算例分析和總結(jié)。2 迭代法的思想2.1 迭代法的基本概念迭代法的由來與線性方程組的求解有較深的關(guān)系，下面我們來舉例說明有線性方程組Ax=b, （2.1）其中A為大型稀疏矩陣（具有較多的零元素），此時直接求解很困難，而用迭代法是比較合理有效的.例如：求解線性方程組8x13x2+2x3=20,4x1+11x2x3=33,6x1+3x2+12x3=36, （2.2）記為Ax=b，其中有A=83241116312, x=x1x2x3, b=203336.該方程組的精確解為x=3，2，1T.線性方程組（2.2）可以恒等變形為下式x1=183x22x3+20,x2=

10、1114x1+x3+33,x3=1126x13x2+36; (2.3)則有x=B0x+f，其中B0=0382841101116123120， f=82233113612.任意選取初始值，例如x(0)=0，0，0T，代入（2.3）式可以得到新的值x(1)=x1(1)，x2(1)，x3(1)T，再將新的x(1)分量代入得到x(2)，反復(fù)利用這個計算程序，即可得到迭代法的一般式x(k+1)=B0 x(k)+f，（k為迭代次數(shù)）.除（2.3）外，線性方程組（2.2）的恒等變形方式還有很多種，不同的恒等變形可得到不同的迭代式，均可寫成如下形式x(k+1)=B0 x(k)+f然而向量序列x(k)并非都可以

11、逐步逼近方程組的解x.2.2 迭代法的收斂條件不難看出若有l(wèi)im kx(k)=x存在，則此迭代法是收斂的（ x為上面方程組的解），否則該迭代法即為發(fā)散的.由上可知，若要研究x(k)的收斂性，可以引入一個誤差向量(k+1)=x(k+1)x,對于線性方程組x=Bx+f,設(shè)有唯一解為x，則x=Bx+f.與x(k+1)=B0 x(k)+f 相減得， (k+1)=B(k)(k=0,1,2),遞推得(k)=B(k1)=Bk0.要考察x(k)的收斂性，就要研究B在什么條件下limk(k)=0 ，此時有定理：迭代公式收斂的充要條件是矩陣的譜半徑B1 13 常見的幾種迭代法3.1 Jacobi迭代法3.1.1

12、Jacobi迭代法的計算公式上面迭代法的思想是從線性方程組的角度來看，我們不妨從矩陣的角度來看.矩陣分解技術(shù)在矩陣中具有廣泛應(yīng)用，它通過對矩陣的分解來構(gòu)造迭代格式.如將矩陣A分解為M+N故此Ax=b即可寫成Mx+Nx=b即Mx =b-Nx當(dāng)M可逆時有x=M1bM1Nx由此可構(gòu)造迭代格式x(k+1)=M1Nxk+M1b其中M1N為迭代矩陣，且M要滿足是可逆的且逼近A.下面我們將討論幾種基于矩陣分解原理的常見迭代法.將線性方程組中的系數(shù)A=(aij)Rnn分成三部分A=a11 a22 ann_0 a210 an1,1an1,20 an1an2an,n10_0a12a1,n1a1n 0a2,n1a2

13、n 0an1,n 0:=D-L-U.設(shè)aii0(i=1,2,n),選取M為A的對角元素部分，即選取 M= D（對角矩陣），A= D-N，即可得到Jacobi迭代法x(0), 初始向量,x(k+1)=Bx(k)+f,k=0,1,其中B=I-D1A=D1（L+U）BJ，f=D1b.其中BJ為雅可比迭代法（Jacobi）的迭代矩陣.由迭代法有公式Dx(k+1)=（L+U） x(k)+b或aiixi(k+1)=j=1i1aiixik+bi, i =1,2,n.因此，對應(yīng)于求解線性方程組Ax=b的雅可比迭代法的計算公式為x(0)=(x10,x20,xn0)T,xi(k+1)=(bij=1jinaijxj

14、k)/aiii=1,2,n;k=0,1,表示迭代次數(shù),雅可比迭代法計算形式比較簡單，每一次只需將算出來的值再代入迭代算式里，循環(huán)反復(fù)即可得到結(jié)果，但也是因為這個原因，迭代次數(shù)往往較多，故在實際生活中已經(jīng)很少被應(yīng)用.3.1.2 雅可比迭代法的收斂性下面我們來討論雅可比迭代法的收斂性.由上面介紹的迭代法的收斂條件可知迭代法的收斂性于它自身的譜半徑有關(guān)，有定理如下定理3.1：設(shè)Ax=b，其中A=D-L-U為非奇異矩陣，且對角矩陣D也非奇異，則解線性方程組雅可比迭代法收斂的充要條件是BJj=1jinaij,i=1,2,n,稱A為嚴格對角占優(yōu)矩陣.(2) 若A的元素滿足aijj=1jinaij,i=1,

15、2,n,且上式至少有一個不等式成立，則稱A為弱對角占優(yōu)矩陣1(可約與不可約矩陣) 設(shè)A=(aij)nn(n2),如果存在置換矩陣矩陣P使PTAP=A11A120A22,其中A11為r階方程，A22為n-r階方程（1r0(i=1,2,n),則解線性方程組Ax=b的雅可比迭代法收斂的充分必要條件是A及2D-A均為正定矩陣，其中D=diag(a11,a22,ann)13.2 Gauss迭代法3.2.1 Gauss迭代法的計算公式類似的，Gauss迭代法也是根據(jù)矩陣分解技術(shù)，將A的下三角部分作為分裂矩陣M，即M= D- L（下三角矩陣），A=M-N，這樣就可以得到Ax=b的Gauss迭代法x(0),

16、初始向量,x(k+1)=Bx(k)+f,k=0,1,其中B=I-（DL)1A=(DL)1UBGs，f=(DL)1b.其中BGs為高斯-塞德爾迭代法（Gauss）的迭代矩陣.由迭代法可構(gòu)造Gauss迭代公式(DL) x(k+1)=U x(k)+b，或?qū)憺镈x(k+1)=L x(k+1)+ U x(k)+b，即aiixi(k+1)=bij=1i1aijxjk+1j=i+1naijxjk,i=1,2,n.于是，高斯-塞德爾迭代法的計算公式為x(0)=(x10,x20,xn0)T,xi(k+1)=(bij=1 i1aijxjk+1j=i+1naijxj(k)/aiii=1,2,n;k=0,1,表示迭代

17、次數(shù),同Jacobi迭代法一樣，Gauss迭代法也是比較傳統(tǒng)的迭代方法，出現(xiàn)了近200年，早已不適合實際運用.3.2.2 高斯-塞德爾迭代的收斂性與雅可比迭代法一樣，解線性方程組的高斯-塞德爾迭代法收斂的充要條件是譜半徑BGs0(i=1,2,n),則解線性方程組Ax=b的高斯-塞德爾迭代法收斂的充分必要條件是A正定53.3 SOR迭代法3.3.1 SOR迭代法的計算公式出現(xiàn)于20世紀70年代的超松弛迭代法（簡稱SOR迭代法）同樣是基于矩陣分解原理，此時迭代法還有矩陣分解技術(shù)已經(jīng)經(jīng)歷了100多年的發(fā)展，在SOR迭代法中的分裂矩陣引入了一個參數(shù)，其分裂矩陣為：M=1(DL),其中0為可選擇的松弛因

18、子.于是，可同上利用迭代法，構(gòu)造的SOR迭代矩陣為LI(DL)1A=(DL)1(1D+U).從而得到SOR迭代法.解線性方程組的SOR方法為x(0), 初始向量,x(k+1)=Lx(k)+f,k=0,1,其中L=(DL)1(1D+U)，f=(DL)1b.由迭代法可構(gòu)造公式(DL) x(k+1)=(1-)D+U) x(k)+b,或Dx(k+1)=D x(k+1)+ (b+Lx(k+1)+U x(k)Dx(k).由此，得到解求解Ax=b的SOR迭代法的計算公式為x(0)=(x10,x20,xn0)T,xi(k+1)=xi(k)+(bij=1 i1aijxjk+1j=i+1naijxj(k)/aii

19、i=1,2,n;k=0,1,表示迭代次數(shù),為松弛因子，3.3.2 SOR迭代法的收斂條件SOR迭代法收斂的充分必要條件與Jacobi，Gauss迭代法一樣，也是譜半徑L1，根據(jù)算式可以看出SOR迭代法的譜半徑L還和松弛因子有關(guān)此時有關(guān)于松弛因子的SOR法的收斂條件有定理如下：定理3.6：設(shè)解線性方程組Ax=b的SOR迭代法收斂，則026定理3.7：對于線性方程組Ax=b，若滿足：（1） A是一個對稱正定矩陣，且A=D-L-U（2） 02那么解Ax=b的SOR迭代法是收斂的6定理3.8：對于線性方程組Ax=b，如果：（1） A是一個嚴格對角占優(yōu)矩陣（或弱對角占優(yōu)不可約矩陣）（2） 01那么解Ax

20、=b的SOR迭代法是收斂的7SOR迭代法的收斂速度與松弛因子有關(guān)，選取不同的，迭代次數(shù)可能相差非常大.而的選取規(guī)律至今沒有很好的理論.3.4 CG方法出現(xiàn)于20世紀的共軛梯度法（簡稱CG方法），既是一種用于解決大型線性方程組的算法，又是處理大型非線性方程組的比較有效的算法. 因為僅需利用一階導(dǎo)數(shù)信息，它既避免了最速下降法收斂性較差的缺點，又避免了牛頓法計算存儲量大的缺點.相較于傳統(tǒng)迭代法，它具有存儲量小，步收斂性，穩(wěn)定性高的優(yōu)點，而且避免的其它參數(shù)的設(shè)置和計算.它是一種對應(yīng)于求一個二次函數(shù)極值的變分方法，其仍然選擇一組搜索方向p(0),p(1), .但它不再是具有正交性的r(0),r(1),

21、方向.考慮線性方程組Ax=b的求解問題，其中A是給定的n維向量.定義一個二次函數(shù)x=12xTAxbTx.可以得知：的梯度為x=Ax-b，其中的Hessian矩陣就是矩陣A，因此，有唯一的極小點x=A1b，從而有Ax=b即x=min x.根據(jù)上式可以看出求解方程組就對應(yīng)于求解二次泛函數(shù)的極小點其算法具體為（1） 任取x(0)Rn,計算r(0)=bAx(0),取p(0)=r(0).（2） 對k=0，1，,計算k=(rk,r(k)(pk,Ap(k)x(k+1)=x(k)+kp(k)r(k+1)=r(k)kAp(k),k=(rk+1,r(k+1)(rk,r(k)p(k+1)=r(k+1)+kp(k)（

22、3）當(dāng)r(k)=0，或（p(k)，Ap(k)）=0，計算停止，此時則x(k)=x.由于A正定，故當(dāng)（p(k)，Ap(k)）=0時，p(k)=0，而(rk,r(k)= (rk,p(k)=0，即rk=0.由于r(k)互相正交，故而在r(0),r(1),r(n)中至少有一個零向量.若r(k)=0，則x(k)=x.所以用CG法求解n維線性方程組時，精確解至多n步即可求得，所以CG法也屬于直接法.但由于舍入誤差的存在，后面的r(k)正交性就很難得到保證33.5 PCG方法共軛梯度法的收斂速度與特征值的分布有一定關(guān)系，如果系數(shù)矩陣的特征值較均勻地分布在一個長區(qū)間上時，此時共軛梯度法的收斂性就會變得較差.而

23、如果我們能令系數(shù)矩陣的特征值分布在一個較小的區(qū)間內(nèi)，則此時算法的收斂性就比較好. 因此，帶有預(yù)條件的預(yù)處理共軛梯度算法誕生了.具體的預(yù)優(yōu)共軛梯度算法如下(1)輸入A，M，b和x0；r0=bAx0,z0=M1r0,p1=z0,0=r0Tz0,k=1.(2)=Apk,k=k1/pkT,xk=xk1+kpk,rk=rk1k,zk=M1rk,k=rkTzk,k=k/k1,k+1=zk+kpk.(3) 若有k0,則輸出xk，結(jié)束；如若不然則令k=k+1返回（2）這就是預(yù)處理的共軛梯度算法，也簡稱為PCG算法.由上式可以得知PCG算法的迭代次數(shù)與預(yù)矩陣M息息相關(guān).要想迭代次數(shù)盡可能的少，預(yù)矩陣M的選取很關(guān)

24、鍵，滿足下列條件時，預(yù)矩陣M即為一個優(yōu)異的預(yù)矩陣：(1) M對稱正定；(2) M的稀疏性與A應(yīng)大致相同；(3) M1A的特征值分布在一個較小的區(qū)間內(nèi)；(4) 具有某種特殊結(jié)構(gòu)，如塊對角，或三角矩陣的乘積等23.6 GMRES方法下面將介紹的一種算法是較為現(xiàn)代的迭代算法GMRES算法，該算法屬于空間逼近法. GMRES算法擁有非常優(yōu)異的高效性和穩(wěn)定性，常常被用于求解大型和對稱不定矩陣和非對稱線性方程組.具體算法如下：（1） 輸入A，b，初始向量x0，最大迭代次數(shù)m及精度要求.（2） r0=bAx0,q1=r0/r02,j=1.（3） ij=qiTAqj,i=1,2,j,qj+1=Aqji=1ji

25、jqi,j+1,j=qj+12.（4） 若滿足j+1,j或j=m，則轉(zhuǎn)至（5）；否則進行步驟qj+1=qj+1j+1,j， j=j+1.返回（3）（5） 對最小二乘問題進行求解min Hjyr0e12：yRj得yj，xj=x0+Qjyj.（6） 若滿足Axjb2,則將xj輸出，結(jié)束;反之，x0=xj，返回步驟（2）.這種算法不會發(fā)生中斷，而且具有節(jié)約存儲量和運算量少等優(yōu)點。而m取多大比較好，至今還沒有理論上的結(jié)果。理論僅有：對于任意的m， GMRES算法具有收斂性，若滿足在系數(shù)矩陣A有正定的對稱部分（即12（AT+A）正定24 數(shù)值算例分析算例4.1考慮下列二維Poisson方程邊值問題u=f

26、 在內(nèi)u=0 在TD上un=g 在TN上其中=(x,y)|0x1,0y1,TD=(x,y)|(x,y) ,y=0或x=1,TN=(x，y)|(x,y),y=1或x=0,f=1,g=1利用分片線性連續(xù)有限元方法離散上述方程可得到相應(yīng)的迭代方程組Ax=b其中h為網(wǎng)格步長，A為對稱正定矩陣.由有限元方法的性質(zhì)可知，A為稀疏矩陣，為了確保解的精度.需要h足夠小，即為A大型稀疏矩陣.此時用直接法求解工作量大.因此接下來，我們用本文的各種迭代法求解AX=b表4.1不同網(wǎng)格步長下Jacobi,Gauss,CG法的迭代步數(shù)不難看出雅可比迭代法的收斂次數(shù)是最大的，故而現(xiàn)在實際生活中幾乎不怎么運用，同樣高斯-塞德

27、爾迭代法的迭代次數(shù)也太大，對于電腦來說過于占用存儲空間，所以也很少應(yīng)用。而共軛梯度算法對于這兩種算法來說，收斂速度要明顯快很多.其次，我們來探討超松弛迭代法以及松弛因子對迭代結(jié)果的影響.下面是取不同值時的迭代次數(shù)(這里的四個矩陣與上面的四個矩陣相同).表4.2不同網(wǎng)格步長下選取不同松弛因子的SOR法的迭代步數(shù)我們可以很容易發(fā)現(xiàn)對于這個步長為123的矩陣，松弛因子取1.7時，迭代次數(shù)最少.而對于另外三個矩陣，都是松弛因子為1.9時迭代次數(shù)最少.其中松弛因子為1.9時，所選的256階矩陣對應(yīng)的迭代次數(shù)比64階矩陣對應(yīng)的迭代次數(shù)還要少，可見松弛因子對于SOR迭代法的影響還是非常大的.此外，實際過程中

28、最佳松弛因子的選擇只能通過實際試驗得到. 在這個算例的最后，我們來探討預(yù)處理的共軛梯度算法.上面我們已經(jīng)討論過了共軛梯度算法，發(fā)現(xiàn)它的迭代速度要比雅可比迭代法，高斯-塞德爾迭代法快很多.盡管如此，由于實際計算時舍入誤差的積累以及rk正交性的逐漸喪失，因此很難得到精確解，所以在實用中效率不高，而預(yù)處理的共軛梯度法的出現(xiàn)很好的解決了這一問題.下面我們來討論預(yù)處理的共軛梯度算法及不同預(yù)條件對迭代結(jié)果的影響.在這里我取了三個預(yù)條件，分別為：B為A的對角陣，B為A的三對角陣以及B為A的上三角陣，同樣的依舊是上面那四個矩陣，下面我們來看對應(yīng)的迭代次數(shù).h預(yù)條件123124125126B為A的對角陣2653

29、104205B為A的三對角陣2959117234B為A的上三角陣29不收斂不收斂不收斂表4.3不同網(wǎng)格步長下選取不同預(yù)條件時PCG法的迭代步數(shù)由表格可知B為A的對角陣是一個較優(yōu)的預(yù)條件，此時的迭代次數(shù)要比沒有預(yù)條件的共軛梯度算法少很多，而B為A的上三角陣則在不適合作為預(yù)條件.B為A的三對角陣雖然也可以作為預(yù)條件但不如B為A的對角陣作為預(yù)條件那么好.而在實際生活中，像高壓斷路器電場，水壩應(yīng)力分析，三維電阻力探測等領(lǐng)域，預(yù)處理的共軛梯度算法及其改進算法都是很好的選擇，當(dāng)然在選擇了好的預(yù)條件的前提下.算例4.2考慮下列二維Helmholfz方程uu=f 在內(nèi)un+u=g 在上其中x,y0x1,0y1，f=