数値計算

第3回 FEMを勉強してみる

今回は 中島 研吾 先生の講演映像の⑤に対応する内容として,連立方程式の数値的解法の話をする.これによって,全体方程式が大規模になってもコンピュータが自動で計算してくれることになる.ICCG法をゴールとして解説するが, 今回に関しても最終的なアルゴリズムだけ知っておけば実装は何とかなる.また,連立方程式の計算はFEMからほぼほぼ独立した問題といってよい.なので,連立方程式をコンピュータが解けることがイメージできるのであれば,先に進んでもらって構わない.最後まで解ける保証がないのに大規模なFEMを勉強するのは気持ち悪いって人向けにこの順番にした(僕がそうだった).

また,今回に関しては全体を通して 藤澤 誠 先生のサイトを参考にしている.

勾配法


  • 係数行列をA, 解ベクトルを x,定数ベクトルを b とすると,一次連立方程式は次のように表現できる:
$$A\boldsymbol{x}=\boldsymbol{b}$$

  • 連立方程式においても微分方程式の近似解と同様に,近似解がどれだけ方程式を成り立たせるかの指標として残差 r (x ) が使えそう:

$$\boldsymbol{r}(\boldsymbol{x})=A\boldsymbol{x}-\boldsymbol{b}$$

  • r (x ) はベクトルを引数とするベクトルである
  • 近似解の信頼度がベクトルで表されるのは都合が悪い
  • そこで近似解ベクトルを引数として0以上のスカラ(値1つだけ)を求める誤差関数 f (x ) を考える 
  • 一番素朴なアイデアとして,r (x ) の大きさが f (x ) として使えそう
  • 誤差関数 f (x ) が小さくなるよう,近似解を x を少しずつずらして,真の解に近づけていく方法を勾配法という
  • 少しづつずらすのは,計算されるずらし方の方向性が場所によって異なり,最初に求めた方向が必ずしも信用できるものではないからである(道案内が「ここを真っすぐ」で済むとは限らないのと一緒.解の世界にも障害物があって,迂回が必要になることはある.)
  • 降下法は一般に,f の定義,方向の決め方,距離の決め方,x 初期値,打切り誤差(終了条件) が必要となり,それによって結果が異なる
  • 方向の決め方については f (x ) の勾配(一番急な方向,傾きの多次元版)の逆方向に進む,最急降下法が代表的である:

$$\boldsymbol{x}_{\bf 次}=\boldsymbol{x}_{\bf 現在}-({\bf移動距離})\times\nabla f(\boldsymbol{x}_{\bf 現在})$$

$$=\boldsymbol{x}_{\bf 現在}-({\bf移動距離})\times\left(\frac {\partial f(\boldsymbol{x}_{\bf 現在}) }{\partial x_{\bf 1,現在}}, \frac {\partial f(\boldsymbol{x}_{\bf 現在}) }{\partial x_{\bf 2,現在}}, … \right)$$

  • 誤差関数,移動距離,初期値の設定方法については様々であり,それによってさらに細分化される(今回は扱わない)

