Convergence of CG Method

共役勾配法(CG法)の収束性

Last Update:2009年05月15日

チェビシェフ多項式による収束性の予測

ここで、 $||\b{x}||_{\b{A}}$ は次のように定義される作用素ノルムである。

$||\b{x}||^2_{\b{A}}=(\b{x},\b{A}\b{x})$

CG法はKrylov部分空間 $\cal{K}_m$ の中で誤差の $\b{A}$ ノルムを最小化する手法であった。

つまり、方程式の解を $\b{x}^*$ 、回反復後の近似解を $\b{x}_m$ とすると、 $\b{x}_m-\b{x}_0$ がKrylov部分空間に含まれることから、

$\b{x}_m=\b{x}_0+q_m(\b{A})r_0$

とかける、但し、多項式 q_m は m-1 次である。

回反復目のCG法の誤差 $\b{e}_m$ の $\b{A}$ ノルム $||\b{e}_m||_{\b{A}}$ は、

$||\b{e}_m||_{\b{A}=||\b{x}^*-\b{x}_m||_{\b{A}}=||\b{x}^*-(\b{x}_0+q_m(\b{A})r_0)||_{\b{A}}=||\b{e}_0-q_m(\b{A})\b{A}\b{e}_0||_{\b{A}}=||(\b{I}-\b{A}q_m(\b{A}))\b{e}_0||_{\b{A}}$

のように初期誤差 $\b{e}_0$ を用いて表すことができる。ここで、 $\b{e}_0=\b{A}\b{e}_0$ を用いた。

q_m はこれを m-1 次の多項式全体 $P_{m-1}$ の中で最小化するので

$||\b{e}_m||_{\b{A}}=||(\b{I}-\b{A}q_m(\b{A}))\b{e}_0||_{\b{A}}=\quad\min_{\small{q\in P_{m-1}}}||(\b{I}-\b{A}q(\b{A}))\b{e}_0||_{\b{A}}$

となる。

ここで簡単のため $\b{I}-\b{A}q(\b{A})=r(\b{A})$ とおく。但し、は r(0)=1 を満たす次の多項式である。これを用いると上の式は

$||\b{e}_m||_{\b{A}}=\quad\min_{\small{r\in P_m \quad r(0)=1}}||r(\b{A})\b{e}_0||_{\b{A}}$

と表すことができる。

固有値による誤差のノルムの評価

ここで、 $\lambda_1,\ldots,\lambda_n$ は $\b{A}$ の固有値であるとし、 $\b{x}_1,\ldots,\b{x}_n$ は固有ベクトルであるとする。 $\xi_1,\ldots,\xi_n$ は $\b{e}_0$ を固有ベクトルの線形結合で表した時の係数であるとする

$||r(\b{A})\b{e}_0||_{\b{A}} = ||\sum_i^n r(\lambda_i)\xi_i\b{x}_i||_{\b{A}} \le |\sum_i^n r(\lambda_i)| \quad ||\sum_i^n\xi_i\b{x}_i||_{\b{A}} \le \max_i \quad |r(\lambda_i)| \quad ||\b{e}_0||_{\b{A}} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \le \max_{\small{\lambda\in[\lambda_{min},\lambda_{max}]}} \quad |r(\lambda)| \quad ||\b{e}_0||_{\b{A}}$

これを用いると回反復時における誤差の $\b{A}$ ノルム $||\b{e}_m||_{\b{A}}$ は次のように評価される。

$||\b{e}_m||_{\b{A}} \le \min_{\small{q\in P_m \quad r(0)=1}} \quad \max_{\small{\lambda\in[\lambda_{min},\lambda_{max}]}} \quad |r(\lambda)| \quad ||\b{e}_0||_{\b{A}}$

つまり、 r(0)=1 を満たす次の多項式の中で、区間 $\lambda\in[\lambda_{min},\lambda_{max}]$ での最大の絶対値を最小値によって、誤差の $\b{A}$ ノルムは抑えられる。

誤差のノルムを抑える係数 $\rho$ （イメージ）

チェビシェフ多項式

チェビシェフ多項式はある区間で絶対値の最大値を最小化するという性質をもった多項式である。

具体的には次のチェビシェフ多項式 C_n(x) は最高次の係数が１であるような任意の次の多項式 p_n(x) の中で

$\max_{\small{|x|<1}} |p_n(x)|$

が最小になるような多項式を $2^{n-1}$ 倍したものと一致する。

チェビシェフ多項式は次のように定義される

$C_k(t)=\left(\begin{array} \cos(k\cos^{-1}t) \qquad \qquad |t|\le1 \\ \cosh(k\cosh^{-1}t) \quad |t|>1 \end{array}$

一見複雑な超越関数に見えるが実はこれは単純な多項式である。

チェビシェフ多項式同士には次のような関係がある。

$C_{k+1}(t)=2tC_k(t)-C_{k-1}(t)$

これを使えば C_1(t)=t 、 C_2(t)=2t^2-1 、 C_3(t)=4t^3-3t 、 C_4(t)=8t^4-8t^2+1 、 $\cdots$ となることがわかる。

チェビシェフ多項式は次のようにも表すことができる。

$C_k(t)=\frac{1}{2}\left[(t+\sqrt{t^2-1})^k+(t+\sqrt{t^2-1})^{-k}\right]$

チェビシェフ多項式の絶対値は区間[-1,1]で最大値１をとることが知られている。

チェビシェフ多項式の変数変換とスケーリング

次のチェビシェフ多項式 C_k は区間[-1,1]の間の絶対値の最大値を最小化する最高次数の係数が $2^{k-1}$ の次の多項式であった。

チェビシェフ多項式の最適化区間を [-1,1] から $[\alpha,\beta]$ に変換するために、 $t=1+2\frac{t'-\beta}{\beta-\alpha}$ という変数変換を用いたものを $\bar{C}_k$ とすると

$\bar{C}_k(t)=C_k(1+2\frac{t-\beta}{\beta-\alpha})$

となる。これにさらに適当なスケーリングを行って、引数が $\gamma$ の時に値１を取るようにした関数 $\hat{C}_k$ とすると

$\hat{C}_k(t)=\frac{\bar{C}_k(t)}{\bar{C}_k(\gamma)}=\frac{C_k(1+2\frac{t-\beta}{\beta-\alpha})}{C_k(1+2\frac{\gamma-\beta}{\beta-\alpha})}$

$\hat{C}_k$ は区間 $[\alpha,\beta]$ の間の絶対値の最大値を最小化する引数 $\gamma$ において１をとる次の多項式である。

$C_k(1+2\frac{t-\beta}{\beta-\alpha})$ の絶対値は区間 $[\alpha,\beta]$ で最大値１をとることから

$\hat{C}_k(t)$ の絶対値は区間 $[\alpha,\beta]$ で最大値 $\frac{1}{|C_k(1+2\frac{\gamma-\beta}{\beta-\alpha})|}$ をとる。つまり、

$\min_{\small{p\in P_k,p(\gamma)=1}} \quad \max_{\small{t\in[\alpha,\beta]}} \quad |p(t)| = \frac{1}{|C_k(1+2\frac{\gamma-\beta}{\beta-\alpha})|}$

CG法の収束性の評価

引数０のとき１をとる次の多項式 $r(\lambda)$ の中で、区間 $[\lambda_{min},\lambda_{max}\]$ において絶対値の最大値の最小値は、上の式の $\alpha=\lambda_{max}$ 、 $\beta=\lambda_{min}$ 、 $\gamma=0$ とおくことで

$\min_{\small{r\in P_k,r(0)=1}} \quad \max_{\small{\lambda\in[\lambda_{min},\lambda_{max}]}} \quad |r(\lambda)| = \frac{1}{|C_m(1+2\frac{-\lambda_{max}}{\lambda_{max}-\lambda_{min}})|} = \frac{1}{|C_m(1+2\frac{\lambda_{min}}{\lambda_{max}-\lambda_{min}})|} = \frac{1}{C_m(1+2\frac{\lambda_{min}}{\lambda_{max}-\lambda_{min}})} = \frac{1}{C_m(1+2\eta)}$

ここで、 $\eta=\frac{\lambda_{min}}{\lambda_{max}-\lambda_{min}}$ とおいた。また、 $C_m(t)>1 \quad ( \quad for \quad t > 1 \quad )$ を用いた。

ここから、 $C_m(1+2\eta)$ を評価する

$C_m(t)=\frac{1}{2}\left[(t+\sqrt{t^2-1})^m+(t+\sqrt{t^2-1})^{-m}\right] \ge \frac{1}{2}(t+\sqrt{t^2-1})^m$ を用いる。

$C_m(1+2\eta)\ge\frac{1}{2}(1+2\eta+\sqrt{(1+2\eta)^2-1})^m = \frac{1}{2}(1+2\eta+2\sqrt{\eta(\eta+1)})^m = \frac{1}{2}\{(\sqrt{\eta}+\sqrt{1+\eta})^2\}^m$

$(\sqrt{\eta}+\sqrt{1+\eta})^2 = \frac{(\sqrt{\lambda_{min}}+\sqrt{\lambda_{max}})^2}{\lambda_{max}-\lambda_{min}} = \frac{\sqrt{\lambda_{min}}+\sqrt{\lambda_{max}}}{\sqrt{\lambda_{max}}-\sqrt{\lambda_{min}}} = \frac{\sqrt{\kappa}+1}{\sqrt{\kappa}-1}$ を用いる。

但し、 $\kappa=\frac{\lambda_{max}}{\lambda_{min}}$ とおいた。 $\kappa$ は正定値対称行列のスペクトル条件数(Spectorum Condition Number)である。

$C_m(1+2\eta) \ge \frac{1}{2}(\frac{\sqrt{\kappa}+1}{\sqrt{\kappa}-1})^m$

$\min_{\small{r\in P_k,r(0)=1}} \quad \max_{\small{\lambda\in[\lambda_{min},\lambda_{max}]}} \quad |r(\lambda)| \le 2(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1})^m$

これを用いると

$||\b{e}_m||_{\b{A}} \le 2(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1})^m||\b{e}_0||_{\b{A}}$

となり、誤差の $\b{A}$ ノルムに対して上のように評価できることがわかる。

このことから条件数が小さいほど収束が早いということがわかる。

反復数と条件数

上の式から、予測される収束率 $\rho_e$ は次のようにかける。

$\rho_e=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$

これは $\kappa$ が十分大きい場合、次のように変形される。

$\rho_e=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}=\frac{(\sqrt{\kappa}+1)-2}{\sqrt{\kappa}+1}=1-\frac{2}{\sqrt{\kappa}+1}\simeq 1-\frac{2}{\sqrt{\kappa}}$

