共軛梯度法實(shí)驗(yàn)報(bào)告

1、本文格式為word版，下載可任意編輯共軛梯度法實(shí)驗(yàn)報(bào)告 數(shù)值代數(shù)試驗(yàn)報(bào)告 一、 試驗(yàn)名稱： 用共軛梯度法解線性方程組。 二、試驗(yàn)?zāi)康模?進(jìn)一步熟識(shí)理解把握共軛梯度法解法思路，提高 matlab 編程力量 。 三、 試驗(yàn)要求： 已知線性方程矩陣，應(yīng)用共軛梯度法在相關(guān)軟件編程求解線性方程 組的解。 四、試驗(yàn)原理： 1共軛梯度法： 考慮線性方程組 ax b 的求解問題，其中 a 是給定的 n 階對(duì)稱正定矩陣， b 是給定的 n 維向量， x 是待求解 的 n 維向量. 為此, 定義二次泛函 t t (x) x ax 2b x . 定理 1 設(shè) a 對(duì)稱正定 , 求方程組 ax b 的解，等價(jià)于求二次

2、泛函 (x) 的微小值點(diǎn) . 定理 1 表明，求解線性方程組問題就轉(zhuǎn)化為求二次泛函 (x) 的微小值點(diǎn)問題 . 求解二次函數(shù)微小值問題，通常似乎盲人下山那樣，先給定一個(gè)初始向量 x 0 ，確定 一個(gè)下山方向 p 0 ，沿著經(jīng)過點(diǎn) x 0 而方向?yàn)?p 0 的直線 x x 0 p 0 找一個(gè)點(diǎn) x x p ， 1 0 0 0 使得對(duì)全部實(shí)數(shù) 有 x p x p ， 0 0 0 0 0 即在這條直線上 x 1 使 (x) 達(dá)到微小. 然后從 x 1 動(dòng)身，再確定一個(gè)下山的方向 p 1 ，沿著 直 線 x x p 再 跨 出 一 步 ， 即 找 到 1 使 得 x 在 x 2 x 1 1 p 1 達(dá)

3、 到 極 小 ： 1 1 x p x p . 1 1 1 1 1 重復(fù)此步驟，得到一串 0 , 1 , 2 ,l 和 p 0 , p 1 , p 2 ,l ， 稱 p k 為 搜尋方向 ， k 為步長(zhǎng). 一般狀況下，先在 x k 點(diǎn)找下山方向 p k ，再在直線 x x p 上確定步長(zhǎng) k 使 k k x p x p k k k k k , 最終求出 x x p . 然而對(duì)不同的搜尋方向和步長(zhǎng)，得到各種不同的算法 . k 1 k k k 1 由此，先考慮如何確定 k . 設(shè)從 x k 動(dòng)身，已經(jīng)選定下山方向 p k . 令 f x p k k t t x p a x p 2b x p k k

4、k k k k 2 t 2 t p ap r p x ， k k k k k 其中 r k b ap k . 由一元函數(shù)極值存在的必要條件有 t t f 2 p ap 2r p 0 k k k k 所確定的 即為所求步長(zhǎng) ，即 k 步長(zhǎng)確定后，即可算出 t r p k k k t p ap k k . x x p . k 1 k k k t 此時(shí)，只要 r p 0 ，就有 k k 即 x x . k 1 k x x x p x 2 k 1 k k k k k t r p k k 2 t 2 t 0 p ap r p k k k k k k t p ap k k 再考慮如何確定下山方向 p k

5、. 易知負(fù)梯度方向是 (x) 減小最快的方向， 但簡(jiǎn)潔分 析就會(huì)發(fā)覺負(fù)梯度方向只是局部最佳的下山方向， 而從整體來(lái)看并非最佳 . 故采納新 的方法尋求更好的下山方向 共軛梯度法 . 下面給出共軛梯度法的詳細(xì)計(jì)算過程 : 給定初始向量 x 0 ，第一步仍選用負(fù)梯度方向?yàn)橄律椒较?，?p 0 r 0 ，于是有 t r r 0 0 ,x x p ,r b ax . 0 t 1 0 0 0 1 0 p ap 0 0 對(duì)以后各步， 例如第 k+1 步（k 1）, 下山方向不再取 r k ，而是在過點(diǎn)由向量 r k 和 p k 1 所張成的二維平面 2 x | x x k r k p k 1 , , r

