初めに

線形な連立方程式$Ax = b$（linear system of equations）を満たす$x$を求めよ，という問題は様々な場面で登場するので，そのような問題を解く方法たちの勉強をしよう．

以下が勉強になる:

Methods

$Ax = b$の解法には，直接法と反復法という分類がある．

直接法: 丸め誤差が無い理想的な条件の下で，有限回の演算で解を与える方法を指す, e.g., Gaussの消去法．
反復法: 適当な漸化式により列を生成し，列の収束で以て解を与える方法を指す．
- 定常法　: 反復過程で$x$以外に変化する変数が無いもの, e.g., Jacobiの反復法．
- 非定常法: 反復過程で$x$以外に変化する変数が有るもの, e.g., 共役勾配法．

最も有名な解法の一つである共役勾配法は，解の更新の形式（$x \leftarrow x + \Delta x$）から反復法に分類されるが，有限回の反復の内に収束することが示されているので，この視点からは直接法と反復法の両者の性質を持つとも言える．

本稿では，反復法を勉強する．

定常法（stationary method）

Jacobi法

$A = D - L - U$と分解する．ただし，$L$は狭義下三角行列，$U$は狭義上三角行列，$D$は対角行列である．

\begin{align}
    A x &= b \\
    (D - L - U) x &= b \\
    D x^{(k+1)} &= (L + U) x^{(k)} + b \\
    \therefore x^{(k+1)} 
    &= D^{-1}(L + U) x^{(k)} + D^{-1}b \\
    &= M_{\text{J}} x^{(k)} + c_{\text{J}} \\
\end{align}

とする．$x^{(k+1)} = M x^{(k)} + c$という形を反復法の一般の表現と呼ぶ．$M$を反復行列と呼び，$M_{\text{J}}$はJacobi法における反復行列である．

真なる解$x^{(*)}$に到達するとき，

\begin{align}
    x^{(*)} = M x^{(*)} + c
\end{align}

を得る（$\lim_{k \to \infty} x^{(k)} = x^{(*)}$）．反復法の一般の表現（$x^{(k+1)} = M x^{(k)} + c$）との差分を取ると，

\begin{align}
    x^{(k+1)} - x^{(*)} 
    &= M (x^{(k)} - x^{(*)}) \\
    e^{(k+1)}
    &= M e^{(k)} \\
    &= M^{k+1} e^{(0)} \\
\end{align}

となり，誤差$e^{(k+1)} = x^{(k+1)} - x^{(*)}$は，1反復前の誤差$e^{(k)}$に対して線形変換$M_{\text{J}} \left( = D^{-1}(L + U) \right)$を作用させたものになっている．$M$のスペクトル半径$\varrho \left( M \right) \left( = \max_{i} | \lambda_{i} | \right)$が$\varrho \left( M \right) \lt 1$を満たすとき（$\lambda_{i}$は$M$の固有値），$x$は収束する（$\lbrace x^{(k)} \rbrace_{k \in \mathbb{N}}$に極限が存在）．これは，今のところ，次のように理解している．$M$の固有値は，$M$が作用するベクトル（今回の文脈では$e^{(k)}$や$e^{(0)}$）をどれだけ引き伸ばすかを意味するので，$M$が誤差が成長させるような引き伸ばしの作用を持っている場合，誤差が成長していき，発散する．スペクトル半径が1より小さいとき，誤差が減衰するから，収束する．また，スペクトル半径が小さい方が，誤差の減衰が速いから，収束が速い．

Damped Jacobi法

Jacobi法を再掲する．

\begin{align}
    x^{(k+1)} &= D^{-1}(L + U) x^{(k)} + D^{-1}b \\
\end{align}

ここで，$\omega = 1$なるパラメータを導入して，

\begin{align}
    x^{(k+1)} 
    &= \omega D^{-1}(L + U) x^{(k)} + \omega D^{-1}b + (1 - \omega) x^{(k)} \\
    &= (\omega D^{-1}(L + U) + (1 - \omega) I) x^{(k)} + \omega D^{-1}b \\
    &= M_{\text{DJ}} x^{(k)} + c_{\text{DJ}} \\
\end{align}

と書き改める．$\omega = 1$を選べば，Jacobi法に一致する．任意の$\omega$を認めることにすれば，それがDambed Jacobi法である．$\omega$を減衰 / 加速 / 重みパラメータなどと呼ぶ．

Gauss-Seidel法

先ほどと同じく，$A = D - L - U$から，

\begin{align}
    A x &= b \\
    (D - L - U) x &= b \\
    (D - L) x^{(k+1)} &= U x^{(k)} + b \\
    \therefore x^{(k+1)} 
    &= (D - L)^{-1} U x^{(k)} + (D - L)^{-1} b \\
    &= M_{\text{GS}} x^{(k)} + c_{\text{GS}} \\
\end{align}

$M_{\text{GS}}$はGauss-Seidel法における反復行列である．

SOR（Successive Over-Relaxation）法