共役勾配法(CG法)


  • FEMで扱うの A は特殊な行列(正定値対称行列)である
  • 正定値対称行列に特化した勾配法共役勾配法(CG法)である
  • CG法は,f の定義,方向の決め方,距離の決め方 をカバーしている
  • 設定しなくてはならないのは初期値と打切り誤差であるが,初期値はすごく適当に決めてもうまくいく
  • とりあえずアルゴリズムだけ先に示す:
  • そのままではよく分からない(アルゴリズムは洗練された手順なので意図が分かりにくくなりがち
  • 何をしているのか,順を追って説明する

① 誤差関数 f (x )

  • まず,共役勾配法では,近似解をx,厳密解をy として誤差を以下のようなものと考える:

$$(\boldsymbol{x}-\boldsymbol{y})\cdot A (\boldsymbol{x}-\boldsymbol{y})$$

  • (補足)「•」 はベクトルの内積を表す.ベクトルを一列の行列と考えるならば,内積 xAx を (x^T)Ax と表すことができるのであるが,今回はベクトルと行列でそれぞれ扱っている対象が全然違うので,ベクトルはなるべくベクトルであることを明確にしたい.行列を太字で表現しないのもそういう意図による.
  • これは「零ベクトル以外の任意のベクトル x に対して xAx > 0 が成り立つような行列A を正定値行列とよぶという定義を利用したものであり,逆にこれが小さくなれば x は零ベクトルに近づくだろうというアイデア
  • 誤差の式を展開する:

$$(\boldsymbol{x}-\boldsymbol{y})\cdot A (\boldsymbol{x}-\boldsymbol{y})=\boldsymbol{x} \cdot A (\boldsymbol{x}-\boldsymbol{y}) \,- \boldsymbol{y} \cdot A (\boldsymbol{x}-\boldsymbol{y}) $$

$$= \boldsymbol{x} \cdot A \boldsymbol{x} \,-\, \boldsymbol{x} \cdot A \boldsymbol{y} \,-\, \boldsymbol{y} \cdot A \boldsymbol{x} + \boldsymbol{y} \cdot A \boldsymbol{y}$$

  • A は対称行列なので y・Ax = x・Ay が成り立つ(試せばすぐ分かる)から

$$ = \boldsymbol{x} \cdot A \boldsymbol{x} \,-\, \boldsymbol{x} \cdot A \boldsymbol{y} \,-\, \color{red}{\boldsymbol{x} \cdot A \boldsymbol{y}} + \boldsymbol{y} \cdot A \boldsymbol{y}$$

$$=\boldsymbol{x} \cdot A \boldsymbol{x}  \,-\, 2\boldsymbol{x} \cdot A \boldsymbol{y} + \boldsymbol{y} \cdot A \boldsymbol{y}$$

  • さらに,y は厳密解なので Ay = bであるから:

$$ =\boldsymbol{x} \cdot A \boldsymbol{x}  \,-\, 2\boldsymbol{x} \cdot \color{red}{\boldsymbol{b}} + \boldsymbol{y} \cdot \color{red}{\boldsymbol{b}}$$

  • x をあれこれと変えてこれを最小化するとき,yb は定数とみなせる
  • よって,(x y )A (x y ) を最小化するx は,以下の関数 f (x ) を最小化する x と同じである:

$$ \color{red}{f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x} \cdot A \boldsymbol{x}  \,-\, \boldsymbol{x} \cdot \boldsymbol{b}}$$

  • 半分にしたのは後々に式が見やすくなるようにするためで,数学的に大した意味はない

② 移動距離 α

  • 移動方向 p が既知であるとき,最適な移動距離を考える
  • 移動後の誤差関数の値はf (x + αp ) であるから,これが最小となるような α を探せばよい
  • とりあえず f が極小(傾きが0)となる α を探せば良しとする:

$$ \frac{\partial f(\boldsymbol{x} + \alpha \boldsymbol{p})}{\partial \alpha}=0$$

  • f (x + αp ) は先に求めた の式の x x + αp を代入するだけ
  • あとはこれをで α 微分して,それをイコールゼロとして解くだけ
  • 大した計算ではないので結果だけ書いておく:

$$ \color{red}{ \alpha=\frac{\boldsymbol{p} \cdot \boldsymbol{r}}{\boldsymbol{p} \cdot  A \boldsymbol{p}}}$$

$$ただし, \color{red}{ \boldsymbol{r} = \boldsymbol{b}-A \boldsymbol{x}}$$

  • r は残差であり,CG法においては移動距離と方向のうまい仲介役になってくれる
  • 残差も近似解の更新に伴って更新する必要がある
  • ただし,更新に伴う変化量は以下のように計算されるので,前の残差にこの量を足すだけで済む:

$$ \boldsymbol{r}_{\bf 次} \,-\, \boldsymbol{r}_{\bf 現在}= \boldsymbol{r}(\boldsymbol{x} + \alpha \boldsymbol{p}) \,-\, \boldsymbol{r}(\boldsymbol{x})$$

$$=(\boldsymbol{b}-A(\boldsymbol{x} + \alpha \boldsymbol{p}) )-(\boldsymbol{b}-A \boldsymbol{x})=-\alpha A \boldsymbol{p}$$

$$ ∴ \boldsymbol{r}_{\bf 次} = \boldsymbol{r}_{\bf 現在} -\alpha A \boldsymbol{p}$$

③ 移動方向 p

  • 移動方向は共役勾配という特殊な勾配に基づいて計算され,不思議なことに普通の勾配よりもうまく方向が求められる
  • 共役勾配と勾配の違いについては Keiさんのブログ ではアニメーション付きで説明されているので直感的に理解できる
  • こうした背景があるのであるが,極力この辺は触れずに説明する(ちょっと無理があるかも)
  • 移動方向は残差をベースに計算される
  • 実際,p の初期値は残差の初期値:r (x_0 ) = b A x_0 に設定する
  • 次ステップの p更新済みの残差に,現在の p に重み β をかけたものを加えて計算する:

