疎行列データのダウンロード(後半のscipyで使える反復法にて使用)
はじめに
$$ % 変数名 \global\def\A{\boldsymbol{A}} \global\def\D{\boldsymbol{D}} \global\def\E{\boldsymbol{E}} \global\def\I{\boldsymbol{I}} \global\def\L{\boldsymbol{L}} \global\def\M{\boldsymbol{M}} \global\def\U{\boldsymbol{U}} \global\def\b{\boldsymbol{b}} \global\def\c{\boldsymbol{c}} \global\def\x{\boldsymbol{x}} \global\def\xk#1{\boldsymbol{x^{({#1})}}} \global\def\y{\boldsymbol{y}} \global\def\P{\boldsymbol{P}} % 数式 \global\def\EQ{\A\x = \b} % n番目の方程式 \global\def\eqn#1{a_{{#1}1} x_{1} + a_{{#1}2} x_{2} + \dots + a_{{#1}n} x_{n} = b_{#1}} % n番目の漸化式 \global\def\rn#1{x^{(k+1)}_{#1} = x^{(k)}_{#1} + \frac{1}{a_{{#1}{#1}}}(b_{#1} - \sum_{i=1}^{n} a_{{#1}i} x^{(k)}_{i})} % 行列・ベクトルの展開表記 \global\def\expA{ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} } \global\def\expL{ \begin{bmatrix} 0 & 0 & \dots & 0 \\ a_{21} & 0 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ a_{n1} & \dots & a_{nn-1} & 0 \end{bmatrix} } \global\def\exL{ \begin{bmatrix} & & \\ \vdots & \ddots & \\ a_{n1} & \dots & \end{bmatrix} } \global\def\expU{ \begin{bmatrix} 0 & a_{12} & \dots & a_{1n} \\ 0 & 0 & \ddots & \vdots \\ \vdots & \ddots & \ddots & a_{n-1n} \\ 0 & \dots & 0 & 0 \end{bmatrix} } \global\def\exU{ \begin{bmatrix} & \dots & a_{1n} \\ & \ddots & \vdots \\ & & \\ \end{bmatrix} } \global\def\expD{ \begin{bmatrix} a_{11} & 0 & \dots & 0 \\ 0 & a_{22} & \dots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & a_{nn} \end{bmatrix} } \global\def\expx{\begin{bmatrix}x_1 \\ x_2 \\ \vdots\\ x_n\end{bmatrix}} \global\def\expxk#1{\begin{bmatrix}x^{({#1})}_1 \\ x^{({#1})}_2 \\ \vdots\\ x^{({#1})}_n\end{bmatrix}} \global\def\exx{\begin{bmatrix}x_1 \\ \vdots\\ x_n\end{bmatrix}} \global\def\expb{\begin{bmatrix}b_1 \\ b_2 \\ \vdots\\ b_n\end{bmatrix}} \global\def\exb{\begin{bmatrix}b_1 \\ \vdots\\ b_n\end{bmatrix}} $$
前回の記事に引き続き、今回は反復法による連立一次方程式(linear system of equations)の解法を見ていきます。
なお、以降を通じて扱う連立一次方程式を$\EQ$と書くことにします。ここで、
$\A = \expA\;,\;\; \b = \expb\;,\;\; \x = \expx\;.$
であり、$\A$と$\b$がそれぞれ既知であるとき、$\x$を求める問題です。
これら連立一次方程式を構成するパーツを
- $\A$: 係数行列(coefficient matrix)
- $\b$: 定数項、定数ベクトル(constant vector)
- $\x$: 解、解ベクトル(solution)
と呼びます。
また、数値計算のプログラムではscipyの疎行列クラスを使用することとします。
ヤコビ法
考え方
2元連立1次方程式で考えます。
$$ % ここだけ定義 行列⇔連立一次方程式の部品 \global\def\Ga{ \begin{bmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{bmatrix} } \global\def\Gx{\begin{bmatrix}x_1 \\ x_2\end{bmatrix}} \global\def\Gb{\begin{bmatrix}b_1 \\ b_2\end{bmatrix}} \global\def\eqi {a_{11} x_{1} + a_{12} x_{2} = b_1} \global\def\eqii {a_{21} x_{1} + a_{22} x_{2} = b_2} $$
$$ \tag{1} {\Ga\Gx = \Gb \Leftrightarrow \begin{cases}\begin{aligned}\eqi\\ \eqii\end{aligned}\end{cases}} $$
$(1)$式を変形して$x_1, x_2$をそれぞれ左辺に持つ式を作ります。
$$ % ここだけ定義 漸化式の前段階 \global\def\eqxi {x_{1} = \frac{1}{a_{11}}(b_1 - a_{12} x_{2})} \global\def\eqxii {x_{2} = \frac{1}{a_{22}}(b_2 - a_{21} x_{1})} $$ $$ \tag{2} \begin{cases}\begin{aligned}\eqxi\\ \eqxii\end{aligned}\end{cases} $$
ここで新たな変数$x^{(k)}_1, x^{(k)}_2\;(k=0,1,\cdots)$を用意し、$x^{(k)}_1, x^{(k)}_2$に関する漸化式を作ります。
$$ % ここだけ定義 漸化式 \global\def\ri {x^{(k+1)}_1 = \frac{1}{a_{11}}(b_1 - a_{12} x^{(k)}_{2})} \global\def\rii {x^{(k+1)}_{2} = \frac{1}{a_{22}}(b_2 - a_{21} x^{(k)}_{1})} $$ $$ \tag{3} \begin{cases}\begin{aligned}\ri\\ \rii\end{aligned}\end{cases} $$
さて、ここで適当な初期値$x^{(0)}_1, x^{(0)}_2$を定めて$k=0,1,\cdots$と順次計算していき、仮に$k \to \infty$で$x^{(k)}_1, x^{(k)}_2$の値がそれぞれ収束したとすると、その値は$(2)$式、すなわち$(1)$式を満たすので解であるということになります。
また、そのとき十分大きな$k$に対する$x^{(k)}_1, x^{(k)}_2$は解の近似値となり得ます。
という訳で、漸化式を順次計算して許容できる誤差で解の近似値を得ようというのがヤコビ法の基本的な考え方です。
実装に向けて
以下のn元連立一次方程式に拡張して考えてみましょう。
$$\begin{cases}\begin{matrix} \eqn{1}\\ \eqn{2}\\ \vdots\\ \eqn{n} \end{matrix}\end{cases} $$
上と同様に漸化式を作るのですが、少し修正を加えて【$x^{(k+1)}_n = x^{(k)}_n + ❓$】のような、前の近似値に修正項を加えていくような形にします。
そうすることで，収束する場合は修正項の大きさも減少していくため情報落ちを起こし， 近似値の変化も小さくなってついには全く変化しなくなることが期待出来，収束の判断がしやすく なる。
そうすると以下の漸化式が得られます。
$$\tag{4}{\begin{cases}\begin{matrix} \rn{1}\\ \rn{2}\\ \vdots\\ \rn{n} \end{matrix}\end{cases}} $$
ここで$x^{(k)}_1, x^{(k)}_2, \cdots, x^{(k)}_n$のベクトルを$\xk{k}$、$\A$の対角行列を$\D$として、$(4)$式を行列で表現すると以下のようになります。
$$\xk{k} = \expxk{k},\;\; D = \expD$$
$$\tag{5}{\xk{k+1} = \xk{k} + \D^{-1}(\b - \A\xk{k})}$$
実装
という訳で実装してみたのが下のjacobi関数です。$\A, \b$を受け取って漸化式の関数を返します。
漸化式の関数は$\x^{(k+1)}$だけではなく残差(residual: $\b - \A\x^{(k)}$)ベクトルのノルムも返すようにしました。
動作確認してみます。
ここではとりあえず初期値は0で始めて、scipy.sparse.linalgのspsolve関数で求めた解程度の精度で反復を打ち切るようにします(残差基準1)。
●諸々の関数定義
せっかく疎行列を使うので10万次元の行列でやってみます。また、参考までにspsolveの実行時間も計っておきます。
残差ベクトルのノルムの推移は以下のようになっています。
spsolveの3倍程度の時間でしょうか？今度は収束の遅い係数行列(対角成分が2.1)で実行してみます。
まぁまぁ時間が掛かった上に、残差ベクトルのノルムの減少が頭打ちになってしまいました。
これはすなわち「収束した」ということで、「収束したとしたらそれは解である」わけなんですが、簡単に実装した限りではこれぐらいの誤差は混入しうるということなんでしょう。
対策としては誤差そのものを減らすような方向と、ある程度の精度で妥協して誤差基準2で収束判定をする方向が考えられます。
このように反復法では問題の係数行列によって実行時間が大きく変動することになります。
ちなみに10万次元の行列を通常のndarrayでインスタンス化しようとすると74.5 GiBのメモリが必要です。疎行列クラスのありがたみが身にしみます。
Gauss-Seidel法
考え方
$$ % _がmdとして変換されるの回避 \global\def\foi {x^{(k+1)}_i} \global\def\foii {x^{(k+1)}_1} \global\def\foiii{x^{(k+1)}_2} \global\def\foiiii{x^{(k+1)}_{i-1}} $$
ヤコビ法にちょっとした修正を加え、$\foi$を求める際に既に求めた$\foii, \foiii, \dots, \foiiii$を使うようにします。
3元連立一次方程式だと以下のような漸化式になります。
$$ % ここだけ定義 漸化式 \global\def\axrc#1#2{a_{#1#2} x^{(k)}_{#2}} \global\def\axk#1#2{a_{#1#2} \color{red}{x^{(k+1)}_{#2}}\color{black}} \global\def\GSri {x^{(k+1)}_1 = x^{(k)}_1 + \frac{1}{a_{11}}\{b_1 - (\axrc{1}{1}\:\:\: + \axrc{1}{2}\:\:\: + \axrc{1}{3})\}} \global\def\GSrii {x^{(k+1)}_2 = x^{(k)}_2 + \frac{1}{a_{22}}\{b_2 - ( \axk{2}{1} + \axrc{2}{2}\:\:\: + \axrc{2}{3})\}} \global\def\GSriii{x^{(k+1)}_3 = x^{(k)}_3 + \frac{1}{a_{33}}\{b_3 - ( \axk{3}{1} + \axk{3}{2} + \axrc{3}{3})\}} $$ $$ \tag{6} \begin{cases}\begin{aligned}\GSri\\ \GSrii\\ \GSriii \end{aligned}\end{cases} $$
$k$時点よりも解に近いであろう$k+1$の値を用いることで、反復回数の減少が期待できます。
実装に向けて
$(6)$式の右辺を行列として見ると、上三角の項$a_{rc}(r \le c)$では$x^{(k)}$が、下三角の項$a_{rc}(r > c)$では$x^{(k+1)}$が使われていることが分かります。
従って漸化式全体を行列で表すと以下のようになります。
$$ \L = \expL,\;\; \D = \expD,\;\; \U = \expU\;\; とする。\\[4em] \begin{aligned} \xk{k+1} & = \xk{k} + \D^{-1}(\b - (\L\xk{k+1}+ (\D+\U)\xk{k})) \\ & = \D^{-1}\b - \D^{-1}\L\xk{k+1} - \D^{-1}\U\xk{k} \\ \xk{k+1}の項を移項\\[1em] \xk{k+1} + \D^{-1}\L\xk{k+1} & = \D^{-1}\b - \D^{-1}\U\xk{k} \\ (\I + \D^{-1}\L)\xk{k+1} & = \D^{-1}(\b - \U\xk{k}) & (7)\\ \\ \xk{k+1} & = (\I + \D^{-1}\L)^{-1}\D^{-1}(\b - \U\xk{k}) \\ & = (\D + \L)^{-1}(\b - \U\xk{k}) & (8)\\ \end{aligned} $$
$\xk{k+1}$の式が得られました。が、$(8)$式で逆行列が出てきてしまいました。というのも、$(7)$式が$\EQ$という連立一次方程式の形になってしまっています。
実装
連立一次方程式を解くための反復法の中でさらに連立一次方程式を解かなければならないとなると本末転倒なのでが、せっかくなのでちょっとズルをして実装してみたのが下のGS_cheat関数です。
どんなズルかと言うと、反復毎の連立一次方程式を$\I + \D^{-1}\L$が下三角行列なのでscipy.sparse.linalg.spsolve_triangularで解くようにしてみました。
行列としてではなく、素朴に$(6)$式を実装したのが下のGS_simple関数です。
動作確認します。
結論を先取りすると、どちらもscipyの疎行列クラスとしてのうまみを活かせずに生のPythonだとめっちゃ遅いです。
なので100次元の行列で実行してみます。
GS_simpleが遅すぎる…。原理を理解するためとはいえ、jacobiの奮闘を見た後ではゲンナリしてしまいます。
それはさておき反復回数が約半分になったのはいいのですが、なぜかGS_cheatとGS_simpleで1回の差があります。
残差ベクトルのノルムを見てみましょう。
ほとんど重なっていますが44回目あたりで差が出ています。ここらへんは計算方法の違いによる誤差といったところでしょう。
SOR法
考え方
Gauss–Seidel法にさらに修正を加え、修正項に収束を早めるようなパラメータ$\omega$を掛けます。
3元連立一次方程式だと以下のような漸化式になります。
$$ % ここだけ定義 漸化式 % \global\def\axrc#1#2{a_{#1#2} x^{(k)}_{#2}} % \global\def\axk#1#2{a_{#1#2} \color{red}{x^{(k+1)}_{#2}}\color{black}} \global\def\w{\color{red}\omega\color{black}{}} \global\def\SORri {x^{(k+1)}_1 = x^{(k)}_1 + \w\frac{1}{a_{11}}\{b_1 - (\axrc{1}{1}\:\:\: + \axrc{1}{2}\:\:\: + \axrc{1}{3})\}} \global\def\SORrii {x^{(k+1)}_2 = x^{(k)}_2 + \w\frac{1}{a_{22}}\{b_2 - ( \axk{2}{1} + \axrc{2}{2}\:\:\: + \axrc{2}{3})\}} \global\def\SORriii{x^{(k+1)}_3 = x^{(k)}_3 + \w\frac{1}{a_{33}}\{b_3 - ( \axk{3}{1} + \axk{3}{2} + \axrc{3}{3})\}} $$ $$ \tag{9} \begin{cases}\begin{aligned}\SORri\\ \SORrii\\ \SORriii \end{aligned}\end{cases} $$
この係数$\omega$を加速パラメータ(acceleration parameter)というようです3。
例によって行列で表現すると以下のようになります。
$$ \begin{aligned} \xk{k+1} & = \xk{k} + \omega\D^{-1}\Big(\b - \big(\L\xk{k+1}+ (\D+\U)\xk{k}\big)\Big) \\ & = (1-\omega)\xk{k} - \omega\D^{-1}\U\xk{k} + \omega\D^{-1}\b - \omega\D^{-1}\L\xk{k+1} \\ \xk{k+1}の項を移項\\[1em] \big(\I + \omega\D^{-1}\L\big)\xk{k+1} & = \big((1-\omega)\I - \omega\D^{-1}\U\big)\xk{k} + \omega\D^{-1}\b & (10)\\ \end{aligned} $$
実装
素朴に実装した方が遥かに簡単なのですが、GS_simpleがあまりにも遅いのでここでもズルをして$(10)$式を元に実装したのが下のSOR_cheatクラスです。
収束の遅い係数行列(その分次元を減らして20次元)を用意して試してみます。
まずは参考としてjacobi法とGauss-Seidel法でやってみましょう。
本題のSOR法を$0 < \omega < 2$の範囲で$0.1$刻みに設定して実行してみます。
🎨描画関数(plot_resis)定義🎨
$\omega = 1$、すなわちGauss-Seidel法の結果を太実線で示します。
それよりも収束が遅ければ点線、早ければ実線で示します。
SOR法($\omega = 1$)とGauss-Seidel法の推移はほぼ一致していることが分かります(下図)。
今回の係数行列の場合では$\omega = 1.5$で最も収束が早く、Gauss-Seidel法の約1/4となりました。
解くべき問題(係数行列)によって最適なパラメータ$\omega$は異なるので、どうにかして事前に見つけておきたいところですが、それがなかなか難しいようです。4
SOR法はGauss-Seidel法だけが対象か
ところで、SOR法の考え方はjacobi法にも適用できるのではないでしょうか？
学術的・実用的にどれだけ意味があるかは置いておいて、以下のSOR_jacobi関数で確認してみました。
結論を先取りして、収束する範囲の$\omega$で反復した結果が以下です。
上で見た通常のSOR法とは収束するパラメータの範囲や最適な値が違えど、何となく雰囲気は似てるような気がします。
なによりjacobi法で頭打ちになった精度が許容範囲まで減少しているケースもあります($\omega = 0.8$)。
最適なパラメータを探すべく、もうちょっと細かく設定してみます。
この結果では、$\omega = 0.98$で設定すれば精度と収束の速さ両取りできそうです。
通常のSOR法ほど大きな改善とは言えませんが、なかなか面白い結果ではないでしょうか。
ヤコビ法の収束条件
ここまで「収束するとしたら」という表現を使ってきましたが、実際に収束するための条件はなんなのでしょうか。
結論から言うと、係数行列$\A$が対角優位行列(優対角行列、Diagonally dominant matrix)であれば収束するらしいです。
$$ | a_{ii}| > \sum_{i \neq j}| a_{ij}|, \;\;\;\; i = 1, \dots, n $$
これは十分条件で、必要十分条件はもう少し緩めることができるらしいです。
ただし、その必要十分条件を解析的に解くのは不可能なので、実用的には十分条件が使われることが多いらしいです。
また、Gauss-Seidel法の収束条件もヤコビ法と同じらしいです。
ここですべて「らしい」というのは私が数学的な証明を理解できていないからなのですが、以下分からないことの記録として幸谷(2021, pp.126-129)(または幸谷(ウェブ教材))を参照して紹介します。
必要十分条件
ヤコビ法の漸化式を行列で表現した式を以下のように書き直します。
$$\begin{aligned} &\begin{aligned} \xk{k+1} & = \xk{k} + \D^{-1}(\b - \A\xk{k}) & (5)式再掲\\ & = (\I - \D^{-1}\A)\xk{k} + \D^{-1}\b & (11)\\ \end{aligned}\\[1em] &\M = \I - \D^{-1}\A,\;\; \c = \D^{-1}\b\;\; とすると、\\[1em] &\begin{cases} \xk{k+1} & = \M\xk{k}\;\;\, + \c\\ \xk{k} & = \M\xk{k-1} + \c\\ \xk{k-1} & = \M\xk{k-2} + \c\\ & \vdots\\ \xk{1} & = \M\xk{0}\;\;\, + \c\\ \end{cases} &(12)\\[1em] &解\xは(12)式を満たすので\\[1em] &\x = \M\x + \c &(13) \end{aligned}$$
ここで、任意の初期値$\xk{0}$に対して(12)式を用いて得られる列$\{\x_i\}^\infty_0$が収束するためには、行列$\M$のすべての固有値の絶対値が1より小さいことが必要十分条件である。
証明(⇒)
$$\begin{aligned} &k回目の反復によって得られる\xk{k}と\xとの差は\\ &\begin{aligned} \x - \xk{k} & = \M\big(\x - \xk{k-1}\big)\\ & = \M\big(\M\big(\x - \xk{k-2}\big)\big)\\ & \vdots \\ & = \M^k\big(\x - \xk{0}\big) \end{aligned}&(14)\\ &であるから\\ & \| \x - \xk{k} \| \le \| \M \|^k \| \x - \xk{0} \| & (15)\\ &を得る。これより\| \M \| \lt 1 であるためには、固有値はすべて1未満でなくてはならない。 \end{aligned}$$
Arakakiコメント：(14)まではまぁ分かる。(15)もそういうもんだと飲み込むとして、「$\| \M \| \lt 1 であるためには～$」のくだりが🤔❓である。$\| \M \| \lt 1$だったら$k \to \infty$のとき$\xk{k}-\x \to 0$で収束するのは分かるが、行列のノルムの性質や固有値との関係がよく分かっていない。
(⇐)
$$\begin{aligned} &仮定より、\| \M \| \lt 1 であるから\\ &\begin{aligned} \xk{k} & = \M\xk{k-1} + \c \\ & = \M\big(\M\xk{k-2} + \c \big) + \c \\ & \vdots \\ & = \M^k\xk{0} + \big(\sum^{k-1}_{i = 1}\M^i\big)\c & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(16) \end{aligned}\\ &を得る。よって、\lim_{k \to \infty}\xk{k} = \big(\I - \M\big)^{-1}\c となる。\;\;\;\; (17) \end{aligned}$$
Arakakiコメント：上述の繰り返しだが「$仮定より、\| \M \| \lt 1 である$」のがまずわからん。(16)式のシグマは$i=0$では？(17)で$\M^k = 0$として消えているが、$\| \M \| \lt 1$なら$\lim_{k \to \infty}\M^k = 0$となるのは当たり前なの？
十分条件
Gerschgorinの定理を用いて十分条件を得ます。
まず$\M$の成分を書き下すと以下のようになります。
$$ \M = \begin{bmatrix} 0 & -\frac{a_{12}}{a_{11}} & \cdots & -\frac{a_{1n}}{a_{11}}\\ -\frac{a_{21}}{a_{22}} & 0 & \ddots & \vdots\\ \vdots & \ddots & \ddots & -\frac{a_{n-1,n}}{a_{n-1,n-1}} \\ -\frac{a_{n1}}{a_{nn}} & \cdots & -\frac{a_{n,n-1}}{a_{nn}} & 0 \end{bmatrix} $$
$$ \begin{aligned} & Gerschgorinの定理 \\ &\;\;\;\; 正方行列A = [a_{ij}](i,j = 1, 2, \dots,n)に対し、複素平面\mathbb{C}において\\ &\;\;\;\; 中心a_{ii}、半径\sum^n_{j=1,j \ne i}|a_{ij}|の円S_i(i=1,2,\dots,n)を定義した時、\\ &\;\;\;\; 全てのAの固有値\lambdaは\cup^n_{i=1}S_i内に存在する。これは任意のAの固有値\lambdaに対し\\ &\;\;\;\; |a_{ii} - \lambda| \le \sum^n_{j=1,j\ne i}|a_{ij}| \\ &\;\;\;\; を満たすi \in \{1,2,\dots,n\}が必ず存在することを意味する。\\ & \\ & \\ & そうすれば任意の\Mの固有値\lambdaに対して\\ & |\lambda| \le \max_i \sum^n_{j=1, j \ne i}\left|\frac{a_{ij}}{a_{ii}}\right| & (18) \\ & だから \\ & \max_j \sum^n_{i=1, i \ne j}\left|\frac{a_{ij}}{a_{ii}}\right| \le 1 & (19)\\ & が収束のための十分条件となる。これは同時に\| \M \|_1 \le 1 となる条件も満たしている。 \\ & また、\| \cdot \|_\inftyを用いて\|\M\|_\infty \le 1 となる条件を満たすためには、\\ & \max_i \sum^n_{j=1, j\ne i}\left|\frac{a_{ij}}{a_{ii}}\right| \le 1 & (20)\\ & となればよい。 \end{aligned}$$
Arakakiコメント：(19)は「行列$\M$のすべての固有値の絶対値が1より小さい」から？なんで不等号と行列の方向が入れ替わるの…？ここでpノルムを指定しだしたが、無印のノルムはなんだったんだ？やっぱり行列ノルムの理解が欠けてる。
SOR法はヤコビ法で収束しない問題にも使える？
ちょっとした疑問なのですが、WikipediaのSOR法のページで「$\omega$が1より小さいとき ガウス＝ザイデル法より収束が遅くなる。ただし、ガウス＝ザイデル法で収束しないような問題には使える」と書いてありました。5
という訳で、検証してみます。
係数行列として、非ゼロ要素が全て1の3重対角行列を用意します。
ヤコビ法とGauss-Seidel法で実行すると以下のようになります。
Gauss-Seidel法では一瞬収束するかのような素振りを見せましたが、結局発散してしまいました。
続いてSOR法($\omega < 1$)で実行してみます。
むむむ…。これは果たして「使える」と言っていいレベルなのだろうか…
もうちょっと粘って最適値を探してみました。今度は発散に転じると反復を打ち切るようにします。
$\omega = 0.2$で得られた解を見てみましょう。
真の解は1,2,...,20なので、15個目ぐらいまではまぁ…？、いや贔屓目に見てもこの精度では「使える」とは言えないでしょう。
という訳で結論：使えない。
scipyで使える反復法
scipy.sparse.linalgパッケージでは、共役勾配法に代表される非定常反復法(クリロフ部分空間法)の様々なアルゴリズムが利用できます。6
これらは前処理によって大幅に反復回数が改善したり、それぞれが得意とする問題が異なったりするようですが、ここではとりあえず触ってみる程度に留めます。
●全てのアルゴリズムを実行する関数(do_all)定義
🐼描画関数(plot_df)定義🐼
色んな係数行列で実行してみます。
なお収束判定のtolは$\frac{\|\b - \A\xk{k}\|}{\|\b\|} < \varepsilon$という形で使われます(相対的な残差基準)。
今回も例によってspsolve関数を基準に設定してみますが、それ以前に反復回数が上限(maxiter)に達したら反復を打ち切ります(do_all関数内参照)。
🌸収束の遅い係数行列(n = 500, 精度：spsolveの1/10程度)
🌸収束の早い係数行列(n = 100,000, 精度：spsolveの1/2程度、ブラウザだとgcrotmkがめっちゃ遅いので除外)
🌸発散する係数行列(n = 100, 精度：spsolveの1/100程度)
各アルゴリズムの得手不得手が出てそうですが、総じてspsolveが最強では…？
実際のデータを係数行列としてやってみます。
bcsstm22
orani678(gmres, lgmres, gcrotmkが遅いので除外)
lns_3937
memplus(gmresが遅いので除外)
結局、前処理しなかったらspsolveがほぼ最強じゃないですか。
係数行列の性格によって内部でアルゴリズムを使い分けたりしてるのかしら🤔
おまけ：幾何的に理解する
jacobi法を2元連立1次方程式で反復して描画すると色んな形になって面白いです。
🎨描画関数定義🎨
「Run」する度に異なる図になります。
パターン別にサンプルを集めてみました。
皆さんはどんな形が好きですか？
自由欄
参考文献
- 渡部 善隆(1995)「連立1次方程式の基礎知識 : およびGaussの消去法の安定性について」九州大学大型計算機センター広報 Vol.28 No. 4
- (2004)加筆修正版：数値解析チュートリアル2004資料
- 高橋 大輔(1996)『数値計算 理工系の基礎数学8』
- 渡部 善隆(1999)「いつ反復計算をやめるべきか? 収束判定基準の設定方法」九州大学大型計算機センター広報 Vol.32 No. 1
- 幸谷 智紀(2021)『Python数値計算プログラミング』
- 幸谷 智紀(ウェブ教材)「Numerical Computation as Software」：同氏のホームページより。上の書籍の元となったウェブ教材。
-
残差基準：$\|b-Ax^{(k)}\| < \varepsilon$ ならば計算を打ち切る。(渡部, 1999)。
↩ -
誤差基準：$\|x^{(k)} - x^{(k+1)}\| < \varepsilon$ ならば計算を打ち切る。(渡部, 1999)
↩ -
英Wikiでは"relaxation factor"とされていたりするなど、表現に若干の揺れがあります。Successive over-relaxation - Wikipedia
↩ -
最適なパラメータ：数学的な導出方法は、例えば幸谷(ウェブ教材)に詳しいが、Arakakiは解説できるほどの知識がありません…。また、今回参照した範囲では実践的な設定方法などを見つけることはできませんでした。
↩ -
種々の方法の簡単な紹介は例えば渡部(1995)を参照
↩