Gauss-Seidel法を再掲する．

\begin{align}
    A x &= b \\
    (D - L) x^{(k+1)} &= U x^{(k)} + b \\
\end{align}

ここで，$\omega = 1$なるパラメータを導入して，

\begin{align}
    (D - \omega L) x^{(k+1)} 
    &= ((\omega - 1) D + \omega U) x^{(k)} + \omega b \\
    \therefore
    x^{(k+1)} 
    &= (D - \omega L)^{-1} ((\omega - 1) D + \omega U) x^{(k)} + \omega (D - \omega L)^{-1} b \\
    &= M_{\text{SOR}} x^{(k)} + c_{\text{SOR}} \\
\end{align}

と書き改める．$\omega = 1$を選べば，$M_{\text{SOR}} = M_{\text{GS}}, \ c_{\text{SOR}} = c_{\text{GS}}$より，Gauss-Seidel法に戻る．$\omega$を自由に選ぶことにすれば，これがSOR（Successive Over-Relaxation，逐次加速緩和）法である．$\omega (\in (0, 2))$を緩和 / 加速 / 重みパラメータなどと呼ぶ．一般的には，$\omega \lt 1$でGauss-Seidel法より遅く，$\omega \gt 1$でGauss-Seidel法より速い．

非定常法（non-stationary method）

最急降下（Steepest Descent）法

SD法（Steepest Descent method，最急降下法）は最も明快な考え方の一つだろう．

連立方程式

\begin{align}
    A x = b
\end{align}

の求解を，2次式

\begin{align}
    \mathcal{J} (x)
    = \frac{1}{2} (x, A x) - (x, b) + c
\end{align}

の最小化問題に置き換える（$c \ge 0$は任意の定数だが，$x$について微分すれば消えてなくなるので，以降は省略する）．$A$は対称で正定値（Symmetric Positive Definite; SPD）な行列なので，$\mathcal{J} (x)$は下向きに凸な（綺麗な）お椀型の関数である．すなわち，$x^{(*)} = \arg \min_{x} \mathcal{J} (x)$がただ一つだけ在って，それは$\partial_{x} \mathcal{J} (x) = 0$を与えるような$x$である（最小二乗法のときの議論に似ている）．

反復法のイメージを忘れないでおくと，$k$回反復したときの近似解$x^{(k)}$を$\Delta x$だけ修正して$k+1$回目の近似解$x^{(k+1)}$を作る．

\begin{align}
    x^{(k+1)}
    &= x^{(k)} + \Delta x \\
    &= x^{(k)} + \alpha^{(k)} p^{(k)}
    \left( = x^{(k)} + \sum_{t=0}^{k} \alpha^{(t)} p^{(t)} \right) \\
\end{align}

$\Delta x = \alpha^{(k)} p^{(k)}$とする．$\alpha^{(k)}$は修正量（スカラー），$p^{(k)}$は修正方向（ベクトル）である¹．$p^{(k)}$で歩く方向を決めておき，$\alpha^{(k)}$でどれだけ歩くか決める，という方針である（line search）．

さて，$p^{(k)}$にはどのような方向を選ぼうか．素朴に考えれば，$p^{(k)} = - \nabla_{x} \mathcal{J}$が思い付く²．

\begin{align}
    \mathcal{J} \left( x^{(k)} \right)
    &= \frac{1}{2} (x^{(k)}, A x^{(k)}) - (x^{(k)}, b) \\
    \therefore
    p^{(k)}
    :=
    - \nabla_{x} \mathcal{J} \left( x^{(k)} \right)
    &= - ( A x^{(k)} - b )
    = r^{(k)} \\
\end{align}

すなわち，残差ベクトルの方向に解を更新せよ，と言っている．解の更新に残差ベクトルを利用する方針は，Jacobi法，Gauss-Seidel法，SOR法いずれにも無かった，非常に興味深い作戦である．

ところで，$k+1$回目の近似解$x^{(k+1)}$はどれくらい小さな$\mathcal{J}$を達成するのだろう？

\begin{align}
    \mathcal{J} \left( x^{(k+1)} \right)
    &= \frac{1}{2} \left( x^{(k+1)}, A x^{(k+1)} \right)
        - \left( x^{(k+1)}, b \right) \\
    &= \frac{1}{2} \left( x^{(k)} + \Delta x, A (x^{(k)} + \Delta x) \right)
        - \left( x^{(k)} + \Delta x, b \right) \\
    &= \frac{1}{2} (x^{(k)}, A x^{(k)}) - (x^{(k)}, b) 
    - (\Delta x, b - A x^{(k)})
    + \frac{1}{2} (\Delta x, A \Delta x) \\
\end{align}

ここで，

\begin{align}
    \mathcal{J} \left( x^{(k)} \right)
    &= \frac{1}{2} (x^{(k)}, A x^{(k)}) - (x^{(k)}, b) \\
    r^{(k)}
    &= b - A x^{(k)} \\
    \Delta x
    &= \alpha^{(k)} p^{(k)} \\
