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

2、算法框架。1.線性方程組迭代解法與Krylov子空間我們考慮迭代法求解線性方程組Ax=b。假定未采用預處理矩陣P，或P矩 陣已經(jīng)隱含在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、得 到。以上發(fā)現(xiàn)自然與Krylov子空間的概念相聯(lián)系起來。Krylov 矩陣：K. = b Ab 血b Aj_1bKrylov子空間：K = b，Ab，Aj-ib的所有線性組合Krylov命名了向量b，Ab，A-ib的全部線性組合構(gòu)成的子空間，并認為在這一子空間中，有比上例中特定元素更與線性方程組的解相接近的元素。共 軛梯度法就是在這一子空間中，每一步迭代都依照某種標準尋求最優(yōu)元素的線性 方程組解法。對于每一步的余項rk=b-Axk,這一標準的不同就衍生出不同的解法:. K.中選取Xj，使得余項rk與子空間K.正交= Conjugate Gradients. %中選取，使得余項rk有最小范數(shù)=

4、 GMRES和MINRES. K.中選取 Xj，使得 rk 與子空間 Kj(AT)正交= BiConjugate Gradients2.共軛梯度法與Krylov子空間通常來說，Krylov基底指其經(jīng)過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+的倍數(shù)。 rk和qk+1的不同僅在于qk+1是規(guī)范化的。這就是共

5、軛梯度法與Krylov子空間的隱 含聯(lián)系。基于這一點，基向量q互相正交即可得出余項r也應互相正交。也即， TOC o 1-5 h z 對 i k, riTrk = 0(5)類似的，rk-1是qk的倍數(shù)，那么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.共軛梯度法與優(yōu)

6、化問題為說明共軛梯度法與優(yōu)化問題的關(guān)系，以下首先偽碼形式給出共軛梯度算法。% 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，定義能量函數(shù)E(x) = 1/2xtAx - xTb，則當A是正定 陣時，最小化能量函數(shù)E(x)等價于求解Ax=b。共軛梯度法在漸次增大的Krylov

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