1、Krylov子空間、優化問題與共軛梯度法自動化富曉鵬工程實踐中經常需要求解大型線性系統KU=F。在很多情況下矩陣K是非常 稀疏的，比如來自偏微分方程的離散化等，此時矩陣中每行僅有較少的非零元素。 面臨這樣的問題，我們首先面對的問題是，應該采用直接消元法還是迭代方法。 對前者來說，為充分利用系數特性，節點重編號是重要的；而對后者來說，適當 的預處理是關鍵。本文將重點放在后一類方法中的一種進行介紹與分析，即共軛梯度法。共軛 梯度法適用于矩陣K為對稱陣的情況，算法本身簡潔高效，且與一些其他的數 學理論、概念相緊密聯系，本文分析了共軛梯度法與Krylov子空間，以及優化 問題之間隱含的聯系，并簡要給出

2、算法框架。1.線性方程組迭代解法與Krylov子空間我們考慮迭代法求解線性方程組Ax=b。假定未采用預處理矩陣P，或P矩 陣已經隱含在A與b中。迭代法求解格式如下： TOC o 1-5 h z P - x = (P 一 A) - x + b(1)為說明問題，我們考慮簡單的迭代格式P=I，并且x1=b。則迭代的最初幾步 為：x = (I - A)b + b = 2b - Ab(2)x3 = (I - A)x2 + b = 3b - 3Ab + A2b(3) 由上面幾個式子可得，以上迭代格式第j步的解Xj是b，Ab，Aj-Ib的線 性組合。當A矩陣稀疏時，這些向量可以采用矩陣向量乘法的稀疏技巧很快

3、得 到。以上發現自然與Krylov子空間的概念相聯系起來。Krylov 矩陣：K. = b Ab 血b Aj_1bKrylov子空間：K = b，Ab，Aj-ib的所有線性組合Krylov命名了向量b，Ab，A-ib的全部線性組合構成的子空間，并認為在這一子空間中，有比上例中特定元素更與線性方程組的解相接近的元素。共 軛梯度法就是在這一子空間中，每一步迭代都依照某種標準尋求最優元素的線性 方程組解法。對于每一步的余項rk=b-Axk,這一標準的不同就衍生出不同的解法:. K.中選取Xj，使得余項rk與子空間K.正交= Conjugate Gradients. %中選取，使得余項rk有最小范數=

4、 GMRES和MINRES. K.中選取 Xj，使得 rk 與子空間 Kj(AT)正交= BiConjugate Gradients2.共軛梯度法與Krylov子空間通常來說，Krylov基底指其經過Arnoldi算法生成的正交基底qq?.,%。在Arnoldi算法中，新的基向量由t=Aqj-1與基向量q1,.,qj-1得到。也即，對 i = 1, 2,j-1, q.Tq.=0(4)對于第k步迭代，由于b屬于子空間K，xk屬于子空間Kk，余項rk=b-Axk 應該在子空間Kk+1中。又由于共軛梯度中rk與Kk正交，rk一定是qk+的倍數。 rk和qk+1的不同僅在于qk+1是規范化的。這就是共

5、軛梯度法與Krylov子空間的隱 含聯系。基于這一點，基向量q互相正交即可得出余項r也應互相正交。也即， TOC o 1-5 h z 對 i k, riTrk = 0(5)類似的，rk-1是qk的倍數，那么dr=rk-rk-1應該與更早的子空間Ki相正交，ik。又因為dx=xi-xi-1屬于子空間虬，那么，dxTdr=0o 將 dr展開，得rk - rk-1 = (b - Axk)- (b AxJ = -A(xk - Xk)(6)那么，x的更新量是“A正交”的或“共軛”的，(xi - xi.1)T A 伉-xk-1)= 0, for i k(7)這也就是共軛梯度法的名字由來。3.共軛梯度法與優

6、化問題為說明共軛梯度法與優化問題的關系，以下首先偽碼形式給出共軛梯度算法。% for linear equation Ax = bd0 = r0 = b, x0 = 0for k = 1,2,., k7 77 maxak=rk-1Trk-1/dk-1TAdk-1(8)xk=xk-1 +ak dk-1(9)rk=rk-1 -akAdk-1(10)Pk=rkTrk/rk-1Trk-1(11)dk=rk +Pk dk-1(12)對于線性方程組Ax=b，定義能量函數E(x) = 1/2xtAx - xTb，則當A是正定 陣時，最小化能量函數E(x)等價于求解Ax=b。共軛梯度法在漸次增大的Krylov

7、 子空間中搜索E(x)的最小值。第一個子空間K1是沿方向d0=b的直線，在這條直線上x=ab使能量函數最 小，即得到 a1： E(ab) = 1/2a2bTAb - abTb 在 a1 = bTb/bTAb 處取極值。這一 a1 就 是共軛梯度法偽碼式(8)中的值。能量函數E(x) = 1/2xtAx - xTb的梯度為Ax-b，若使用梯度下降法，則使用 負梯度方向rk = b-Axk即可。然而梯度下降法局部收斂性能較好，全局收斂性能 卻較差。上述偽碼采用的(12)的更新方法是保證搜索方向d1 = r1 +p1 d0與原方向 依式(7)的含義“A正交”。依照更新方向和更新步長，即可得到下一步的值xk = xk-1 +ak dk-1。以此迭代求解，直到滿足停止條