雑記：反復法の話 - クッキーの日記

おかしい点がありましたらご指摘いただけますと幸いです。

※ キャラクターは架空のものです。

f:id:cookie-box:20190101155733p:plain:w60

連立一次方程式の数値解法にヤコビ法ってありますよね。

f:id:cookie-box:20190101160814p:plain:w60

連立一次方程式を解くアルゴリズムの中でも反復法といわれるものだね。

$Ax=b$ の右辺を $A$ の対角成分だけ取り残して移項して $Dx = - (L + U)x + b$ とし、さらに $D^{-1}$ をかけて $x = - D^{-1}(L + U)x + D^{-1}b$ として、これを漸化式に見立て $x^{(k+1)} = - D^{-1}(L + U)x^{(k)} + D^{-1}b$ としてしまう。適当な初期値 $x^{(0)}$ から出発して、もし値がほぼ変化しなくなったらそれが解である…。

ヤコビ法（ $A = L + D + U$ ， $D$ は対角行列， $L$ は狭義下三角行列， $U$ は狭義上三角行列．）

$x^{(k+1)} = - D^{-1}(L + U)x^{(k)} + D^{-1}b$

手続きは理解できますが、なぜこうしなければならなかったのでしょう。下三角行列 $L$ と上三角行列 $U$ に取り残されてしまった対角行列 $D$ が不憫ではないですか？

そ、そうかな…。

それに、反復法といいますけど、別に反復をしたくてしているわけでもないはずなんです。連立一次方程式を解きたいのなら、いま手元にある適当な $x$ に $\Delta x$ を足して本当の解との誤差を埋めたいのですよね。であれば、 $x$ から進むべき $\Delta x$ は理想的には $\Delta x = -x + A^{-1}b = A^{-1}(-Ax+b)$ であるはずです。

そりゃ本当の解が $A^{-1}b$ だからね。

ええ、 $A^{-1}$ が求まるなら反復する必要なんてないわけです。 $A^{-1}$ を求めることを回避して $\Delta x$ に近いものを得たいんです。そこで、 $\Delta x = A^{-1} (-Ax +b) \fallingdotseq D^{-1} (-Ax + b)$ という雑な近似をしたわけです。「逆行列求めるの面倒だから対角成分だけの行列の逆行列でいいや」ということですね。こうすれば、

ヤコビ法にたどり着くまで（？）

$\begin{split} x + \Delta x &= x + A^{-1}( -Ax + b) \\ &\fallingdotseq x + D^{-1} (-Ax +b) \\ &= x - D^{-1} (L + D + U) x + D^{-1}b \\ &= - D^{-1} (L + U) x + D^{-1}b \end{split}$

と、ヤコビ法に一致するわけです。副部長、どうですか？この私の推理―。

いや、推理とかじゃないからね。

ガウス・ザイデル法は下三角行列と対角行列を取り残して $(D+L)x = - Ux + b$ としてこれを漸化式に見立てたものなんですね。

ガウス・ザイデル法

$x^{(k+1)} = - (D+L)^{-1}Ux^{(k)} + (D+L)^{-1}b$

ヤコビ法と下三角行列の部分が異なりますが、これは、「 $x$ の $i$ 番目の成分を更新するとき、 $i-1$ 番目までの成分は既に更新後のものが計算されているはずだから、もう更新後のものを使ってしまおう」というものなのですよね。なので、ヤコビ法よりも小さいステップで収束することが期待されますが、 $i-1$ 行目の更新を終えないと $i$ 行目の更新ができないために並列計算できなくなっているデメリットもあります。

ヤコビ法	$Dx^{(k+1)} = - (L + U)x^{(k)} + b$
ガウス・ザイデル法	$Dx^{(k+1)} = - L x^{(k+1)} - Ux^{(k)} + b$

「ガウス・ザイデル法にたどり着くまで」は「ヤコビ法にたどり着くまで」の $\fallingdotseq$ のところで $A^{-1}$ を $D^{-1}$ に取り換えるのではなく $(D+L)^{-1}$ に取り換えるという違いだけなので省略しますね。

ガウス・ザイデル法にたどり着くまで（？） ヤコビ法と同様なので省略。

$x + \Delta x \fallingdotseq - (D+L)^{-1}Ux + (D+L)^{-1}b$

手抜きだな…。

SOR法は、一度ガウス・ザイデル法に照準を定めておいて、そちらの方向へ加速パラメータ $\omega$ で近づくというものですね。

ヤコビ法	$Dx^{(k+1)} = - (L + U)x^{(k)} + b$
ガウス・ザイデル法	$Dx^{(k+1)} = - L x^{(k+1)} - Ux^{(k)} + b$
SOR法	$\begin{cases} D\tilde{x}^{(k+1)} = - L x^{(k+1)} - Ux^{(k)} + b \\ x^{(k+1)} = x^{(k)} + \omega (\tilde{x}^{(k+1)} - x^{(k)}) \end{cases}$

ただ $\omega$ で加速させる部分に $x^{(k+1)}$ が含まれていますから、 $x^{(k+1)}$ について解くとややこしい形になりますね。ただ $\omega=1$ のときはガウス・ザイデル法に一致します。SOR法をそのように決めたので当然ですね。

SOR法

$\begin{split} Dx^{(k+1)} &= Dx^{(k)} + \omega (- L x^{(k+1)} - Ux^{(k)} + b - Dx^{(k)}) \\ (D + \omega L) x^{(k+1)} &= \bigl( (1 - \omega)D - \omega U \bigr) x^{(k)} + \omega b \\ x^{(k+1)} &= (D + \omega L)^{-1} \bigl( (1 - \omega)D - \omega U \bigr) x^{(k)} + \omega (D + \omega L)^{-1} b \end{split}$

「SOR法にたどり着くまで」は「ヤコビ法にたどり着くまで」において $\Delta x$ を $\omega$ で加速させて、 $\fallingdotseq$ のところで $A^{-1}$ を $(D + \omega L)^{-1}$ に取り換えると得られますね。単純に $(D + L)^{-1}$ に取り換えるのではない点に注意です。下三角行列 $L$ は既に $\omega$ で加速して更新されているから、といったイメージでしょうか。

SOR法にたどり着くまで（？）

$\begin{split} x + \omega \Delta x &= x + \omega A^{-1} (-Ax + b) \\ &\fallingdotseq x + \omega (D + \omega L)^{-1} (-Ax + b) \\ &= x - \omega (D + \omega L)^{-1} (D + L + U)x + \omega(D+\omega L)^{-1}b \\ &= (D + \omega L)^{-1} \bigl( D - \omega D - \omega U \bigr)x + \omega(D+\omega L)^{-1}b \end{split}$

つづかない