BiCGSTAB Method

双共役勾配安定化法,BiCGSTAB法

Last Update:2009年05月15日

概要

Bi-CGSTAB法はH.A. van der Vorstによって提案された[1}非対称行列用の反復法である。

BiCG法における回の反復における残差 $\b{r}_j'$ とすると、BiCGSTAB法とはこのBiCG法の残差 $\b{r}_j'$ に以下のような多項式 $\psi_j$ を用いることで残差 $\b{r}_j$ を減少させ、BiCG法を安定化した方法である。

$\b{r}_j=\ch_j(\b{A})\b{r}_j'$

但し、 $\ch_j$ は次のような漸化式で表される多項式である。

$\ch_{j+1}(t)=(1-\omega_j t)\ch_j(t)$

つまり、回目の反復では残差

$||\b{r}_{j+1}||_2=||(\b{I}-\omega_j\b{A})\ch_{j-1}(\b{A})\b{r}_j'||_2$

を最小化するように $\omega_j$ を選ぶ。

基本的な手続き

BiCGSTAB法の手続きで難しいのは、いかにしてBiCG法で用いる残差ベクトルや解ベクトルなどを明示的に持たずにBiCG法に残差の最小化の機能を追加するのかということである。クリロフ部分空間法であるBiCG法では残差などのベクトルは $\b{A}$ の多項式と初期残差を使って表現できる。これを利用してBiCGSTAB法では初期残差からBiCG法における係数を計算する。

BiCGSTAB法における残差ベクトル $\b{r}$ 、解の更新ベクトル $\b{p}$

BiCG法における回反復における残差を $\b{r}_j'$ 、解の更新ベクトルを $\b{p}_j'$ とする。次のように多項式 $\phi_j$ と $\psi_j$ を用いて初期残差 $\b{r}_0$ でこれらを表すとする。つまり、

$\b{r}_j'=\phi_j(\b{A})\b{r}_0 \qquad \qquad \qquad$ $\b{p}_j'=\psi_j(\b{A})\b{r}_0$

この多項式 $\psi$ と $\phi$ は次の漸化式を満たす。

$\phi_{j+1}(t)=\phi_j(t)-\alpha_j t \psi_j(t)$

$\psi_{j+1}(t)=\phi_{j+1}(t)+\beta_j\psi_j(t)$

多項式 $\psi$ と $\chi$ を用いると残差は次のように書ける。

$\b{r}_j=\chi_j(\b{A})\phi_j(\b{A})\b{r}_0$

ここで、 $\chi_j \phi_j$ が出てきたので、これがどのように更新されるのか詳しく調べる。

$\chi_{j+1} \phi_{j+1} = (1-\omega_j t)\chi_j(t)\phi_{j+1} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad = (1-\omega_j t)(\chi_j\phi_j-\alpha_j t \chi_j \psi_j)$

$\chi_j \phi_j$ を更新するためには $\chi_j\psi_j$ が必要であることがわかる。そこで $\chi_j\psi_j$ についても、どのように更新されるのかを調べてみる。

$\chi_{j+1}\psi_{j+1}=\chi_{j+1}(\phi_{j+1}+\beta_j\psi) \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad = \chi_{j+1}\phi_{j+1}+\beta_j(1-\omega_j t)\chi_j\psi_j$

このことから、 $\chi\phi$ と $\chi\psi$ があれば多項式を更新できるということがわかる。そこでBiCGSTAB法におけるベクトル $\b{p}$ を次のように定義する。

$\b{p}_j=\chi_j(\b{A})\psi_j(\b{A})\b{r}_0$

これらを用いて上の漸化式をまとめると、

$\b{r}_{j+1}=(\b{I}-\omega_j\b{A})(\b{r}_j-\alpha_j\b{A}\b{p}_j)$

$\b{p}_{j+1}=\b{r}_{j+1}+\beta_j(\b{I}-\omega_j\b{A})\b{p}_j$

となる。

次はこのベクトル $\b{r}$ 、 $\b{p}$ から、BiCG法の係数 $\alpha$ 、 $\beta$ を作る方法を考える。

