Algebraic Multigrid Methods

AMG法、代数的マルチグリッド法、代数的多重格子法

Last Update:2009年05月15日

概要

AMG法はMultigrid法の一種の方法であり、直交格子だけでなく、非構造格子や一般的な疎行列も解けるように幾何マルチグリッド法を拡張したものである。マルチグリッド法と同じく、多くの問題では反復数がメッシュサイズに依存しないという特徴があり、現在もっとも高速な連立一次方程式の解法の1つといえる。

基本的な手続き

AMGが幾何的マルチグリッド法と違う点は、

行列の成分をもとに粗いグリッドの点を決定すること
行列の成分をもとに粗いグリッドから細かいグリッドへ更新量を射影するProlongation行列などを作ること

である。AMG法では行列を解く準備として、

Coarse Grid/Fine Gridを決める(F/C Splitting)
Prolongation行列を作る
Restriction行列,CoarseGridのオペレータを作る

ProlongagionとRistrictionとCoarseGridオペレータ

Coarse Grid CorrectionはSmootherの収束が遅いような誤差を落とすことが目的である。この誤差を取り除くためにはProlongationオペレータの値域にこの誤差が含まれていなければならない。よってこのSmootherの収束が遅いような誤差がどのような性質を持っているのかを調べてそのような誤差を表現できるようにProlongationのオペレータを作成する。Smootherの収束が遅いような誤差は次のような性質を持っている。

残差が小さい
強い接続上で急激に変化しない

以下この上の性質について詳しく述べる

残差が小さい誤差

古典的反復法における残差は次のように更新される

$\b{e}_{k+1}=(\b{I}-\b{Q}^{-1}\b{A})\b{e}_k$

古典的反復法で収束が遅い場合は $\b{e}_{k+1}\simeq\b{e}_k$ であるから、

$(\b{I}-\b{Q}^{-1}\b{A})\b{e}_k\simeq\b{e}_k$

$\b{e}_k-\b{Q}^{-1}\b{r}_k\simeq\b{e}_k$

$\b{Q}^{-1}\b{r}_k\simeq0$

つまり、 $\b{r}_k\simeq0$

Smootherで収束が遅い誤差は残差が小さいということがいえる。

CoarseGridCorrectionではこのようなSmootherの効率がわるい、残差が小さいという性質をもった誤差を効率よく小さくすることを目標としている。そのため、残差が小さいという性質は非常に重要で、この単純な性質をもった誤差が値域にはいるようにProlongationのオペレータを決め、CoarseGridCorrectionを行う。

このような小さい残差を生むような誤差は代数的に滑らかな誤差(Algebraically Smooth error)という。問題によって代数的に滑らかな誤差でも幾何的には振動している可能性がある。またこのように小さな残差を生むような誤差の別の呼び方として値が０に近いことからNear Kernel errorともいう。

残差が小さい誤差（もう少し正確に）

上の議論は大雑把で、小さい残差とはいっても何と比較して小さいのかわからない。もっと正確には、２つの誤差ノルムに次の関係が成り立つとき残差が小さいという。

$(Ae,D^{-1}Ae)<<(e,Ae)$ つまり、 $(r,D^{-1}r)<<(e,r)$

成分で書くと

$\sum_{i=1}^n\frac{r_i^2}{a_{ii}}<<\sum_{j=1}^ne_jr_j$

上の式は和で表されているが、平均的な性質として

$|r_i|<<a_{ii}|e_j|$

となることが期待される。

強い接続に沿って誤差の変化が小さいという性質

さて、残差が小さいつまり、 $(Ae,D^{-1}Ae)<<(e,Ae)$ であるような誤差はどのような性質を持つのかを調べてみよう。

次の関係が成り立つ。

$(Ae,e)=(D^{-1/2}A^{1/2}A^{1/2}e,D^{1/2}e)=((D^{-1}A)^{1/2}A^{1/2}e,D^{1/2}e)\\ \qquad\qquad \le((D^{-1}A)^{1/2}A^{1/2}e,(D^{-1}A)^{1/2}A^{1/2}e)^{1/2} (D^{1/2}e,D^{1/2}e)^{1/2}\\ \qquad\qquad =(Ae,D^{-1}Ae)^{1/2} (e,De)^{1/2}$

ここではシュワルツの不等式を用いた。これを利用して、 $(Ae,D^{-1}Ae)<<(e,Ae)$ である場合は (e,Ae)<<(e,De) であることがわかる。

成分を用いてこれを書くと次のとおりになる。