\end{align}

に注意すれば，

\begin{align}
    \mathcal{J} \left( x^{(k)} + \alpha^{(k)} p^{(k)} \right)
    &= \mathcal{J} \left( x^{(k)} \right)
    - \alpha^{(k)} (p^{(k)}, r^{(k)})
    + \frac{(\alpha^{(k)})^2}{2} (p^{(k)}, A p^{(k)}) \\
\end{align}

より，$\alpha^{(k)}$の2次式になっている．かつ右辺第3項の符号から，$\mathcal{J} \left( x^{(k)} + \alpha^{(k)} p^{(k)} \right)$は$\alpha^{(k)}$について下に凸である．これより，最も効果的な歩幅は，$\partial_{\alpha^{(k)}} \mathcal{J} \left( x^{(k)} + \alpha^{(k)} p^{(k)} \right) = 0$を与える$\alpha^{(k)}$であると分かる（そのような$\alpha$が，1反復の間に，最も素早く$\mathcal{J}$を減少させる，これも最小二乗法の議論に似ている）．即ち，

\begin{align}
    0
    =
    \frac{\partial}{\partial \alpha^{(k)}}
    \left( \mathcal{J} \left( x^{(k)} + \alpha^{(k)} p^{(k)} \right) \right)
    &= - (p^{(k)}, r^{(k)})
    + \alpha^{(k)} (p^{(k)}, A p^{(k)}) \\
    \therefore
    \alpha^{(k)}
    &= \frac{(p^{(k)}, r^{(k)})}{(p^{(k)}, A p^{(k)})} \\
\end{align}

を歩幅に選ぶ．

以上より，initial guess$x^{(0)}$を用意しておいて，

$x^{(k)} \leftarrow x^{(0)}$
$r^{(k)} \leftarrow b - A x^{(k)}$
$p^{(k)} \leftarrow r^{(k)}$
$\alpha^{(k)} \leftarrow (p^{(k)}, r^{(k)}) / (p^{(k)}, A p^{(k)})$
$x^{(k+1)} \leftarrow x^{(k)} + \alpha^{(k)} p^{(k)}$
$r^{(k+1)} \leftarrow r^{(k)} - \alpha^{(k)} A r^{(k)}$
$\text{if norm } (r^{(k+1)}) \lt \epsilon \text{ break, else set $k \leftarrow k + 1$, go back to step 3 and repeat}$

を実行する（大雑把な書き方だが）．

共役方向（Conjugate Direction）法

SD法は明快だが，効率が悪い．これは，各反復で降下する方向が常に直交しているためである．例えば，楕円状の$\mathcal{J}$などを考えれば，解の軌跡が振動するから，この方法は中々苦労しそうだとイメージできる³．降下方向$p$を考え直そう．CD法（Conjugate Direction method，共役方向法，例えば，Hestenes1980）を導入する．

先ず，誤差$s^{(k)} = A^{-1} b - x^{(k)} = x^{(*)} - x^{(k)} \left( \Leftrightarrow A s = b - A x = r \right)$を用意する．議論の出発点は，

\begin{align}
    \left( s^{(k+1)}, A p^{(k)} \right) 
    &= \left( A s^{(k+1)}, p^{(k)} \right) \quad (\because A = A^{\top}) \\
    &= \left( r^{(k+1)}, p^{(k)} \right) \\
    &= \left( r^{(k)} - A \Delta x, p^{(k)} \right) \\
    &= \left( r^{(k)} - \alpha^{(k)} A p^{(k)}, p^{(k)} \right) \\
    &= \left( r^{(k)}, p^{(k)} \right) 
        - \alpha^{(k)} \left( A p^{(k)}, p^{(k)} \right) \\
    &= \left( r^{(k)}, p^{(k)} \right) 
        - \frac{(p^{(k)}, r^{(k)})}{(p^{(k)}, A p^{(k)})} \left( A p^{(k)}, p^{(k)} \right) \\
    &= 0 \\
\end{align}

なる等式である．定義から，$s^{(k+1)}$は，$x^{(k+1)}$から見てどこに$x^{(*)}$が居るかを示している．一方，$A p^{(k)}$は既知である．ここから（$s^{(k+1)}$が具体的に分かるかは置いておいて⁴），既知量$A p^{(k)}$に直交する方向を探索すると何か良いことがありそうだ，と感じる．

ここで，用語の定義を記しておく．

直交
2つの非ゼロなベクトル$u$と$v$が$\left( u, v \right) = 0$を満たすとき，$u$と$v$は直交するとか，直交であるなどとという．

共役
2つの非ゼロなベクトル$u$と$v$が$\left( u, A v \right) = 0$を満たすとき（$u^{\top} A v = 0$と書くことの方が多い？），$u$と$v$はAについて共役とか，$A$-共役 / $A$-直交であるなどという．$A$が$I$ならば，所謂通常の直交に戻る．