$\beta$ の導出

まず始めに $\beta$ を導入する。

$\rho_j$ を次のように定めると、

$\rho_j=(\b{r}_j',\b{r}^*_j')=$\phi_j(\b{A})\b{r}_0,\phi_j(\b{A}^T)\b{r}_0$$

BiCG法より、 $\beta$ は $\rho$ を使って $\beta=\rho_{j+1}/\rho_j$ のように定義された。

しかし、 $\rho_j$ を作る時に必要なBiCG法のベクトル $\b{r}_j'$ 、 $\b{r}^*_j'$ は用いることができない。そこで、次のような値 $\tilde{\rho}_j$ を考える。

$\tilde{\rho}_j=(\b{r}_j,\b{r}^*_0)$

この $\tilde{\rho}_j$ は次のように変形される。

$\tilde{\rho}_j=$\chi_j(\b{A})\phi_j(\b{A})\b{r}_0,\b{r}^*_0$=$ \phi_j(\b{A})\b{r}_0, \chi_j(\b{A}^T)\b{r}^*_0 $ = $\b{r}_j',\chi_j(\b{A}^T)\b{r}^*_0$$

ここで、BiCG法は残差 $\b{r}_j'$ が $\b{r}^*_0$ と $\b{A}^T$ で作られるクリロフ部分空間 $\cal{K}_j(\b{A}^T,\b{r}^*_0)=span\{\b{r}^*_0,\b{A}^T\b{r}^*_0,\ldots,(\b{A}^T)^{j-1}\b{r}^*_0\}$ が直交するという条件の中で解を探すものであった。 $\b{a}$ についての多項式 $\chi_j$ の中で次数がもっとも高い次の項以外の項は $\b{r}_j'$ と内積を計算すると０になることがわかる。よってこの性質を利用すると $\tilde{\rho}$ は $\chi_j$ の最高次数の係数を $\eta_j$ とすると次のように書くことができる。

$\tilde{\rho}_j=$ \b{r}_j', (\b{A}^T)^j\b{r}^*_0 $ \eta_j$

同じように、 $\rho_j$ についてもこの性質を利用して表示してみる。