$$ \boldsymbol{p}_{\bf 次} = \boldsymbol{r}_\color{red}{\bf 次} + \beta \boldsymbol{p}$$

  • β がかかるのは残差でなく,更新前の移動方向なことに注意(α と関係性が逆)
  • 残差が更新済みであることにも注意
  • 後は β の決め方が問題になる
  • 唐突ながら,以下のことが成り立つ:

$$ A\boldsymbol{p} \cdot \left( \boldsymbol{y}\, -\, \boldsymbol{x}_{\bf 次}\right) = 0$$

クリックで証明表示
  • x_次 = x + αp を代入:

$$ A\boldsymbol{p} \cdot \left( \boldsymbol{y}\, -\, \boldsymbol{x}_{\bf 次}\right) =A\boldsymbol{p} \cdot \left\{ \boldsymbol{y} \,-\, (\boldsymbol{x} + \alpha \boldsymbol{p})\right\}$$

  • 内積の前後は交換可能:

$$=\left\{ \boldsymbol{y} \,-\, (\boldsymbol{x} + \alpha \boldsymbol{p})\right\} \cdot A\boldsymbol{p}$$

  • ベクトルを行列として考えたら結合則(積をとる順番はどちらが先でもよい)が成り立つので:

$$=\left\{ \boldsymbol{y} \,-\, (\boldsymbol{x} + \alpha \boldsymbol{p})\right\} A\cdot \boldsymbol{p}$$

  • また,A が対称行列のとき (xA)^T = (A^T)(x ^T) = Ax ^T であり,現在はベクトルの転置は気にしないことにしている(列ベクトルも行ベクトルも同じとして扱っている)ので,xA = Ax より:

$$=A\left\{ \boldsymbol{y} \,-\, (\boldsymbol{x} + \alpha \boldsymbol{p})\right\}\cdot \boldsymbol{p}$$

  • Aを分配する:

$$=\left\{A \boldsymbol{y} \,-\, A(\boldsymbol{x} + \alpha \boldsymbol{p})\right\}\cdot \boldsymbol{p}$$

  • Ay = b より:

$$=\left\{ \boldsymbol{b} \,-\, A(\boldsymbol{x} + \alpha \boldsymbol{p})\right\}\cdot \boldsymbol{p}$$

$$=\left( \boldsymbol{b} \,-\, A\boldsymbol{x} \,-\, \alpha A  \boldsymbol{p}\right)\cdot \boldsymbol{p}$$

  • b – Ax = r より:

$$=\left( \boldsymbol{r}  \,-\, \alpha A \boldsymbol{p}\right)\cdot \boldsymbol{p}$$

$$=\boldsymbol{r} \cdot \boldsymbol{p}  \,-\, \alpha A \boldsymbol{p} \cdot \boldsymbol{p}$$

$$=\boldsymbol{r} \cdot \boldsymbol{p}  \,-\, \alpha  \boldsymbol{p} \cdot A \boldsymbol{p}$$

  • α の式を代入すると:

$$=\boldsymbol{r} \cdot \boldsymbol{p}  \,-\, \frac{\boldsymbol{p} \cdot \boldsymbol{r}}{\boldsymbol{p} \cdot  A \boldsymbol{p}} \, \boldsymbol{p} \cdot A \boldsymbol{p}$$

$$=\boldsymbol{r} \cdot \boldsymbol{p}  \,-\, \boldsymbol{p} \cdot \boldsymbol{r} = 0$$

  • ここで,yx_次 は,次の更新でx が厳密解にたどり着けるような移動 (α_次)(p_次) である
  • ただし,y は当然ながらわからないので,これらを直接等号で繋いでも p_次 は求まらない
  • そこで,理想的な移動 yx_次 が満たしている上の式を,(α_次)(p_次) も満たしているとする

$$ A\boldsymbol{p} \cdot \alpha_{\bf 次} \boldsymbol{p}_{\bf 次} = 0$$

$$\Longleftrightarrow \; A\boldsymbol{p} \cdot \boldsymbol{p}_{\bf 次} = 0$$

$$\Longleftrightarrow \; A\boldsymbol{p} \cdot (\boldsymbol{r}_{\bf 次} + \beta \boldsymbol{p}) = 0$$

$$\Longleftrightarrow \;  (\boldsymbol{r}_{\bf 次} + \beta \boldsymbol{p}) \cdot A\boldsymbol{p}= 0$$

$$\Longleftrightarrow \;  \boldsymbol{r}_{\bf 次} \cdot A\boldsymbol{p} + \beta (\boldsymbol{p} \cdot A\boldsymbol{p})= 0$$