6、內(nèi)找出訪函數(shù) 下降最快的方向作為新的下山方向 p k . 考慮 在 2 上的限制： , ( x k r k p k ) 1 t (x r p ) a(x r p ) k k k 1 k k k 1 t 2b (x r p ) . k k k 1 計(jì)算 關(guān)于 , 的偏導(dǎo)得： t t t 2 r ar r ap r r , k k k k 1 k k t t 2 r ap p ap , t 其中最終一式用到了 r p 1 0 ，這可由 r k 的定義直接驗(yàn)證 . 令 k k 1 k 1 k 1 k k 0 ， 即知 在 2 內(nèi)有唯一的微小值點(diǎn) x% x r p ， k 0 k 0 k 1 2 其中

7、 0 和 0 滿意 t t t r ar r ap r r 0 k k 0 k k 1 k k , t t r ap p ap 0 k k 1 0 k 1 k 1 0. 由于 r k 0 必有 0 0 ，所以可取 1 0 p x% x r p k k k k 1 0 0 作為新的下山方向 . 明顯，這是在平面 2 內(nèi)可得的最佳下山方向 . 令 0 k ，則可 1 得 t r ap k k 1 . k 1 t p ap k 1 k 1 t 注：這樣確定的 p k 滿意 p ap 1 0 ，即 p k 與 p k 1 是相互共軛的 . k k 0 總結(jié)上面的爭(zhēng)論，可得如下的計(jì)算公式： t r p

8、k k ， x k 1 x k k p k ， k t p ap k k r b ax ， k 1 k 1 t r ap k 1 k ， p k 1 r k 1 k p k . k t p ap k k 在實(shí)際計(jì)算中，常將上述公式進(jìn)一步簡(jiǎn)化，從而得到一個(gè)形式上更為簡(jiǎn)潔而且對(duì)稱 的計(jì)算公式 . 首先來(lái)簡(jiǎn)化 r k 1 的計(jì)算公式： r 1 b ax 1 b a(x p ) r ap . k k k k k k k k 由于 ap k 在計(jì)算 k 是已經(jīng)求出，所以計(jì)算 r k 1 時(shí)可以不必將 x k 1 代入方程計(jì)算，而是 從遞推關(guān)系 r k 1 b k ap k 得到. 再來(lái)簡(jiǎn)化 k 和 k

9、的計(jì)算公式 . 此處需要用到關(guān)系式 t t t r r 1 r p 1 r 1 p 0, k 1,2,k . k k k k k k 從而可導(dǎo)出 1 t t r r r , ， k 1 k 1 k 1 1 1 k t t t p ap p r r p r k k k k k 1 k k 1 t 1 t k k r r p r r . k k k 1 k 1 k k k k 由此可得 t t r r r r k k 1 1 . k k , . ， k t k t p ap r r k k k k 從而有求解對(duì)稱正定方程組的共軛梯度法算法如下： x 初值 0 r b ax ； k 0 0 0 wh

10、ile r k 0 k k 1 if k 1 p r 0 0 3 else t t k r k r k r k r k 2 1 1 2 2 p r p k 1 k 1 k 2 k 2 end t t k r k r k p k ap k 1 1 1 1 1 x x p k k 1 k 1 k 1 r r ap k k 1 k 1 k 1 end x x k 注：該算法每迭代一次僅需要使用系數(shù)矩陣 a 做一次矩陣向量積運(yùn)算 . 定理 2 由共軛梯度法得到的向量組 r 和 p i 具有如下基本性質(zhì)： i t （1） p r 0， 0 i j k; i j t （2） r r 0 ， i j ， 0

11、 i, j k; i j t （3） p ap 0， i j ， 0 i, j k; i j （4） span r ,k , r k span p ,k , p k (a, r ,k 1) ， 0 0 0 其中 k ( a,r , k 1) span r , ar ,k , a r ， 0 0 0 0 通常稱之為 krylov 子空間. 下面給出共軛梯度法全局最優(yōu)性定理： 定理 3 用共軛梯度法計(jì)算得到的近似解 x k 滿意 x min x : x x (a,r ,k) k 0 0 或 x x x x x x a r k ， * min * : 0 ( , 0 , ) k a a 其中 t x

12、 x ax ， x * 是方程組 ax b 的解， ( a, r 0 ,k) 是由所定義的 krylov 子空間. a 定理 2 表明，向量組 r 0 ,k ,r k 和 p 0 ,k , p k 分別是 krylov 子空間 ( a, r 0 ,k 1) 的正 交基和共軛正交基 . 由此可知，共軛梯度法最多 n 步便可得到方程組的解 x * . 因此， 理論上來(lái)講，共軛梯度法是直接法 . 然而實(shí)際使用時(shí)，由于誤差的消失，使 r k 之間的正交性很快損失，以致于其有 限步終止性已不再成立 . 此外，在實(shí)際應(yīng)用共軛梯度法時(shí)，由于一般 n 很大，以至于 迭代 o n 次所耗費(fèi)的計(jì)算時(shí)間就已經(jīng)使用戶