再掲すると，

\begin{align}
    \left( s^{(k+1)}, A p^{(k)} \right)
    = 0 \\
\end{align}

そう言えば，SD法から脱却するために，降下方向を改めるのだった．

先ず，いまの場所を$x^{(k)}$として，ここからの降下方向$p^{(k)}$を考える．適当なベクトル$t^{(k)}$を用意し，これを降下方向として（$p^{(k)} = t^{(k)}$），$x^{(k+1)} (= x^{(k)} + \alpha^{(k)} p^{(k)})$を作る．

次に，新しい場所$x^{(k+1)}$からの降下方向$p^{(k+1)}$を考える．適当なベクトル$t^{(k+1)}$を用意し，これを降下方向としよう（$p^{(k+1)} = t^{(k+1)}$）として，ちょっと考える．$\left( s^{(k+1)}, A p^{(k)} \right) = 0$を思い出し，$t^{(k+1)}$を$p^{(k)}$と$A$-共役な方向（$A$について直交する方向）に射影する．

\begin{align}
    p^{(k+1)}
    = t^{(k+1)} 
    - \frac{\left( t^{(k+1)}, A p^{(k)} \right)}{\left( p^{(k)}, A p^{(k)} \right)} p^{(k)}
\end{align}

ここで，念の為確認しておくと，

\begin{align}
    \left( p^{(k+1)}, A p^{(k)} \right)
    &= \left( t^{(k+1)} - \frac{\left( t^{(k+1)}, A p^{(k)} \right)}{\left( p^{(k)}, A p^{(k)} \right)} p^{(k)}, A p^{(k)} \right) \\
    &= \left( t^{(k+1)}, A p^{(k)} \right)
        - \frac{\left( t^{(k+1)}, A p^{(k)} \right)}{\left( p^{(k)}, A p^{(k)} \right)} 
            \left( p^{(k)}, A p^{(k)} \right) \\
    &= 0
\end{align}

より，$p^{(k+1)}$と$p^{(k)}$は$A$-共役になっている．

続いて，$x^{(k+2)}$からの降下方向$p^{(k+2)}$を考える．再度，適当なベクトル$t^{(k+2)}$を用意し，これを降下方向としよう（$p^{(k+2)} = t^{(k+2)}$）として，考え直す．$t^{(k+2)}$を$p^{(k)}$と$A$-共役，かつ，$p^{(k+1)}$と$A$-共役な方向に射影する．

\begin{align}
    p^{(k+2)}
    &= t^{(k+1)} 
    - \frac{\left( t^{(k+2)}, A p^{(k+1)} \right)}{\left( p^{(k+1)}, A p^{(k+1)} \right)} p^{(k+1)}
    - \frac{\left( t^{(k+1)}, A p^{(k)} \right)}{\left( p^{(k)}, A p^{(k)} \right)} p^{(k)} \\
\end{align}

これは，実は0回目の反復から続いていた．つまり，

\begin{align}
    p^{(k+1)}
    = t^{(k+1)} 
    - \sum_{i = 0}^{k} \frac{\left( t^{(k+1)}, A p^{(i)} \right)}{\left( p^{(i)}, A p^{(i)} \right)} p^{(i)}
\end{align}

が一般化された降下方向$p^{(k+1)}$の式である．即ち，常に，今までに登場した降下方向$p$と$A$-共役な方向を探し続ける（$\left( p^{(k+1)}, A p^{(i)} \right) = 0 \text{ for } i \le k$）．これにより，効果的に探索範囲が絞られてゆく．その結果，1回の更新ごとに解が存在する範囲が1次元ずつ狭まってゆき，$K$回反復したとき，解が存在する範囲は残り0次元となっている．即ち，厳密解に到達している．

以上より，適当なinitial guess$x^{(0)}$の下で，

$x^{(k)} \leftarrow x^{(0)}$
$r^{(k)} \leftarrow b - A x^{(k)}$
$p^{(k)} \leftarrow r^{(k)}$
$\alpha^{(k)} \leftarrow (p^{(k)}, r^{(k)}) / (p^{(k)}, A p^{(k)})$
$x^{(k+1)} \leftarrow x^{(k)} + \alpha^{(k)} p^{(k)}$
$r^{(k+1)} \leftarrow r^{(k)} - \alpha^{(k)} A p^{(k)}$
$\text{if norm } (r^{(k+1)}) \lt \epsilon \text{ break}$
$\text{take } t^{(k+1)} \text{ so is linearly independent of } \lbrace t^{(i)} \rbrace_{i = 0}^{k}$
$p^{(k+1)} \leftarrow t^{(k+1)} - \sum_{i = 0}^{k} \left( t^{(k+1)}, A p^{(i)} \right) / \left( p^{(i)}, A p^{(i)} \right) p^{(i)}$
$\text{set $k \leftarrow k + 1$, go back to step 4 and repeat}$

を実行する．

