Preconditioned CG Method

前処理つき共役勾配法(CG法)

Last Update:2009年05月15日

前処理つきCG法(PCG法)

前節ではCG法の収束が行列 $\b{A}$ の固有値分布に強く依存することが分かった。

ここで、解くべき方程式 $\b{A}\b{x}=\b{b}$ を行列 $\b{P}$ を使って次のように変形できたとする。

$\tilde{\b{A}}\tilde{\b{x}}=\tilde{\b{b}}$

但し、 $\tilde{\b{A}}={\b{P}^{-1}\b{A}\b{P}^{-T}$ 、 $\tilde{\b{x}}=\b{P}^T\b{x}$ 、 $\tilde{\b{b}}=\b{P}^{-1}\b{b}$ である。ここで、 $\b{P}$ は特異でない行列である。

$\b{A}$ が正定値対称の場合、 $\tilde{\b{A}}$ も正定値対称となる。そこで連立一次方程式 $\tilde{\b{A}}\tilde{\b{x}}=\tilde{\b{b}}$ にCG法を適応することを考える。

もしも、 $\tilde{\b{A}}$ が $\b{A}$ よりも条件数が小さいようなよい固有値の分布をしているとすると、より早く解を求めることができることが分かる。

ここで行列 $\b{M}=\b{P}\b{P}^T$ とおく。以下、行列 $\b{P}$ 、 $\tilde{\b{A}}$ や $\tilde{\b{b}}$ を用いることなく、 $\b{M}$ の逆行列 $\b{M}^{-1}$ のみを用いて、連立一次方程式 $\tilde{\b{A}}\tilde{\b{x}}=\tilde{\b{b}}$ にCG法を適応しているのと同一になるようにCG法のアルゴリズムを書き直す。

$\b{M}=\b{A}$ の時、 $\tilde{\b{A}}=\b{I}$ となり、変形した連立一次方程式 $\tilde{\b{A}}\tilde{\b{x}}=\tilde{\b{b}}$ は解かずとも解が求められる。

連立一次方程式 $\tilde{\b{A}}\tilde{\b{x}}=\tilde{\b{b}}$ にCG法した際の残差 $\tilde{\b{r}}$ 、探索方向ベクトル $\tilde{\b{p}}$ とする。

$\tilde{\b{r}}=\tilde{\b{b}}-\tilde{\b{A}}\tilde{\b{x}}=\b{P}^{-1}\b{b}-{\b{P}^{-1}\b{A}\b{P}^{-T}\b{P}^T\b{x}=\b{P}^{-1}(\b{b}-\b{A}\b{x})=\b{P}^{-1}\b{r}$

$(\tilde{\b{r}},\tilde{\b{r}})=(\b{P}^{-1}\b{r},\b{P}^{-1}\b{r})=(\b{r},\b{P}^{-T}\b{P}^{-1}\b{r})=(\b{r},\b{M}^{-1}\b{r})$

$(\tilde{\b{p}},\tilde{\b{A}}\tilde{\b{p}})=(\tilde{\b{p}},\b{P}^{-1}\b{A}\b{P}^{-T}\tilde{\b{p}})=(\b{P}^{-T}\tilde{\b{p}},\b{A}\b{P}^{-T}\tilde{\b{p}})$

よって、 $\b{z}=\b{M}^{-1}\b{r}$ 、 $\b{p}=\b{P}^{-T}\tilde{\b{p}}$ とおくと、CG法の係数は次のように表すことができる。

$\alpha_k=\frac{(\tilde{\b{r}}_k,\tilde{\b{r}}_k)}{(\tilde{\b{p}}_k,\tilde{\b{A}}\tilde{\b{p}}_k)}=\frac{(\b{r}_k,\b{z}_k)}{(\b{p}_k,\b{A}\b{p}_k)$

$\beta_k=\frac{(\tilde{\b{r}}_{k+1},\tilde{\b{r}}_{k+1})}{(\tilde{\b{r}}_k,\tilde{\b{r}}_k)}=\frac{(\b{r}_{k+1},\b{z}_{k+1})}{(\b{r}_k,\b{z}_k)}$

CG法の関係式から、

$\Delta\tilde{\b{x}}_k=\tilde{\b{x}}_{k+1}-\tilde{\b{x}}_k=\alpha_k\tilde{\b{p}}_k=\alpha_k\b{P}^T\b{p}_k$

$\Delta\tilde{\b{r}}_k=\tilde{\b{r}}_{k+1}-\tilde{\b{r}}_k=-\tilde{\b{A}}\Delta\tilde{\b{x}}_k=-\alpha_k\tilde{\b{A}}\tilde{\b{p}}_k=-\alpha_k\b{P}^{-1}\b{A}\b{p}_k$

$\tilde{\b{p}}_{k+1}=\tilde{\b{r}}_{k+1}+\beta_k\tilde{\b{p}}_k$

が成り立つ。よってこれから、解 $\b{x}$ ,残差 $\b{r}$ ,ベクトル $\b{p}$ の更新について求めることができる。

$\Delta\b{x}_k=\b{P}^{-T}\Delta\tilde{\b{x}}_k=\alpha_k\b{P}^{-T}\b{P}^T\b{p}_k=\alpha_k\b{p}_k$

$\Delta\b{r}_k=\b{P}\Delta\tilde{\b{r}}_k=-\alpha_k\b{P}\b{P}^{-1}\b{A}\b{p}_k=-\alpha_k\b{A}\b{p}_k$

$\b{p}_{k+1}=\b{P}^{-T}\tilde{\b{p}}_{k+1}=\b{P}^{-T}\tilde{\b{r}}_{k+1}+\beta_k\b{P}^{-T}\tilde{\b{p}}_k=\b{P}^{-T}\b{P}^{-1}\b{r}_{k+1}+\beta_k\b{P}^{-T}\b{P}^T\b{p}_k\\ \qquad \qquad \qquad \quad = \b{M}^{-1}\b{r}_{k+1}+\beta_k\b{p}_k=\b{z}_{k+1}+\beta\b{p}_k$

よって、CG法のアルゴリズムは次のようになる。

アルゴリズム(PCG法)

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0$ , $\b{z}_0=\b{M}^{-1}\b{r}_0$ , $\b{p}_0=\b{z}_0$
2.   For $k=0,1,\ldots,m\quad$ Do:
3.      $\alpha_k=\frac{(\b{r}_k,\b{z}_k)}{(\b{p}_k,\b{A}\b{p}_k)}$
4.      $\b{x}_{k+1}=\b{x}_k+\alpha_k\b{p}_k$
5.      $\b{r}_{k+1}=\b{r}_k-\alpha_k\b{A}\b{p}_k$
6.      $\b{z}_{k+1}=\b{M}^{-1}\b{r}_{k+1}$
7.      $\beta_k=\frac{(\b{r}_{k+1},\b{z}_{k+1})}{(\b{r}_k,\b{z}_k)}$
8.      $\b{p}_{k+1}=\b{z}_{k+1}+\beta\b{p}_k$
9.   $End \quad Do$

********************************

$\b{M}$ は前処理行列と呼ばれる。実際に用いられるのは前処理行列 $\b{M}$ そのものではなく、前処理行列 $\b{M}$ は連立一次方程式 $\b{M}\b{z}=\b{r}$ を解いて $\b{z}$ を求めることさえできればよい。 $\b{M}$ を解くことが $\b{A}$ を解くこととあまり難しさが変わらなければ本末転倒である。つまり前処理行列 $\b{M}$ は $\b{A}$ より簡単に解くことができる必要がある。前処理行列 $\b{M}$ には、 $\b{A}$ と十分に近い $\b{M}\simeq\b{A}$ という性質が求められる。 $\b{M}=\b{A}$ の時、 $\tilde{\b{A}}=\b{I}$ となり、わずか１回の反復で収束する。つまり $\b{M}\simeq\b{A}$ ならば、 $\tilde{\b{A}}\simeq\b{I}$ となりよい固有値分布が得られ収束も早くなると考えられる。またCG法に限って言えば、 $\b{M}=\b{P}\b{P}^T$ と表されていたので、前処理行列 $\b{M}$ は正定値対称でなくてはならない。逆に $\b{M}$ が正定値対称の場合は $\b{M}=\b{P}\b{P}^T$ となるような実行列 $\b{P}$ が存在することが知られている。以上をまとめると、

前処理行列 $\b{M}$ は $\b{M}\b{z}=\b{r}$ を解いて $\b{z}$ を簡単に求めることができる、 $\b{A}$ を近似した正定値対称の行列

ということがいえる。CG法向けの前処理行列が対称である前処理として代表的なものは

対角スケーリング(Diagonal Scaling)
不完全LU分解(ILU分解、Incomplete LU Fractarization)
SSOR法(Symetric Successive Over Relaxation)
点ヤコビ法(Point Jacobi Method)
マルチグリッド法(Multigrid Method)

があげられる。不完全LU分解を前処理として使ったPCG法は特にICCG法と呼ばれよく使われる。

サンプルプログラム

次のサンプルプログラムは次のようなポアソン方程式から生じるような線形一次方程式を

「０レベルの不完全ILU分解を前処理としたCG法」によって解くプログラムである。

DOWNLOAD	ILU_CG.zip

OS:WindowsXP, Windows2000
開発環境:VC++ 2005
言語:C++

このプログラムもまた実行すると収束の履歴がconv_history.datに保存され、gnuplotで見ることができる。

Preconditioned CG Method

目次

前処理つきCG法(PCG法)

アルゴリズム(PCG法)

サンプルプログラム

前処理なしのCG法との比較

参考にしたもの

Books