$\rho_j = (\b{r}_j',\phi_j(\b{A}^T)\b{r}_0) = $ \b{r}_j', (\b{A}^T)^j\b{r}^*_0 $\gamma_j$

但し、 $\gamma_j$ は多項式 $\phi_j$ の最高次数の係数である。上の２式を見比べると $\rho_j$ は次のように作ることができることがわかる。

$\rho_j=\frac{\gamma_j}{\eta_j}\tilde{\rho}_j$

$\eta_j$ と $\gamma_j$ はそれぞれ $\chi_j$ と $\phi_j$ の最高次数の係数であったことを考えると、 $\eta_j$ と $\gamma_j$ は明らかに次の漸化式を満たす

$\eta_{j+1}=-\omega_j\eta_j \qquad \qquad \gamma_{j+1}=-\alpha_j\gamma_j$

よって、 $\beta_j$ は次のようになる。

$\beta_j \quad = \quad \frac{\rho_{j+1}}{\rho_j} \quad = \quad \frac{\eta_j}{\eta_{j+1}} \times \frac{\gamma_{j+1}}{\gamma_j} \times \frac{\tilde{\rho}_{j+1}}{\tilde{\rho}_j} \quad = \quad \frac{\alpha_j}{\omega_j}\frac{\tilde{\rho}_{j+1}}{\tilde{\rho}_j}$

$\alpha_j$ の導出

つぎにBiCG法の係数 $\alpha_j$ について調べる。BiCG法によると $\alpha_j$ は次のように定義された。

$\alpha_j=\frac{(\b{r}_j',\b{r}^*_j')}{(\b{A}\b{p}_j',\b{p}^*_j')}$

分子は前の議論によると、 $(\b{r}_j',\b{r}^*_j')=\rho_j=\frac{\gamma_j}{\eta_j}\tilde{\rho}$ であり計算できる値である。そこで分母の $(\b{A}\b{p}_j',\b{p}^*_j')$ を取り出して詳しく調べる。

$(\b{A}\b{p}_j',\b{p}^*_j')=$\b{A}\b{p}_j',\phi_j(\b{A}^T)\b{r}^*_0 $$

また、BiCG法のアルゴリズムから $\b{A}\b{p}_j'=-\frac{1}{\alpha_j}(\b{r}_{j+1}'-\b{r}_j')$ が成り立つ。

ここでも $\beta$ の導出の際と同様に、BiCG法は残差 $\b{r}_j'$ とクリロフ部分空間 $\cal{K}_j(\b{A}^T,\b{r}^*_0)=span\{\b{r}^*_0,\b{A}^T\b{r}^*_0,\ldots,(\b{A}^T)^{j-1}\b{r}^*_0\}$ を直交させるという条件を用いる。

$\b{r}_j'\bot\cal{K}_j(\b{A}^T,\b{r}^*_0)$

$\b{r}_{j+1}'\bot\cal{K}_{j+1}(\b{A}^T,\b{r}^*_0)\ni\cal{K}_j(\b{A}^T,\b{r}^*_0)$

よって次が成り立つ。

$\b{A}\b{p}_j^*=-\frac{1}{\alpha_j}(\b{r}_{j+1}'-\b{r}_j')\bot\cal{K}_j(\b{A}^T,\b{r}^*_0)$

以上から $\phi_j(\b{A}^T)\b{r}^*_0$ は次数がであるため、最高次の項以外は $\b{A}\b{p}_j^*$ との内積を計算すると０になる。 $\phi_j(\b{A}^T)\b{r}^*_0$ の最高次の係数を $\gamma_j$ とおいていたので、これを用いると

$(\b{A}\b{p}_j',\b{p}^*_j') = $ \b{A}\b{p}_j',(\b{A}^T)^j\b{r}^*_0 $\gamma_j$

同様の議論から次がいえる。

$(\b{A}\b{p}_j,\b{r}^*_0) = $ \b{A}\chi_j(\b{A})\psi_j(\b{A})\b{r}_0,\b{r}^*_0 $ = $ \b{A}\psi_j(\b{A})\b{r}_0,\chi_j(\b{A}^T)\b{r}^*_0 $ = $ \b{A}\b{p}_j',(\b{A}^T)^j\b{r}^*_0 $\eta_j$

上の２式から次のようになる。

$(\b{A}\b{p}_j',\b{p}^*_j') = \frac{\gamma_j}{\eta_j}(\b{A}\b{p}_j,\b{r}^*_0)$

よって $\alpha_j$ を求めることができる。

$\alpha_j=\frac{\frac{\gamma_j}{\eta_j}\tilde{\rho}}{\frac{\gamma_j}{\eta_j}(\b{A}\b{p}_j,\b{r}^*_0)} = \frac{\tilde{\rho}}{(\b{A}\b{p}_j,\b{r}^*_0)}$

$\omega$ の導出

続いて $\omega$ の導出する。これは、残差の２乗ノルムを最小化するように選ばれる。

回目の反復における残差の２乗ノルムは以下のように表された。

$||\b{r}_{j+1}||_2=||(\b{I}-\omega_j\b{A})\ch_{j-1}(\b{A})\b{r}_j'||_2$

簡単のために以下のようにベクトル $\b{s}_j$ を定義する。

$\b{s}_j=\ch_{j-1}(\b{A})\b{r}_j'=\b{r}_j-\alpha_j\b{A}\b{p}_j$

これを用いると、残差の２乗ノルムの２乗は次のようにかける

$||\b{r}_{j+1}||^2_2=||(\b{I}-\omega_j\b{A})\b{s}_j||^2_2=$(\b{I}-\omega_j\b{A})\b{s}_j,(\b{I}-\omega_j\b{A})\b{s}_j$ = (\b{A}\b{s}_j,\b{A}\b{s}_j)\omega_j^2-2(\b{s}_j,\b{A}\b{s}_j)\omega_j+(\b{s}_j,\b{s}_j)$

残差の２乗ノルムが最小であるという条件から

$\frac{\partial}{\partial\omega_j}||\b{r}_{j+1}||^2_2 = 2(\b{A}\b{s}_j,\b{A}\b{s}_j)\omega_j-2(\b{s}_j,\b{A}\b{s}_j)=0$

となる。よって $\omega_j$ は次のようになる。

$\omega_j=\frac{(\b{s}_j,\b{A}\b{s}_j)}{(\b{A}\b{s}_j,\b{A}\b{s}_j)}$

解の更新

残差 $\b{r}$ は次のように更新された

$\b{r}_{j+1} = (\b{I}-\omega_j\b{A})(\b{r}_j-\alpha_j\b{A}\b{p}_j) = \b{s}_j-\omega_j\b{A}\b{s}_j=\b{r}_j-\alpha_j\b{A}\b{p}_j-\omega_j\b{A}\b{s}_j$

つまり、残差の増分は以下のように表される。

$\b{r}_{j+1}-\b{r}_j=-\b{A}(\alpha\b{p}_j+\omega_j\b{s}_j)$

また解の増分と残差の増分の関係から以下がいえる。

$\b{r}_{j+1}-\b{r}_j=-\b{A}(\b{x}_{j+1}-\b{x}_j)$

上の２式を比較すると解の更新の関係が得られる。

$\b{x}_{j+1}=\b{x}_j+\alpha\b{p}_j+\omega_j\b{s}_j$

アルゴリズム（BiCGSTAB法)

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0$ , Choose $\b{r}_0^*$ arbitrary;
2.   $\b{p}_0=\b{r}_0$
3.   For $j=1,2,\ldots,m\quad$ Do:
4.      $\alpha_i=\frac{(\b{r}_j,\b{r}_0^*)}{(\b{A}\b{p}_j,\b{r}_0^*)}$
5.      $\b{s}_j=\b{r}_j-\alpha_j\b{A}\b{p}_j$
6.      $\b{\omega}_j=\frac{(\b{s}_j,\b{A}\b{s}_j)}{(\b{A}\b{s}_j,\b{A}\b{s}_j)}$
7.      $\b{x}_{j+1}=\b{x}_j+\alpha_j\b{p}_j+\omega_j\b{s}_j$
8.      $\b{r}_{j+1}=\b{s}_j-\omega_j\b{A}\b{s}_j$
9.      $\beta_j=\frac{(\b{r}_{j+1},\b{r}^*_0)}{(\b{r}_j,\b{r}^*_0)}\times\frac{\alpha_j}{\omega_j}$
10.     $\b{p}_{j+1}=\b{r}_{j+1}+\beta_j(\b{p}_j-\omega_j\b{A}\b{p}_j)$
11. End Do

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