上記のように，CD法は，$K$元の方程式に対して高々$K$回で厳密解に届くが，メモリの要求が大きく（$p^{(0)}, p^{(1)}, \cdots, p^{(K-1)}$を格納する必要がある），途中の$x^{(k)} (k \lt K)$が近似解となっている保証が無い．これより，直接法に対する優位性（メモリ容量の少なさ）と，反復法としての有用性（途中で打ち切っても近似解を得る）が弱い．

共役勾配（Conjugate Gradient）法

CG法（Conjugate Gradient method，共役勾配法，Hestenes&Stiefel1952）は，現代で最もポピュラーな反復法の一つと言える⁵（SPDなら事実上最強？）．

やはり，降下方向にもう少しの工夫が欲しい．具体的には，次の降下方向である

\begin{align}
    p^{(k+1)}
    = t^{(k+1)} 
    - \sum_{i = 0}^{k} \frac{\left( t^{(k+1)}, A p^{(i)} \right)}{\left( p^{(i)}, A p^{(i)} \right)} p^{(i)}
\end{align}

をもう少し簡略化したい．このために，適当に選んでいた$t$の代わりに$r$を選ぶ．$r^{(k)}$は，

\begin{align}
    \mathcal{J} \left( x^{(k)} \right)
    &= \frac{1}{2} (x^{(k)}, A x^{(k)}) - (x^{(k)}, b) \\
    \therefore
    \frac{\partial}{\partial x^{(k)}} \left( \mathcal{J} \left( x^{(k)} \right) \right)
    &= - \left( b - A x^{(k)} \right) \\
    &= - r^{(k)} \\
\end{align}

より，$\mathcal{J}$のローカルな勾配として定まる（実は最急降下法はこれだった）．したがって，次の降下方向を

\begin{align}
    p^{(k+1)}
    = r^{(k+1)} 
    - \sum_{i = 0}^{k} \frac{\left( r^{(k+1)}, A p^{(i)} \right)}{\left( p^{(i)}, A p^{(i)} \right)} p^{(i)}
\end{align}

に定める．これだけではCD法とあまり変わらないように見えるが，$\left( r^{(k+1)}, A p^{(i)} \right) = 0 \left( \text{for } i \le k-1 \right)$を示すことができる（らしい）．これより，

\begin{align}
    p^{(k+1)}
    &= r^{(k+1)} 
    - \sum_{i = 0}^{k} \frac{\left( r^{(k+1)}, A p^{(i)} \right)}{\left( p^{(i)}, A p^{(i)} \right)} p^{(i)} \\
    &= r^{(k+1)} 
    - \frac{\left( r^{(k+1)}, A p^{(k)} \right)}{\left( p^{(k)}, A p^{(k)} \right)} p^{(k)}
\end{align}

が結論付けられる．過去の$p$を覚えておく必要が無くなり，嬉しい．また，最急降下法のエッセンスを取り込んで$r$を利用することにしている（ローカルな最急降下方向）ので，途中で打ち切ってもある近似になっていると考えられる．

また，

\begin{align}
    \alpha^{(k)}
    &:= \frac{\left( r^{(k)}, p^{(k)} \right)}
            {\left( p^{(k)}, A p^{(k)} \right)}
     = \frac{\left( r^{(k)}, r^{(k)} \right)}
            {\left( p^{(k)}, A p^{(k)} \right)} \\
    \beta^{(k)}
    &:= - \frac{\left( r^{(k+1)}, p^{(k)} \right)}
            {\left( p^{(k)}, A p^{(k)} \right)}
     = \frac{\left( r^{(k+1)}, r^{(k+1)} \right)}
            {\left( r^{(k)}, r^{(k)} \right)} \\
\end{align}

を示すことができる（らしい）．

以上より，適当なinitial guess$x^{(0)}$の下で，

$x^{(k)} \leftarrow x^{(0)}$
$r^{(k)} \leftarrow b - A x^{(k)}$
$p^{(k)} \leftarrow r^{(k)}$
$\alpha^{(k)} \leftarrow (r^{(k)}, r^{(k)}) / (p^{(k)}, A p^{(k)})$
$x^{(k+1)} \leftarrow x^{(k)} + \alpha^{(k)} p^{(k)}$
$r^{(k+1)} \leftarrow r^{(k)} - \alpha^{(k)} A p^{(k)}$
$\text{if norm } (r^{(k+1)}) \lt \epsilon \text{ break}$
$\beta^{(k)} \leftarrow (r^{(k+1)}, r^{(k+1)}) / (r^{(k)}, r^{(k)})$
$p^{(k+1)} \leftarrow r^{(k+1)} + \beta^{(k)} p^{(k)}$
$\text{set $k \leftarrow k + 1$, go back to step 4 and repeat}$

を実行する．

なお，

\begin{align}
    r^{(k)} &= b - A x^{(k)} \\
    r^{(k+1)} &= b - A x^{(k+1)} \\
    \therefore
    r^{(k+1)}
    &= r^{(k)} - A \left( x^{(k+1)} - x^{(k)} \right) \\
    &= r^{(k)} - \alpha^{(k)} A p^{(k)} \\