例えば $\kappa$ が a^2 倍されたとしたときに収束率 $\rho_e'$ がもとの収束率 $\rho_e$ からどのように近似的に求められるかを考える。

$\rho_e'=1-\frac{2}{\sqrt{a^2\kappa}}=1-\frac{1}{a}\frac{2}{\sqrt{\kappa}}\simeq(1-\frac{2}{\sqrt{\kappa}})^{\frac{1}{a}}={\rho_e}^{\frac{1}{a}}$

ここで $\kappa$ が十分大きいという近似を用いた。

$\kappa$ が a^2 倍されると収束率がだいたい $\frac{1}{a}$ 乗されることがわかる。これは同じ収束を得るまでに倍の反復数を要するということを意味する。

メッシュサイズとの関係

ラプラス演算子を有限要素法離散化した場合、メッシュサイズと条件数 $\kappa$ はだいたい次のような関係にあることが知られている。

$\kappa\propt\frac{1}{h^2}$

このことから例えばメッシュサイズが倍されるとだいたい次のことがいえる。

hが倍　 $\rightarrow$ 条件数 $\kappa$ が 1/a^2 倍　 $\rightarrow$ 収束率が乗　 $\rightarrow$ 反復数が 1/a 倍

つまり、メッシュサイズが半分になると反復数がだいたい２倍になるということがわかる。

