next up previous
Next: B. 余裕があれば以下のようなことにチャレンジして欲しい Up: 応用数理実験 連立 次方程式に対する CG Previous: 4 MATLAB, Octave でやってみる?

A. CG 法の計算手順について補足

以下、$A$$N$ 次の正値対称行列, $b\in\R^N$ で、ともに与えられている ものとする。

記号の定義
$x_k$ $x^\ast \DefEq A^{-1}b$ の近似と考えるとき、

\begin{displaymath}
r_k\DefEq b-A x_k=A(x^\ast -x_k)
\end{displaymath}

残差 (residual) と呼ぶ。

以下現れるベクトル $p_i$, $r_i$ については、

\begin{displaymath}
\mbox{(\textbf{直交性})}\quad (r_i,r_j) =0\quad\mbox{($i\ne j$)},
\end{displaymath}


\begin{displaymath}
\mbox{(\textbf{共役直交性})}\quad
\langle p_i,p_j\rangle =0\quad\mbox{($i\ne j$)}
\end{displaymath}

が成り立つ。

CG 法にはいくつかのバージョンがあるが、大きく

(I)
2 項漸化式 (Hestenes 版に代表される)
(II)
3 項漸化式 (Rutishauser 版に代表される)
の二つに分類される。この講義では 2 項漸化式版を取り上げる。

CG 法のアルゴリズム

  初期ベクトル ${\bf x_{0}}$ をとる$;\;$  目標とする相対残差 $\varepsilon$ を決める$;$ 

$\Vector{r}_{0} := \Vector{b} - A{\bf x_{0}}; \; {\bf p_{0}}:=\Vector{r}_{0};$
for $k:=0,1,\cdots$ until $\Vert\Vector{r}_{k}\Vert\le \varepsilon\Vert\Vector{b}\Vert \;$ do
    begin
    $\alpha_k := \dsp\frac{(\Vector{r}_{k},{\bf p_{k}})}{({\bf p_{k}},A{\bf p_{k}})};$
${\bf x_{k+1}}:={\bf x_{k}}+\alpha_k {\bf p_{k}};$
$\Vector{r}_{k+1}:=\Vector{r}_{k}-\alpha_k A{\bf p_{k}};$
    $\beta_k:=-\dsp\frac{(\Vector{r}_{k+1},A{\bf p_{k}})}{({\bf p_{k}},A{\bf p_{k}})};$
${\bf p_{k+1}}:=\Vector{r}_{k+1}+\beta_k {\bf p_{k}};$
end

細かい工夫だが、次のアルゴリズムは演算回数が少ない。

CG 法のアルゴリズム (高橋秀俊版)

  初期ベクトル ${\bf x_{0}}$ をとる$;\;$  目標とする相対残差 $\varepsilon$ を決める$;$ 

$\Vector{r}_{0} := \Vector{b} - A{\bf x_{0}};\quad \beta_0 := 1 / (\Vector{r}_{0},\Vector{r}_{0}); \quad {\bf p_{0}}:=\beta_0 \Vector{r}_{0};$
for $k:=0,1,\cdots$ until $\Vert\Vector{r}_{k}\Vert\le \varepsilon\Vert\Vector{b}\Vert \;$ do
    begin
    $\alpha_k := \dsp\frac{1}{({\bf p_{k}},A{\bf p_{k}})};$
${\bf x_{k+1}}:={\bf x_{k}}+\alpha_k {\bf p_{k}};$
$\Vector{r}_{k+1}:=\Vector{r}_{k}-\alpha_k A{\bf p_{k}};$
    $\beta_k:=\dsp\frac{1}{(\Vector{r}_{k+1},\Vector{r}_{k+1})};$
${\bf p_{k+1}}:={\bf p_{k}}+\beta_k \Vector{r}_{k+1};$
end


next up previous
Next: B. 余裕があれば以下のようなことにチャレンジして欲しい Up: 応用数理実験 連立 次方程式に対する CG Previous: 4 MATLAB, Octave でやってみる?
Masashi Katsurada
平成16年11月28日