\end{align}

を利用して残差を評価している．

Results

各手法を自作で実装した．CG法は，scipy.sparse.linalg.cgと比較して正しく作れた（であろう）ことを確認している．

Problem 1

1次元の棒（rod）におけるLaplace方程式のDirichlet問題を考える（と言うか，2階のODE?）．

\begin{alignat}{2}
    - \frac{d^2}{d x^2} u &= 0 \quad &&\text{ } \quad (x \in (0, 1)) \\
    u &= 0 \quad &&\text{ } \quad (x = 0) \\
    u &= 1 \quad &&\text{ } \quad (x = 1) \\
\end{alignat}

0度の棒の右端を1度に熱し，その温度が棒全体に伝播するまで十分な時間待つ．このとき，左端は0度に固定する．解は線形な温度分布となる．

未知数の個数（方程式の元）$N$を変えながら各手法で求解し，残差$r^{(k)} = b - A x^{(k)}$の$L^{\infty}$ノルムが$\epsilon \left( = 10^{-9} \right)$より小さくなったときに収束判定した．$L^{\infty}$ノルムを選んだのは，$L^{1}$や$L^{2}$では，ベクトルの次元に応じてノルムが大きくなるためである．最大の反復回数は$4N$とした．

$N$	$N=16$	$N=32$	$N=64$
$r$
$u$

残差の収束過程（1行目）から，CG法以外の方法が収束判定を満たす前に最大反復回数に達している．また，解の分布（2行目）から，各方法が解をSmoothにはしているものの，CG法以外の方法は線形な解を表現するに達していない．ただし，最大の反復回数を大きくすれば，全ての方法が収束判定を満たし，解関数が同一となることを確認している．Jacobi法とSD法が与える解が全く重なっているのが興味深い．CG法では，反復回数が未知数の個数$N$に等しくなったとき，残差ノルムが$\mathcal{O} (10^{-16})$に落ちた．有限回の反復で収束する性質が確認できる．

ところで，Jacobi法やGauss-Seidel法は，収束しておらず，解も線形になっていない（そうなるまで待ち切れなかった，というのが正しいか？）．しかし，「解が滑らかになっている」ことに注目したい．残差ノルムの動きを見ると，初めの数十回の反復の間，Jacobi/Gauss-Seidel法は残差を急速に減少させるものの，その後の残差の減少傾向は非常に緩やかになる．これは，残差場の高周波成分が急速に平滑化（スムージング）され，低周波成分がしぶとく残っていることを意味している．この性質を積極的に活用しよう，という考えがマルチグリッド法（例えば，Brandt1977, Stüben&Trottenberg1982, Trottenberg+2000）である．すなわち，数~数十回のJacobi / Gauss-Seidel反復を行い，空間を粗視化する．すると，先程までの細かい空間で低周波だった残差が，粗い空間では高周波になっている．では，この空間でJacobi / Gauss-Seidel反復を実行すると，また低周波数成分が残るので，再度粗視化しよう，，，これを繰り返す．注意として，Aliasingが起こるのを避けるため，各レベルでJacobi / Gauss-Seidel反復を行う必要がある（急激な粗視化は不安定）．ただし，各レベルでのJacobi / Gauss-Seidel法は，厳密求解に到達していない．この意味で，各レベルで用いたJacobi / Gauss-Seidel法は，"ソルバー"ではなく，"スムーザー"と呼ばれる（ここの議論は，この動画が分かり易い）．また，マルチグリッド法は，最も粗い空間でもJacobi / Gauss-Seidel法を用いる（ここでは，"ソルバー"として）方法であるが，これを，例えばCG法などの非定常ソルバーに置き換えるものはマルチグリッド前処理付きCG法（Multi-Grid preconditioned Conjugate Gradient method; MGCG method, e.g., Tatebe1993）などと呼ばれる．

Problem 2

ShermanのOil reservoir simulation challenge matrices（Sherman1984）を考える．3次元有限差分シミュレーションから用意された問題である．対象な問題である，Sherman1を取り上げる．Matrix MarketやSparse Matrix Collectionから取得できる．

$A$	$b$

条件数は，

$\kappa_{1} = 2.26 \times 10^4$
$\kappa_{2} = 1.56 \times 10^4$
$\kappa_{\infty} = 2.26 \times 10^4$

である．戸川1977に従うと，悪条件に分類される（はず）．

$r$	$x$

正直，$x$の物理的な意味が良く分かっていないので（3次元にreshapeしないと元の現象が分からないのだろうが，問題設定も良く分かっていないので），残差ノルムの収束過程だけ見る（こういうことをするのは良くない気がするが，論文ではないので大らかに，，，）．大雑把に悪条件の括りに入っている（多分？）が，CG法は1,000次元の問題を283回の反復の間に収束させている．戸川1977には，