13、無(wú)法接受了 . 因此，實(shí)際上將共軛梯度 法作為一種迭代法使用， 而且通常是 r 是否已經(jīng)很小及迭代次數(shù)是否已經(jīng)達(dá)到最大 k 允許的迭代次數(shù) k max 來(lái)終止迭代 . 從而得到解對(duì)稱正定線性方程組的有用共軛梯度 法，其算法如下： x 初值 4 k 0; r b ax; t r r while b and k k max 2 k k 1 if k 1 p r else %; p r p end t ; ; ap p x x p t r r ; % ; r r end 算法中，系數(shù)矩陣 a 的作用僅僅是用來(lái)由已知向量 p 產(chǎn)生向量 ap ，這不僅可以 充分利用 a 的稀疏性，而且對(duì)某些供應(yīng)矩陣 a

14、 較為困難而由已知向量 p 產(chǎn)生向量 ap 又非常便利的應(yīng)用問題是非常有益的。 2算例 10 1 2 3 4 x 1 12 運(yùn)用共軛梯度法求解下述方程，并解釋你所觀看到的結(jié)果 1 9 -1 2 - 3 x 27 2 2 -1 7 3 - 5 x 3 14 3 2 3 12 -1 x 17 解：已知共軛梯度法的 matlab 程序代碼如下所示： 4 4 - 3 - 5 -1 15 x 12 functionx,n=conjgrad(a,b,x0) 5 %采納共軛梯度法求線性方程組 ax=b 的解 %線性方程組的系數(shù)矩陣： a %線性方程組中的常數(shù)向量： b %迭代初始向量： x0 %線性方程組的

15、解： x %求出所需精度的解實(shí)際的迭代步數(shù)： n r1=b-a*x0; p=r1; n=0; for i=1:rank(a) if(dot(p,a*p)1.0e-50) break; end alpha=dot(r1,r1)/dot(p,a*p); 5 x=x0+alpha*p; r2=r1-alpha*a*p; if(r21.0e-50) break; end belta=dot(r2,r2)/dot(r1,r1); p=r2+belta*p; n=n+1; end 用共軛梯度法求解，在 matlab 命令窗口中輸入求解程序： a=10 1 2 3 4 2 4 2 1 1;1 9 -1 2

16、-3 3 -3 -1 9 2;2 -1 7 3 -5 4 -5 7 -1 3;3 2 3 12 -1 5 -1 3 2 4;4 -3 -5 -1 15 1 3 2 3 2;2 3 4 5 1 10 4 5 3 4;4 -3 -5 -1 3 4 9 1 -3 2;2 -1 7 3 2 5 1 7 -1 1;1 9 -1 2 -3 3 -3 -1 12 10;1 2 3 4 2 4 2 1 10 15; b=12;-27;14;-17;12;12;-27;14;-17;12; x0=0;0;0;0;0;0;0;0;0;0; x,n=conjgrad(a,b,x0) 輸出的計(jì)算結(jié)果為： x = 8.

17、5669 -7.1981 -7.3479 -13.9388 1.1303 11.5188 -26.8966 -2.2388 0.0829 13.2786 輸出的迭代次數(shù)為 n = 10 共軛梯度法理論上只要迭代 n 步，就能得出方程組的解，但是由于存在計(jì)算誤差， 即分?jǐn)?shù)向小數(shù)轉(zhuǎn)化時(shí)存在舍入誤差，很難保證在第 n 步時(shí)得到精確解。 6 3 總結(jié) 本文首先給出最小二乘問題的定義，隨后從盲人下山法開頭爭(zhēng)論了共軛梯度法 的詳細(xì)推導(dǎo)過程及其相關(guān)性質(zhì)與算法 . 繼而重點(diǎn)給出正則化方法的有用共軛梯度算 法并舉例進(jìn)行檢驗(yàn) . 最終，需要說明雖然正則化方法是求一般線性方程組 ax b， m n a r m n 且 a 列滿秩的最小二乘解的一種方法且簡(jiǎn)潔易行，但是也有很多不 足之處，如 m n 時(shí)一般無(wú)解； t a a 形成時(shí)運(yùn)算量大， a 中某些信息會(huì)丟失；當(dāng) a病 態(tài)時(shí)其收斂性速度由于 t 2 2 (a a) 2 (a) 很大變得特別之慢等，故為了避開正則化方 法的缺點(diǎn)，還可運(yùn)用殘量微小化方法或殘量正交化方法等更好的方法來(lái)解決此類問 題. 在實(shí)際的科學(xué)與工程問題中，經(jīng)常將問題歸結(jié)為一個(gè)線性方程組的求解問