以上がBICGSTAB法の基本的な手続きである。6,7,8行目が残差最小化のためのステップであり、4,5,9,10行目がBiCG法に相当するステップであると解釈できる。

サンプルプログラム

次のプログラムは次のような移流拡散問題の例

移流拡散問題の例

から生じる行列をBiCGSTAB法で解くプログラムである。

DOWNLOAD	BiCGSTAB.zip

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

結果を表示する場合は次のようにすればよい。

結果の表示

前処理つきBiCGSTAB法

解くべき方程式 $\b{A}\b{x}=\b{b}$ を前処理行列を使って次のように変形する。

$\b{A}\b{M}^{-1}\b{u}=\b{b}, \qquad \qquad \b{x}=\b{M}^{-1}\b{u}$

これは $\b{M}^{-1}$ が $\b{A}$ の右から掛けられているので右前処理と呼ばれる。

係数行列が $\b{A}$ から $\b{A}\b{M}^{-1}$ に変わっている。この $\b{A}\b{M}^{-1}$ の性質が単位行列に近ければ元の方程式よりも早い収束が期待される。

この前処理を使ってBiCGSTAB法で解を求める方法は基本的にBiCGSTAB法で $\b{A}$ の部分を $\b{A}\b{M}^{-1}$ に変えて $\b{u}$ を求めた後に最後に１回解く事で $\b{u}$ を求めるということだが、少し変形すると $\b{u}$ を明示的に持たずに直接 $\b{x}$ の増分を求めることが可能である。

