Conjugate gradient method

數值分析在解線性方程最美麗的演算法莫過於 Conjugate gradient method,簡稱 CG method。這個演算法在假設矩陣 \(A \in M_n(R) \) 是對稱正定的前提下,可以在最多 n 步的演算中將 \(Ax=b\) 的解給解出。CG method 在結構上也非常美麗,第一、每一次遞迴都增加一個新的垂直方向的基底,因此當完整的基底被建構出來之後,自然就求得答案。第二、每一次的遞迴誤差就縮小一次,因此有機會在n 步之內就解到期待中的逼近解。以下是詳細地 CG method 的介紹。

再次假設矩陣 \(A\) 是一個對稱正定的矩陣,雖然一般的線性問題 \(Ax=b\) 當中的矩陣 \(A\) 未必是對稱正定,但是在 \(A\) 為可逆的情況下, \(A^TAx=A^Tb\) 當中的 \(A^TA\) 是滿足對稱正定的條件。因此,在不考慮計算量的情況下,CG method 還不失一般性。由於 \(A\) 是對稱正定矩陣,我們可以透過 \(A\) 在 \(R^n\) 上定義一個新的內積 \(<u,v>_A= u^T Av \)。

若對應這內積 \(<\cdot, \cdot>_A\) 存在某一個正交基底 \(\{p_i\}\) ,假設 \(x^*\) 為 \(Ax=b\) 的解,我們可以把 \(x^*\) 展開成 \(x^* = \sum \alpha_i p_i \) ,因此 \( b = Ax^* = \sum \alpha_i A p_i\)。利用 \(\{p_i\}\) 是正交基底的特性,我們有 \( p_k^Tb = \alpha_k p_k^TAp_k \) 的關係式,這個方程式可以讓我們解出係數 \(\alpha_i\) ,當有了\(\alpha_i\)之後我們就有解 \(x^*\)。

接下來的問題是如何找到正交基底 \(\{p_i\}\),這答案很顯然離不開 Gram-Schmidt 的正交演算法。上述的求解過程並沒有對 \(\{p_i\}\) 有過多的限制,因此選擇 \(\{p_i\}\) 的自由度還蠻高的。我們可以任意選定一個 \(p_0\) 然後產生新的亂數向量,利用 Gram-Schmidt 演算法依序把 \(\{p_i\}\) 建構出來,有了\(\{p_i\}\) 之後 \(Ax=b\) 的解就可以解出來了。但是 Gram-Schmidt 是一個很貴的演算法,如果能把計算 \(\{p_i\}\) 的過程稍加修改順帶算出最後的解,並且控制誤差,這不是兩全其美?Conjugate gradient 的設計正是如此。

例如,令 \(p_0 = r_0 = b - A x_0\) ,這裡的 \(x_0\) 是猜測解的起始值,許多人會直接令 \(x_0=0\)。 有了 \(p_0\) 就可以先計算 \(\alpha_0\) ,\(\alpha_0 = \frac{p_0^Tb}{p_0^TAp_0}\)。既然選擇 \(p_1\) 時用Gram-Schmidt 演算法我們可以選取一個亂數向量,這麼大的自由度不好好利用也實在浪費,因此,我們可以期待有 \(x_{k+1} = x_k + \alpha_k p_k\) 的結構。這樣由 \(x_0=0\) 出發的演算法,可以保證在 \(n\) 個步驟之後一定結束。有了新的 \(x_{k+1}\) 就可以決定新的殘量( residue)\(r_{k+1} = b-Ax_{k+1}\) ,再利用這個新的殘量,利用 Gram-Schmidt 演算法,計算正交於其他現有 \(p_i\) 的 \(p_{k+1}\)。重複這個步驟,得到 \(\alpha_{k+1}\) 更新新的 \(x_k\) ... 直到算出最後的答案。由於每次更新 \(x_{k+1}\) 時就多一個線性獨立的方向 \(p_k\) 當然對應的殘量 \(r_{k+1}\) 也會隨著遞迴步驟增加而減少。

實際上 CG method 在計算 \(p_{k+1}\) 時並不會真的用 Gram-Schmidt 方法 將 \(r_{k+1}\) 扣除現有 \(p_i\) 的方式來作,因為這樣的算法成本很高。我們可以很簡單地證明 \(p_{k+1} = r_{k+1} + \beta_k p_k\) ,其中 \(\beta_k = \frac{r_{k+1}^Tr_{k+1}}{r_k^T r_k}\)。從上述得知,CG method 的結構真的非常優美,也充分利用到各種縮減計算量的好處。在我們已知 \(x\) 的近似解時,我們會令 \(p_0=x_0\),這樣依照同樣的演算方法,計算到殘量夠小時便停止,就有機會在 \(n\) 步之內就將答案給解出來。