$$\Longleftrightarrow \;  \color{red}{\beta = -\frac{\boldsymbol{r}_{\bf 次} \cdot A\boldsymbol{p}}{\boldsymbol{p} \cdot A\boldsymbol{p}}}$$

  • これによって β が求まったので移動方向が決定できる
  • また,2つ目の式について順番を入れ替えると:
  • 以上より,移動方向,移動距離の更新式ができたので,解を繰り返し更新することができる
  • 終了条件はf で判断した方がよさそうであるが,実際には残差(移動距離と考えてもよい)のノルムで判定する
  • これはf がどれほどの値に落ち着くかはまちまちであり,成果(解の精度)を要求するよりも,アルゴリズムにどれぐらいがんばってもらうかを決めておく方が,収束しやすいからである
  • これでアルゴリズムは成立するが,もう少し工夫できるところがある
  • β の導出における2つ目の式について順番を入れ替えると:

$$\boldsymbol{p}_{\bf 次} \cdot A\boldsymbol{p} = 0$$

  • これは p の更新前と更新後が A に関して共役であることを意味する
  • 数学的帰納法を用いることで,このアルゴリズムにおける各ステップで求まる p は互いに共役であることが導かれる(ただし同一の p の場合は除く)
  • さらに以下のことが次々と導かれる(証明略):

$$\boldsymbol{p} \cdot \boldsymbol{r}_{\bf 次} = 0$$

$$\boldsymbol{p}_i \cdot \boldsymbol{r}_{k} = 0 \;(i < k)$$

$$\boldsymbol{r}_i \cdot \boldsymbol{r}_{j} = 0 \;(i \neq j)$$

$$\boldsymbol{p} \cdot \boldsymbol{r} = \boldsymbol{r} \cdot \boldsymbol{r}$$

  • これによってα,β の式が簡単になる:

$$\alpha= \frac{\boldsymbol{p} \cdot \boldsymbol{r} } {\boldsymbol{p} \cdot A\boldsymbol{p}}= \color{red}{\frac{\boldsymbol{r} \cdot \boldsymbol{r} } {\boldsymbol{p} \cdot A\boldsymbol{p}}}$$

$$\beta = -\frac{\boldsymbol{r}_{\bf 次} \cdot A\boldsymbol{p}}{\boldsymbol{p} \cdot A\boldsymbol{p}}=-\frac{\boldsymbol{r}_{\bf 次} }{\boldsymbol{p} \cdot A\boldsymbol{p}} \cdot \frac{\boldsymbol{r}_{\bf 次} – \boldsymbol{r}_{\bf 現在} }{-\alpha}$$

$$=\frac{\boldsymbol{r}_{\bf 次} \cdot \boldsymbol{r}_{\bf 次} \,-\, \boldsymbol{r}_{\bf 次} \cdot \boldsymbol{r}_{\bf 現在} }{\boldsymbol{p} \cdot A\boldsymbol{p}} \frac{1}{\alpha}$$

$$=\frac{\boldsymbol{r}_{\bf 次} \cdot \boldsymbol{r}_{\bf 次} }{\boldsymbol{p} \cdot A\boldsymbol{p}} \frac{\boldsymbol{p} \cdot A\boldsymbol{p}}{\boldsymbol{r} \cdot\boldsymbol{r}} =\color{red}{\frac{\boldsymbol{r}_{\bf 次} \cdot \boldsymbol{r}_{\bf 次} } {\boldsymbol{r} \cdot \boldsymbol{r} }}$$

  • ループ内の計算順序と,各変数の計算に用いる変数の関係をまとめると上のようになる
  • αxr,条件判定,βp の順で計算する

まとめ


  • 勾配法は誤差関数が小さくなるように近似解を少しずつ更新する方法である
  • 勾配法には,誤差関数,移動方向の計算方法,移動距離の計算方法,初期値,打切り誤差の設定が必要である
  • 共役勾配法(CG法)は正定値対称行列を係数行列とする連立一次方程式を効率的に解くことができる手法である
  • CG法は初期値と打切り誤差以外は明確に設定されている
  • CG法の導出は高度であるが,アルゴリズムは洗練されて単純なため実装はかなり容易である

ちょっとしたコメント


  • 本当はもう少し続けるはずだったんですが,長くなってしまったので今回はここまでにします
  • CG法では正定値対称行列に特化していますが,まだFEMに表れる行列の特徴を網羅できていません
  • FEMの行列のもうひとつの特徴はスパース行列であることです
  • ということで次回はスパース正定値対称行列に特化した手法を紹介します

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です