$\hat{\b{p}}=\b{M}^{-1}\b{p}$ 、 $\hat{\b{s}}=\b{M}^{-1}\b{s}$

とおくと、 $\Delta\b{x}=\b{M}^{-1}\Delta\b{u}=\b{M}^{-1}(\alpha\b{p}+\omega\b{s})=\alpha\hat{\b{p}}+\omega\hat{\b{s}}$ というように解の増分を求めることができる。 $\hat{\b{p}}$ 、 $\hat{\b{s}}$ を用いることで $\b{u}$ を明示的に持たなくてもよいことがわかる。また明らかに

$\b{A}\b{M}^{-1}\b{p}=\b{A}\hat{\b{p}}$ 、 $\b{A}\b{M}^{-1}\b{s}=\b{A}\hat{\b{s}}$

である。残差については

$\b{r}=\b{b}-\b{A}\b{x}=\b{b}-\b{A}\b{M}^{-1}\b{u}$

であるから $\b{u}$ に対する残差はそのまま $\b{x}$ に対する残差でもある。よって反復計算しながら $\b{u}$ を求める際に得られる残差の大きさは $\b{x}$ に対する残差の大きさでもある。右前処理は反復しながら解の残差の大きさがわかるので残差の大きさを判断して計算を打ち切るかどうか決める際に便利である。

アルゴリズム（右前処理つきBiCGSTAB法）

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0$ , Choose $\b{r}_0^*$ arbitrary;
2.   $\b{p}_0=\b{r}_0$
3.   For $j=1,2,\ldots,m\quad$ Do:
4.      $\hat{\b{p}}_i=\b{M}^{-1}\b{p}$
5.      $\alpha_i=\frac{(\b{r}_j,\b{r}_0^*)}{(\b{A}\hat{\b{p}}_j,\b{r}_0^*)}$
6.      $\b{s}_j=\b{r}_j-\alpha_j\b{A}\hat{\b{p}}_j$
7.      $\hat{\b{s}}_j=\b{M}^{-1}\b{s}_j$
8.      $\b{\omega}_j=\frac{(\b{s}_j,\b{A}\hat{\b{s}}_j)}{(\b{A}\hat{\b{s}}_j,\b{A}\hat{\b{s}}_j)}$
9.      $\b{x}_{j+1}=\b{x}_j+\alpha_j\hat{\b{p}}_j+\omega_j\hat{\b{s}}_j$
10.     $\b{r}_{j+1}=\b{s}_j-\omega_j\b{A}\hat{\b{s}}_j$
11.     $\beta_j=\frac{(\b{r}_{j+1},\b{r}^*_0)}{(\b{r}_j,\b{r}^*_0)}\times\frac{\alpha_j}{\omega_j}$
12.     $\b{p}_{j+1}=\b{r}_{j+1}+\beta_j(\b{p}_j-\omega_j\b{A}\hat{\b{p}}_j)$
11. End Do

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

以上が右前処理つきBICGSTAB法の基本的な手続きである。

サンプルプログラム

次のプログラムは次のような移流拡散問題の例

移流拡散問題の例

から生じる行列を「ILU分解による右前処理つきBiCGSTAB法」で解くプログラムである。

DOWNLOAD	ILU_BiCGSTAB.zip

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

結果を表示する場合は次のようにすればよい。

結果の表示

n.umetani@gmail.com

BiCGSTAB Method

概要

基本的な手続き

BiCGSTAB法における残差ベクトル $\b{r}$ 、解の更新ベクトル $\b{p}$

$\beta$ の導出

$\alpha_j$ の導出

$\omega$ の導出

解の更新

アルゴリズム（BiCGSTAB法)

サンプルプログラム

前処理つきBiCGSTAB法

アルゴリズム（右前処理つきBiCGSTAB法）

サンプルプログラム

参考にしたもの

Books

Paper