サンプルプログラムでの確認

サンプルプログラムではメッシュサイズを容易に変化させることができる。

main関数の中の最初の宣言で問題を解くための正方形の領域での四角形メッシュを定義している。

field::CElemAry_SquarePow2 ea(6,1.0,0.5,0.5);

第一の初期化変数はデフォルトで６になっているこれは一辺を 2^6 つに分割することを意味している。

第二の初期化変数は1.0となっているがこれは一辺の長さが1.0であることを意味している。

つまりデフォルトでメッシュサイズは 1/2^6 である。第一初期化変数を変化させてメッシュサイズと反復数が上のような反比例の関係に従うかどうか確認してみよう。

メッシュサイズ
反復数	47	94	188	376	758

このように単純な問題においては正確に上の性質を満たしていることがわかる。

収束の超線形性

CG法の収束は超線形性を示す。収束が超線形であるとは、対数グラフで残差の変化を表したときにグラフが直線的ではなく徐々に加速して始めより急な収束を示すことをいう。

次のグラフは前のサンプルプログラムの収束の様子である。

青い線は始めの収束の様子を表している。緑の線のようにあきらかに収束が加速しているのがわかる。このような超線形の収束はJacobi法やGauss-Sidel法や古典的などの古典的な反復法ではみられず、Krylov部分空間法ならではのものである。

上で行ったチェビシェフ多項式による収束の予測では線形の収束しか予測できないので、実際のところ上の式で得られる収束よりも現実の方がずっと早い収束が得られることが多い。

参考にしたもの

Books

Made by Nobuyuki UMETANI 梅谷　信行

n.umetani@gmail.com