$(e,Ae)=\sum_i\sum_ja_{ij}e_ie_j=\sum_i\sum_j\frac{1}{2}\{(-a_{ij})(e_i-e_j)^2+a_{ij}e_i^2+a_{ij}e_j^2\}\\ \qquad\qquad=\sum_i\sum_j\frac{1}{2}(-a_{ij})(e_i-e_j)^2+\sum_i(\sum_ja_{ij})e_i^2\\ \qquad\qquad<<\sum_ia_{ii}e_i^2$

ここで、多くの場合に成り立つ $\sum_ja_{ij}\simeq0$ の関係を用いる。また上式は和について成り立つ式であるが、平均的な性質として各iについても成り立つとすると次のとおりになる。

$\sum_j\frac{1}{2}|a_{ij}|(e_i-e_j)^2 << a_{ii}e_i^2$

つまり以下が成り立つ

$\sum_j\frac{|a_{ij}|}{a_{ii}}\frac{(e_i-e_j)^2}{e_i^2} << 1$

これは $a_{ij}$ が大きければ $e_i\simeq e_j$ が成り立つことを意味している。

つまり、誤差は強い接続上で急激に変化しないという性質を持っていることがわかる。

収束の遅い誤差の例

大きさ１×１の矩形領域における次のような異方性のあるポアソン方程式(Poisson's Equation)

$\{ \begin{array}{ll} \frac{\partial^2\phi}{\partial x^2}+1000\frac{\partial^2\phi}{\partial y^2}=0 \qquad ( \quad in \quad \Omega = \{(x,y):0.0<x,y< 1.0 \} \quad \\ \phi=0 \qquad ( \quad on \quad \partial\Omega \quad ) \end{array}$

において乱数を初期値としてGauss Sidel法を用いたSmootherを300反復させた後の誤差の様子である。

ｘの方向よりもｙの方向の拡散係数が1000倍大きく、領域の内側の全ての点はｙの方向に強い接続、ｘの方向に弱く接続している。ｘ方向には解がランダムに振動しているのに対してｙ方向には滑らかに変化している。Smootherで収束の遅い誤差は強い接続上で急激に変化しないということがわかる。

強い接続上で誤差の変化が少ないのでこの強い接続の方向に点を間引き、補間によって間引かれた点の値を表現することができる。このことを利用してF/C SplittiongやProlongationオペレターの作成を行う。

Strong-Connection と Week-Connection

$\b{A}\b{e}\simeq 0$ であった。このとき行列の成分 $a_{ik}$ が大きい場合、 e_k の値は e_i の値に大きく影響する。逆に $a_{ik}$ が小さい場合、 e_k の値が e_i の値に及ぼす影響は小さくなる。

kがiに大きく影響するとき(つまり、 $a_{ik}>>1$ )はkはiにStrong Connectionであるつまり強く接続しているといい、
kがiにあまり影響していないときはkはiにWeek Connectionであるつまり弱く接続しているという。

F/C SplitingやProlongationオペレター作成時ではこの強い接続、弱い接続をはっきり区別する必要がある。そこで、i行目の非対角成分 $a_{ik}\quad(\quad i\ne k\quad)$ の相対的な大きさから判断して、どの成分kがiに大きく影響するのかを調べ、強い接続、弱い接続に分ける。

$\b{A}$ がＭ行列(対角が正で非対角が負)の場合は一般的にStrong Connectionは次のようにして定義される。

$-a_{ij}>\theta \max_{k\ne i}\{-a_{ik}\}$

$\theta$ は $(\quad 0<\theta\le1\quad)$ であるパラメータで0.25とされることが多い。

Prolongation

マルチグリッド法におけるCoarce Grid CorrectionはSmootherが収束しにくいような誤差を落とすためのものであった。

このような収束しにくい誤差の性質として、

残差が小さい
強い接続間で誤差の変化が小さい

ということが分かった。このような誤差を表現できるようなProlongationのオペレータを作ることが目標である

もちろん粗いグリッドの点の値はそのまま反映されるので、 $i\in C$ のとき e_i=e_i である。以下では $i\in F$ のとき、 e_i どのようにして粗いグリッドの点から補間されるのかを調べよう。

Smootherで収束しにくい誤差は残差が小さかった。

$(\b{A}\b{e})_i\simeq0$

成分ごとに書き出すと次のようになる。

$a_{ii}e_i\simeq-\sum_{j\ne i}a_{ij}e_j$

つまり e_i の値はに接続している点全ての値 e_j から知ることができる。以下ではこの e_i に接続している点の値について述べる。 e_i に接続している点は次の3種類ある。

粗いグリッドの点
に弱く接続している細かいグリッドの点
に強く接続している細かいグリッドの点

細かいグリッドの点に接続している３種類の点

①：粗いグリッドの点 ②：iに弱く接続している細かいグリッドの点 ③：iに強く接続している細かいグリッドの点

に接続している粗いグリッドの点は値 e_j をそのまま使えばいい。

に弱く接続している細かいグリッドの点

$k\in F^w_i$ の場合は補間を行わない。次のような値をとることでCoarseGridCorrectionの安定性を増すことができる。
$\qquad\qquad\qquad e_k=e_i$

に強く接続している細かいグリッドの点

$k\in F^s_i$ の場合は行列の係数の重みを用いた平均として計算する。
$\qquad\qquad\qquad e_k=\frac{\sum_{j\in C_i}a_{kj}e_j}{\sum_{l\in C_i}a_{kl}}$
Smootherで収束しにくい成分は強い接続上で変化が少ないのでこのように強い接続間の補完によって値を近似することができる。このような補間はもしも e_j が定数である場合に e_k が同じ値をとるという性質を満たす。このような定数モードの補間ができることは重要である。
分母が０にならないという条件から、全てのiへ強く影響している点kとiと隣接する細かい節点lは必ず一つは粗い節点に隣接していなければならないということがわかる。通常この性質が満たされるようにF/C Splittingを行う

粗いグリッドから細かいグリッドへの補間

節点の周りの点の値が分かったので節点の値を求めてみよう。

節点iについての残差が小さいということから

$a_{ii}e_i\simeq-\sum_{j\ne i}a_{ij}e_j$

が成り立つのであった。周囲の点を上の3パターンに分けると次のようになる。

$a_{ii}e_i\simeq-\sum_{j\in C_i}a_{ij}e_j-\sum_{k\in F^s_i}a_{ik}e_k-\sum_{k\in F^w_i}a_{ik}e_k$

これに周囲の節点の値を代入して、

$a_{ii}e_i\simeq-\sum_{j\in C_i}a_{ij}e_j-\sum_{k\in F^s_i}a_{ik}\frac{\sum_{j\in C_i}a_{kj}e_j}{\sum_{l\in C_i}a_{kl}}-\sum_{k\in F^w_i}a_{ik}e_i$

これらを変形してまとめると

$(a_{ii}+\sum_{k\in F^w_i}a_{ik})e_i\simeq-\sum_{j\in C_i}\{a_{ij}-\sum_{k\in F^s_i}$\frac{a_{ik}a_{kj}}{\sum_{l\in C_i}a_{kl}$\}e_j$

$e_i\simeq-\sum_{j\in C_i}$\frac{a_{ij}-\sum_{k\in F^s_i}\frac{a_{ik}a_{kj}}{\sum_{l\in C_i}a_{kl}} }{a_{ii}+\sum_{k\in F^w_i}a_{ik}}$e_j$

AMGにおける粗いグリッドから細かいグリッドへの補間式
$w_{ij}=-\frac{a_{ij}-\sum_{k\in F^s_i}\frac{a_{ik}a_{kj}}{\sum_{l\in C_i}a_{kl}} }{a_{ii}+\sum_{k\in F^w_i}a_{ik}}$

$\b{P}=\[W\\I\]$

RistrictionとCoarse Grid Operator

が正定値対称行列の場合はRistrictionの行列は次のようにProlongationの転置(が複素数行列の場合はHarmetian)を用いるとノルムの意味での誤差ノルムの最小化ができる。このような方法は有限要素法との類推からGalerkin法とも呼ばれる

R=P^T

$A_c=\b{R}\b{A}\b{P}$

が正定値対称行列ではない場合にはRistrictionにProlongationの転置を用いる方法は必ずしも最適とはいえない。様々なRistrictionを決定する方法が提案されているが、問題依存であることが多い。

F/C Spliting

AMG法ではコースグリッドを行列の係数を元に作成しなければならない。

細かいグリッド上で粗いグリッドの点の集合を、それ以外の点をとする。

また細かいグリッドにおいて点と接続している点、つまり $a_{ik}\ne0$ の集合を N_i とする。

細かいグリッド上で点に接続している粗いグリッドの点の集合を C_i というふうに書く。つまり $C_i= N_i \cup C$

同様に細かいグリッド上で点に接続している細かいグリッドの点の集合を F_i というふうに書く。つまり $F_i= N_i \cup F$

粗いグリッドの点から細かいグリッドの点に補間ができるためには次の条件を満たす必要がある。

iを細かいグリッドの点とし、jがiに強く接続している点（つまり $j\in S_i$ )とするとjは粗い点 $j\in C$ か、最低１つは粗い点に強く接続( $S_j\cup C\ne \phi$ )していなければならない

計算効率のことを考えるとできるだけ少ない粗いグリッドの点で細かいグリッドへの補間ができたほうがいいと考えられる。

上の条件を満たしつつできるだけ少ないグリッドで計算を行うためには

粗いグリッド同士は強く接続していない

という性質を満たとよい。この粗いグリッド同士の接続がないという条件はグリッドの数を減らして計算量を小さくするという目的のもので、決定的に重要な条件ではない。

収束について

ここではマルチグリッド法の収束について述べる。解くべき行列は正定値対称であるとする。

また、Ristriction I^H_h はProlongation I^h_H の転置であるとする。つまり、 I^H_h=(I^h_H)^T

この場合明らかにコースグリッドのオペレータ A^H=I^H_h A I^h_H も正定値対称である。

ノルム ||.||_0 は次のようなノルムであるとする。 ||e||_0=(e,De) 、但しはの対角成分

ノルム ||.||_1 はオペレーターノルムであるとする。つまり、 ||e||_1=(e,Ae)

また、ノルム ||.||_2 を次のように定義する $||e||_2=(Ae,D^{-1}Ae)$

Prolongationオペレータとコースグリッドコレクションオペレータの値域はオペレータノルムについて直交するという性質について

Prolongationオペレータ I^h_H とコースグリッドコレクションオペレータ T^h の値域がオペレータノルム ||.||_1 について直交する、つまり $R(I^h_H)\bot_A R(T^h)$ ということを示す。

v^H 、 e^h を任意のベクトルとする。

$(I^h_Hv^h,AT^he^h)=(I^h_Hv^h,A\{I-I^h_H(I^H_hAI^h_H)^{-1}I^H_hA\}e^h)=(I^h_Hv^h,Ae^h)-(I^h_Hv^h,AI^h_H(I^H_hAI^h_H)^{-1}I^H_hAe^h)\\ \qquad\qquad\qquad=(I^h_Hv^h,Ae^h)-(v^h,(I^H_hAI^h_H)(I^H_hAI^h_H)^{-1}I^H_hAe^h)=(I^h_Hv^h,Ae^h)-(v^h,I^H_hAe^h)\\ \qquad\qquad\qquad=(I^h_Hv^h,Ae^h)-(I^h_Hv^h,Ae^h)=0$

コースグリッドコレクションはプロロンゲーションオペレータの値域にオペレータノルムにおいて直交する部分空間への射影ということがわかる。

以上より

Te=e-I^h_Hv^H

の両辺のAノルムを計算すると

||e||^2_1=||Te||^2_1+||I^h_Hv^H||^2_1

が成り立つ

δを用いた収束の予測

、

はそれぞれグリッド

におけるスムーザーとコースグリッドコレクションのオペレータであるとする。
各グリッド

について以下の式が成り立つとすると
$||G^me^m||^2_1\le||e^m||^2_1-\delta||T^me^m||^2_1$
ポストスムージングのみのVcycleマルチグリッド法について、誤差のオペレータノルムの収束率が $\sqrt{1-\delta}$ で抑えられる。

以下の式は少し複雑なので、具体的な例を用いて物理的な意味を考えよう。

$||G^me^m||^2_1\le||e^m||^2_1-\delta||T^me^m||^2_1$

$\delta\simeq1$ つまりマルチグリッド法の収束が早いと仮定すると、

コースグリッドコレクションが効率的でない、つまり $||T^me^m||^2_1\simeq||e^m||^2_1$ の場合の誤差については、 $\delta\simeq 1$ より明らかにが成り立つ。つまり、マルチグリッド法が収束している場合は、コースグリッドコレクションが効率的でない誤差をスムーザーが効率的に小さくすることができる。
スムーザーが効率的でない、つまり $||G^me^m||^2_1\simeq||e^m||^2_1$ の場合の誤差については、 $\delta||T^me^m||^2_1<<||e^m||^2_1$ である。ここで $\delta$ が１に近いと仮定したのでが成り立つ。つまり、マルチグリッド法が収束している場合は、スムーザーが効率的でない誤差をコースグリッドコレクションが効率的に小さくすることができる。

以上からδが１に近いほどコースグリッドコレクションとスムーザーの能力を足し合わせた収束が早いということが言える。

2cycle multigrid 法における証明

上式において e^h に T^he^h を代入すると

$||G^hT^he^h||^2_1\le||T^he^h||^2_1-\delta||T^hT^he^h||^2_1$

コースグリッドコレクション T^h は射影オペレータであったから、

T^hT^h=T^h 、 $||T^he_h||_1\le||e^h||_1$ が成り立つ。

これを使うと次のように変形できる

$||G^hT^he^h||^2_1\le||T^he^h||^2_1-\delta||T^hT^he^h||^2_1=(1-\delta)||T^he^h||^2_1\le(1-\delta)||e^h||^2_1$

よって次のようにオペレータノルムでの収束が次のように評価できる。

$||G^hT^he^h||_1\le\sqrt{1-\delta}||e^h||_1$

δの値について

上のδの式をスムーザーとコースグリッドコレクションのオペレータについての２つの分離した式に分けられる。

$||G^me^m||^2_1\le||e^m||^2_1-\alpha||e^m||^2_2$

$||T^me^m||^2_1\le\beta||T^me^m||^2_2$
が成り立つ場合は $\delta=\frac{\alpha}{\beta}$ である。

上の式スムーザーについての式をsmoothing assumptionといい、下のコースグリッドコレクションについての式をapproximation assumptionという。

2grid Multigrid法での簡単な証明

上の式に e^h=T^h e^h を代入すると

$||G^hT^he^h||^2_1\le||T^he^h||^2_1-\alpha||T^he^h||^2_2=||T^he^h||^2_1-\frac{\alpha}{\beta}||T^he^h||^2_1=(1-\frac{\alpha}{\beta})||T^he^h||^2_1\le(1-\frac{\alpha}{\beta})||e^h||^2_1$

$\delta$ を用いた収束の式と見比べると $\delta=\frac{\alpha}{\beta}$ であることがわかる。

$\sqrt{1-\delta}$ がポストスムージングのみのVサイクルマルチグリッド法の収束半径の上限であったので $\alpha$ が大きく、 $\beta$ が小さいほど収束が早い

approximation assumptionについて

2cycle multigrid法では、任意の e^h

について $\min_{v^H}||e^h-I^h_Hv^H||_0^2\le\beta||e^h||^2_1$ が成り立つとき
$||T^he^h||^2_1\le\beta||T^he^h||^2_2$ が成り立つ

証明

$e^h\in R(T^h)$ とすると $R(I_H^h)\bot_AR(T^h)$ から、任意の e^H に対して以下が成り立つ

$||e^h||^2_1=(Ae^h,e^h)=(Ae^h,e^h-I^h_He^H)\\ \qquad\qquad = (D^{1/2}D^{-1/2}A^{1/2}A^{1/2}e^h,e^h-I^h_He^H)\\ \qquad\qquad = (D^{-1/2}A^{1/2}A^{1/2}e^h,D^{1/2}(e^h-I^h_He^H))\\ \qquad\qquad = ((D^{-1}A)^{1/2}A^{1/2}e^h,D^{1/2}(e^h-I^h_He^H))\\ \qquad\qquad \le ||(D^{-1}A)^{1/2}A^{1/2}e^h||\quad ||D^{1/2}(e^h-I^h_He^H)||\\ \qquad\qquad = ||e^h||_2\quad ||e^h-I^h_He^H||_0$

ここでShwarzの不等式、 $D^{-1}A>0$ であることを用いた。

任意の任意の e^h について $\min_{v^H}||e^h-I^h_Hv^H||_0^2\le\beta||e^h||^2_1$ が成り立つとき、

$||e^h||^2_1=||e^h||_2\quad ||e^h-I^h_He^H||_0\le||e^h||_2\sqrt{\beta}||e^h||_1$

両辺2乗して

$||e^h||^2_1\le\beta||e^h||^2_2$

が成り立つ。以上は $e^h\in R(T^h)$ について成り立っていたので e^h:=T^he^h を代入して最終的に

$||T^he^h||^2_1\le\beta||T^he^h||^2_2$

が成り立つ。

参考にしたもの

LINK

An Algebraic Mutigrid Tutorial: http://www.llnl.gov/CASC/linear_solvers/talks/AMG_TUT_PPT/sld001.htm

BOOK

Multigrid Methods

Stephen F.McCormic　著

PAPER

[1] J.W.Ruge and K.Stuben, algebraic multigrid (AMG), in Mutigrid Methods, S. F. McCormic, ed., Frontiers in Appl. Math. 3, SIAM, Philadelphia, 1987, pp.73-130

Made by Nobuyuki UMETANI 梅谷　信行

n.umetani@gmail.com