以前、有限体積法のところで、共役勾配法について説明しました。
共役勾配法には様々な種類がありますが、今回はいくつかの異なる共役勾配法について解説します。
目次
共役勾配法
共役勾配法（Conjugate Gradient Method）は線形方程式を解くための解法として広く用いられています。詳細は上記の記事を参照してください。おさらいのため再度アルゴリズムを書いておきます。
線形方程式系、
$${\bf A}{\bf x} = {\bf b} \tag{1}$$
に対して、共役勾配法のアルゴリズムは以下のようになります。
-
適当な初期値${\bf x}_0$を決める。
-
${\bf r}_0 = {\bf b} – {\bf A} {\bf x}_0$、${\bf d}_0 = {\bf r}_0$とする。
-
収束するまで以下の計算を繰り返す。① $\alpha_k=({\bf d}_k, {\bf r}_k)/({\bf d}_k, {\bf A}{\bf d}_k) $② ${\bf x}_{k+1} = {\bf x}_k + \alpha_k {\bf d}_k$③ ${\bf r}_{k+1} = {\bf r}_k – \alpha_k{\bf A} {\bf d}_k$④ 残差 ${\bf r}_{k+1}$ で収束判定。収束なら終了。⑤ $\beta_k= ({\bf r}_{k+1}, {\bf r}_{k+1})/({\bf r}_k, {\bf r}_k)$⑥ ${\bf d}_{k+1} = {\bf r}_{k+1} + \beta_k {\bf d}_k $
Pythonによるプログラムの例は以下の記事を参照してください。
前処理付き共役勾配法
収束までの反復回数を減らし計算時間を短縮するために、前処理付き共役勾配法（Pre-conditioned Conjugate Gradient Method）が考えられています。
(1)式に対して、ある行列${\bf P}$の逆行列${\bf P}^{-1}$を左から掛けて、
$${\bf P}^{-1}{\bf A}{\bf x} ={\bf P}^{-1} {\bf b} \tag{2}$$
としても求める結果は変わりません。そこで(1)式の代わりに(2)式を解きます。(2)式を解くメリットは行列の条件数に関係します。
条件数とは方程式を解くときにどれだけ誤差が生じやすいかを示す指標です。例えば(1)式において、行列${\bf A}$やベクトル${\bf b}$を少しだけ変えたときに、解${\bf x}$が大きく変化する場合は条件数が大きく、変化が小さい場合は条件数が小さいと言えます（数学的には条件数は行列の固有値のばらつきを表します）。コンピュータによる数値計算で条件数が大きい問題を解くと、丸め誤差などの小さな誤差によって解が大きく変動し収束しにくくなります。つまり、行列の条件数を小さくすることで安定に計算できるようにするのです。
そこで、ある行列${\bf P}$を用いて、行列${\bf A}$よりも行列${\bf P}^{-1}{\bf A}$の方が条件数が小さくなるようにします。この行列${\bf P}$を前処理行列といいます。大規模な行列を解く場合には、よく前処理付きの解法が使われます。
※前処理行列の求め方は次の記事で紹介しています。
前処理付き共役勾配法のアルゴリズムは以下のようになります。
- 適当な初期値${\bf x}_0$を決める。
- ${\bf r}_0 = {\bf b} – {\bf A} {\bf x}_0$、${\bf d}_0 = {\bf P}^{-1}{\bf r}_0$とする。
- 収束するまで以下の計算を繰り返す。
① $\alpha_k=({\bf r}_k, {\bf P}^{-1}{\bf r}_k)/({\bf d}_k, {\bf A}{\bf d}_k) $
② ${\bf x}_{k+1} = {\bf x}_k + \alpha_k {\bf d}_k$
③ ${\bf r}_{k+1} = {\bf r}_k – \alpha_k{\bf A} {\bf d}_k$
④ 残差 ${\bf r}_{k+1}$ で収束判定。収束なら終了。
⑤ $\beta_k= ({\bf r}_{k+1}, {\bf P}^{-1}{\bf r}_{k+1})/({\bf r}_k, {\bf P}^{-1}{\bf r}_k)$
⑥ ${\bf d}_{k+1} = {\bf P}^{-1}{\bf r}_{k+1} + \beta_k {\bf d}_k $
BiCG法
共役勾配法（CG）は正定値対称行列しか扱うことができません。流体解析などで解かれる行列は非対称な形をしています。そこで、非対称行列に対応した共役勾配法として双共役勾配法（BiCG：Bi-conjugate Gradient）という方法があります。
以下は前処理付きBiCGのアルゴリズムです。
- 適当な初期値${\bf x}_0$を決める。
- ${\bf r}_0 = {\bf b} – {\bf A} {\bf x}_0$、${\bf r}_0^* ={\bf r}_0$、 ${\bf d}_0 = {\bf P}^{-1}{\bf r}_0$、${\bf d}_0^* = {\bf P}^{-1}{\bf r}_0^*$とする。
- 収束するまで以下の計算を繰り返す。
① $\alpha_k=({\bf r}_k^*, {\bf P}^{-1}{\bf r}_k)/({\bf d}_k^*, {\bf A}{\bf d}_k) $
② ${\bf x}_{k+1} = {\bf x}_k + \alpha_k {\bf d}_k$
③ ${\bf r}_{k+1} = {\bf r}_k – \alpha_k{\bf A} {\bf d}_k$
④ ${\bf r}_{k+1}^* = {\bf r}_k^* – \alpha_k{\bf A}^{\rm T} {\bf d}_k^*$
⑤ 残差 ${\bf r}_{k+1}$ で収束判定。収束なら終了。
⑥ $\beta_k= ({\bf r}_{k+1}^*, {\bf P}^{-1}{\bf r}_{k+1})/({\bf r}_k^*, {\bf P}^{-1}{\bf r}_k)$
⑦ ${\bf d}_{k+1} = {\bf P}^{-1}{\bf r}_{k+1} + \beta_k {\bf d}_k $
⑧ ${\bf d}_{k+1}^* = {\bf P}^{-1}{\bf r}_{k+1}^* + \beta_k {\bf d}_k^* $
BiCGSTAB法
BiCGを安定化させた方法として安定化双共役勾配法（BiCGSTAB：Bi-conjugate Gradient Stabilized）という方法もあります。汎用の流体解析ソフトなどでも使われています。当サイトの流体解析ツールCATCFDzeroでも行列解法にBiCGSTABを採用しています。
以下は前処理付きBiCGSTABのアルゴリズムです。
- 適当な初期値${\bf x}_0$を決める。
- ${\bf r}_0 = {\bf b} – {\bf A} {\bf x}_0$、${\bf r}_0^* = {\bf r}_0$、${\bf d}_0 = {\bf r}_0$とする。
- 収束するまで以下の計算を繰り返す。
① $\alpha_k=({\bf r}_0^*, {\bf r}_k)/({\bf r}_0^*, {\bf A}{\bf P}^{-1}{\bf d}_k) $
② ${\bf s}_k = {\bf r}_k – \alpha_k {\bf A}{\bf P}^{-1}{\bf d}_k$
③ $\omega_k = ({\bf A}{\bf P}^{-1} {\bf s}_k, {\bf s}_k) / ({\bf A}{\bf P}^{-1} {\bf s}_k, {\bf A}{\bf P}^{-1} {\bf s}_k)$
④ ${\bf x}_{k+1} = {\bf x}_k + \alpha_k {\bf P}^{-1}{\bf d}_k + \omega_k {\bf P}^{-1}{\bf s}_k$
⑤ ${\bf r}_{k+1} = {\bf s}_k – \omega_k{\bf A} {\bf P}^{-1}{\bf s}_k$
⑥ 残差 ${\bf r}_{k+1}$ で収束判定。収束なら終了。
⑦ $\beta_k= \alpha_k / \omega_k \times ({\bf r}_0^*,{\bf r}_{k+1})/({\bf r}_0^*,{\bf r}_{k})$
⑧ ${\bf d}_{k+1} = {\bf r}_{k+1} + \beta_k ( {\bf d}_k -\omega_k {\bf A}{\bf P}^{-1}{\bf d}_k)$