共役勾配法によれば有限回（$n$元の方程式の場合，高々$n$回）で厳密解を得ることができるが，実際には，$n$よりももっとはるかに少ない回数で十分な精度の解が得られる場合が少なくない．

という記述があるが，これは確からしい．toleranceが$10^{-9}$なのは少し優しい気がしたので，$10^{-15}$にして再計算したが，$536$回で収束した．やはり，元の数より小さい．

終わりに

興味深い内容であった．

本当はKrylov部分空間というものを勉強して，何故Krylov部分空間法という括りがあるのかまで勉強したかったが，宿題に残しておく．

Appendices

A. 取り零し

共役勾配法に取り零しがある．

A-1. 残差の直交性

$\left( r^{(k+1)}, r^{(k)} \right) = 0$を示す．

$\left( r^{(k+1)}, r^{(k)} \right) = 0$から，最急降下法が，楕円状の$\mathcal{J}$を振動する様子が改めて良く分かる．

（宿題）

A-2. 残差と降下方向の共役性

$\left( r^{(k+1)}, A p^{(i)} \right) = 0 \left( \text{for } i \le k-1 \right)$を示す．

（宿題）

A-3. 係数の書き換え

$\alpha$，$\beta$を示す．

（宿題）

B. 条件数

正則な行列$A$に対して，

\begin{align}
    \kappa_{p} = \lVert A \rVert_{p} \lVert A^{-1} \rVert_{p}
\end{align}

を$p$ノルムに関する$A$の条件数（condition number）とか，$A$の$p$-条件数などと呼ぶ．特に，

\begin{align}
    \kappa_{2} = \lambda_{\text{max}} / \lambda_{\text{min}}
\end{align}

のことを単に$\kappa$や条件数などと呼ぶこともある．ただし，$\lambda_{i}$は$A$の固有値である．収束率と固有値には密接な関係がある（vanderSluis&vanderVorst1986）．

条件数$\kappa$は，収束性と密接に関係しており，1に近いほど都合が良い．$\kappa_{p}$については未だ良く分かっていないのだが，話を$\kappa_{2}$に限れば，これは最大固有値と最小固有値の比率だから，各軸方向への曲率の比のうち，最も大きいものを示している．即ち，$\kappa_{2}$は$\mathcal{J}$が形成する楕円の細長さである．こう考えれば，そりゃあ条件数が1に近いほど問題が解き易いところまでは解釈できる．

戸川1977によると，凡その目安として$\kappa_{2} \gt 10^{4}$程度から悪条件という括りに入るらしい．

C. 都合の良い緩和パラメータはどの程度か？

Damped Jacobi法とSOR法の緩和パラメータ$\omega$の値の範囲を検討する．最適な$\omega$を事前に知る方法もあるらしいが，それができる問題も限られているらしいので，適当に都合の良い範囲を数値的に検討する．

C-1. 非対称

Gilbert Strang's Final Lectureから，

\begin{align}
    A &=
    \begin{bmatrix}
        2 &  3 &  4 \\
        4 & 11 & 14 \\
        2 &  8 & 17 \\
    \end{bmatrix}
    , \
    b = (19, 55, 50)^{\top}
    , \
    x^{*} = (4, 1, 2)^{\top} \\
\end{align}

を考える．条件数は，

$\kappa_{1} = 62.50$
$\kappa_{2} = 34.57$
$\kappa_{\infty} = 39.77$

である．

Damped Jacobi	SOR

Damped Jacobi法の傾向として，$\omega \gt 1$では，全てスペクトル半径が1を超え，発散した．杉原&室田2009は$\omega \in (0, 1]$が定義と言っているので，ルールを破ってしまったのだと思う．$\omega = 1$（即ち，Jacobi法）が発散したのは興味深い（勝手にタフな方法だと思っていたが，そういうわけでもないか？）．これより，Damped Jacobi法は安定化の作用があると見ることもできる（？）．$\omega \in (0, 1)$の範囲では，$\omega = 0.8$がスペクトル半径が最も小さく，速い．今回の問題では，$\omega = 0.8$が最も都合が良いだろうか．

SOR法の傾向として，$\omega \lt 1$でGauss-Seidel法（$\omega = 1$）より遅い．$\omega \gt 1$でGauss-Seidel法より速い．が，大きければ大きいほど良い訳ではなく，$\omega = 1.8$では寧ろ遅い．$\omega = 1.6$も速くはあるが，振動があり，安定性に不安が残る．今回の問題では，$\omega = 1.4$程度が安定性と高速性のバランスが良さそうだ（スペクトル半径も最も小さい）．

C-2. 対称

ShermanのOil reservoir simulation challenge matrices（Sherman1984）を考える．

再掲だが，条件数は，

$\kappa_{1} = 2.26 \times 10^4$
$\kappa_{2} = 1.56 \times 10^4$
$\kappa_{\infty} = 2.26 \times 10^4$

である．戸川1977に従うと，悪条件に分類される（はず）．

Damped Jacobi	SOR

Damped Jacobi法は，先ほどと同じく，$\omega \gt 1$ではスペクトル半径が1を超えて発散した．$\omega = 0.6, \ 0.8, \ 1.0$程度が速くはあるものの，そもそも収束性が悪すぎて使い物にならない．

SOR法も，先ほどと同様に，$\omega \lt 1$でGauss-Seidel法より遅い．$\omega \gt 1$でGauss-Seidel法より速い．今回は，$\omega$が大きければ大きいほど速い．

2つの実験を合わせて，えいやっと

DJ: $\omega = 0.6 \sim 0.8$
SOR: $\omega = 1.4 \sim 1.6$

周辺を，「ある程度versatileなもの」という緩いルールで選んだ（本文ではこの中から選んだ）．

D. CG法のVariants

$A$が非対称なとき，短い再帰の内に残差ベクトルを直交化することが難しいため，CG法には都合が悪い（Faber&Manteuffel1984）．

GMRES法（Generalized Minimal Residual method，一般化最小残差法，Saad&Schultz1986）は長い再帰を利用して残差ベクトルを直交化するが，メモリ消費が大きい．

非対称の問題には，

BiCG法（Bi-Conjugate Gradient method，双共役勾配法，Fletcher1976）
CGS法（Conjugate Gradient Squared method，自乗共役勾配法，Sonneveld1989）
BiCGSTAB法（Bi-Conjugate Gradient STABilized method，安定化双共役勾配法，vanderVorst1992）

などが提案されてきた（まあ，もっともっと沢山あるわけだが（BiCGSTAB-(2), BiCGSTAB-($l$), etc.），，，）．どれが良いか？というのが気になるが，Gutknecht1993やdePillis1998が数値実験を報告している．

E. ベクトルの射影

非ゼロのベクトル$u$と$v$を用意する．一般には，$u$と$v$は直交していない．

良く知られた射影のルール（正式な名前を知らない，，， ~~そしてHelmholtzの分解の定理に似ている気がするのは気のせいだろうか，，，~~ いや，違うな）から，

\begin{align}
    u &= u_{\parallel v} + u_{\perp v}
\end{align}

なる分解が可能である．ここで，$u_{\parallel v}$は$u$の内$v$に平行な成分，$u_{\perp v}$は$u$の内$v$に垂直な成分である．なお，

\begin{align}
    u_{\parallel v}
    &= \left( u, \hat{v} \right) \hat{v} \\
    &= \frac{\left( u, v \right)}{\| v \|} \frac{v}{\| v \|} \\
    &= \frac{\left( u, v \right)}{\left( v, v \right)} v \\
\end{align}

である．

これを利用すれば，

\begin{align}
    u_{\perp v}
    &= u - u_{\parallel v} \\
    &= u - \frac{\left( u, v \right)}{\left( v, v \right)} v \\
\end{align}

などとして，直交なベクトルを取り出すことができる（この関係は，例えば，Yu+2020などで利用されている）．

F. 気になること

関数

\begin{align}
    \mathcal{J} (x)
    = \frac{1}{2} (x, A x) - (x, b)
\end{align}

の$x$を$h$だけ動かすとどうなるか？（$x = x^{(*)}$とする）

\begin{align}
    \mathcal{J} (x + h)
    &= \frac{1}{2} ((x+h), A (x+h)) - ((x+h), b) \\
    &= \frac{1}{2} (x, A x)
        + \frac{1}{2} (x, A h) 
        + \frac{1}{2} (h, A x) 
        + \frac{1}{2} (h, A h) 
        - (x, b) 
        - (h, b) \\
    &= \mathcal{J} (x)
        + (h, A x - b) + \frac{1}{2} (h, A h) \\
\end{align}

$A$が対象なので，$(x, A h) = (h, A x)$を用いた．また，$A$は正定値なので，$(h, A h) \gt 0$．したがって，$h \neq 0$に対して，

\begin{align}
    \mathcal{J} (x + h) \gt \mathcal{J} (x)
\end{align}

である．

このタイミングで，ここまでベクトルをボールド体にせずに書いてきたことを酷く後悔したが，書き直すのも酷く面倒に感じてしまった． ↩
機械学習における勾配降下法を思い出せ． ↩
機械学習の視点からは，QianのMomentum法（Qian1999），Nesterovの加速法（Nesterov1983）などがremedyとして発展してきた（機械学習におけるこの辺りの話題は，Ruder2016がすっきり纏まっている）． ↩
まあ，分からない．分かるなら，その方向に1歩歩くだけで済む． ↩
歴史的には，1950年代に華々しくデビューしたが，丸め誤差に対する敏感さが指摘されてからは暫く話題が落ち着いた．1960~70年代に非線形最適化問題や疎行列問題に対する解法として再度注目されるようになり（例えば，Reid1972），1980年代からは前処理に支えられ現在の地位を確立した（らしい）． ↩