Gauss-Jordanの消去法（Pythonで数値計算プログラムを書き直そうシリーズ）

Pythonで数値計算プログラムを書き直そうシリーズ
　間違っていたら教えてください（土下座）

##はじめに
　Gauss-Jordan法は、線形連立方程式（連立1次方程式）の最も基礎的な数値解法です。やっていることは中学で習った加減法ですが、いまいちプログラムの見通しは悪いし、計算速度も非常に遅いです。for文が遅いPythonでプログラムを組むメリットもほぼ無いですが、ポイントなどを確認できればと思い一応書きました。

　ネットで見つけた参考例もいくつか貼っておきます。
C言語でのガウス法による線形連立方程式の求解
ガウス・ジョルダン法

##完成したプログラム

gauss_jordan.py

#Gauss-Jordan法による線型連立方程式の解法プログラム
import numpy as np  #NumPyライブラリ
import matplotlib.pyplot as plt  #データ可視化ライブラリ
import numba  #最適化ライブラリ


#解きたい連立方程式の係数行列
mat_A = np.array([
    [ 1.0, 2.0],
    [ 4.0, 5.0]
])
vec_b = np.array(
    [ 3.0, 6.0]
)


#Gauss-Jordan法（係数行列、係数ベクトル、途中計算表示、微小量）
@numba.jit  #JITで最適化
def gauss_jordan(mat_A, vec_b, process=False, eps=1e-10):
    #拡大係数行列を用意
    vec_b = np.reshape(vec_b, (-1,1))  #列ベクトルに変換
    mat_Ab = np.concatenate((mat_A, vec_b), axis=1)  #横方向に結合
    num_vec_x = mat_Ab.shape[0]  #未知数の数

    #行列の確認
    if (process==True):
        print("pre matrix = ")
        print(mat_Ab)

    #各行でpivot選択して掃き出す
    for pivot in range(num_vec_x):
        #pivot列で成分の絶対値が最大の行を探す
        row_max = pivot
        val_max = mat_Ab[pivot,pivot]
        for row in range(pivot+1, num_vec_x):
            if (val_max < abs(mat_Ab[row,pivot])):
                row_max = row
                val_max = abs(mat_Ab[row,pivot])

        #pivotが小さかったら（係数行列が非正則だったら）例外を発生させる
        if(abs(val_max) < eps):
            print("error: pivot (=", val_max, ") is too small (<=", eps, ").")
            if(abs(mat_Ab[pivot,num_vec_x]) < eps):
                raise Exception("連立方程式は不定です。")
            else:
                raise Excaption("連立方程式は不能（解なし）です。")
            #quit()
        elif(process==True):
            print("pivot = ",val_max)

        #pivotの行と入れ替え
        if (pivot != row_max):
            mat_Ab[pivot,:], mat_Ab[row_max,:] = mat_Ab[row_max,:], mat_Ab[pivot,:].copy()

        #pivot行をpivotで割る(mat_Ab[pivot,pivot]=1にする)
        mat_Ab[pivot, :] = mat_Ab[pivot, :]/val_max

        #掃き出し操作で、mat_Ab[pivot,pivot]より下の係数を0にする
        for row in range(0, num_vec_x):
            #pivot行をmat_Ab[row,pivot]倍してrow行から引く
            if row != pivot:
                pivot_row = mat_Ab[row, pivot]*mat_Ab[pivot, range(pivot, num_vec_x+1)]
                mat_Ab[row, range(pivot, num_vec_x+1)] -= pivot_row

            #val_pivot = mat_Ab[row,pivot]
            #for col in range(pivot, num_vec_x+1):
            #    if (row != pivot):
            #        mat_Ab[row,col] -= val_pivot*mat_Ab[pivot,col]

        if (process==True):
            print(mat_Ab)

    #行列の確認
    if (process==True):
        print("Post matrix = ")
        print(mat_Ab)

    #解の取り出し
    return mat_Ab[:, num_vec_x]


#解の検算
def solution_check(mat_A, vec_b, vec_x):
    print('vec_b = ', vec_b)
    print('Calculated vec_b = ', np.dot(mat_A, vec_x))


#メイン実行部
if (__name__ == '__main__'):
    #Gauss-Jordan法で線型連立方程式の解を計算
    vec_x = gauss_jordan(mat_A, vec_b, process=True)
    print("vec_x = ", vec_x)
    print()

    #解の検算
    solution_check(mat_A, vec_b, vec_x)

　主なプログラムはgauss_jordan関数に書きました。引数として、係数行列$\boldsymbol{A}$、係数ベクトル$\boldsymbol{b}$を渡すと、方程式の解を配列で返します。また、計算途中の結果を表示するかどうか選べるようにしました。あと速度は気にしないつもりでしたが、せっかく関数化しているのでnumbaで高速化しました。「import numba」と「@numba.jit」の2行しか増えてないので大目に見てください。
　個人的にモジュールとして呼び出し可能にしたかったので、そういう構造になっています。モジュールとして使うテストプログラムは以下に書きました。mat_A1、vec_b1の組では解が計算されますが、mat_A2、vec_b2の組ではエラーになり、不定解を持つと表示されます。

gauss_jordan_test.py

import numpy as np  #数値計算用モジュール
import gauss_jordan


mat_A1 = np.array([
    [ 2.0, 1.0, 1.0],
    [ 1.0, 2.0, 1.0],
    [ 1.0, 1.0, 2.0]
])
vec_b1 = np.array(
    [ 7.0, 8.0, 9.0]
)
vec_x1 = gauss_jordan.gauss_jordan(mat_A1,vec_b1)
print("vec_x1 = ", vec_x1)
gauss_jordan.solution_check(mat_A1, vec_b1, vec_x1)  #解の検算
print()


mat_A2 = np.array([
    [ 1.0, 2.0, 3.0],
    [ 5.0, 6.0, 7.0],
    [ 9.0,10.0,11.0]
])
vec_b2 = np.array(
    [ 4.0, 8.0,12.0]
)
vec_x2 = gauss_jordan.gauss_jordan(mat_A2,vec_b2, process=True)
print("vec_x2 = ", vec_x2)
gauss_jordan.solution_check(mat_A2, vec_b2, vec_x2)  #解の検算
print()

##注意すべきと思う点

###係数行列と係数ベクトルを渡すと未知数ベクトルを返す関数を作る
　最初に書くことではないかもしれないですが、書籍などで見るプログラムは汎用性のないものが多すぎるように思います。よくあるのは、プログラムの冒頭でa[][]のような配列を拡大係数行列として宣言し、その配列を直接操作しているプログラムです。これだと他の係数で計算させたい時に不便です。また酷い本では、a[0][0]=1,a[0][1]=2,…のように成分を1個ずつ代入して行列を作ってるような例もありました。プログラムが冗長になるし、可読性も悪いです。
　今回のプログラムでは、scipy.linalgのように係数行列mat_Aと係数ベクトルvec_bを渡すと、未知数ベクトルvec_xの配列が返ってくる仕様にしました。Pythonだと行列の次数も行の長さから取得できて便利です。今回のプログラム中では結局、np.concatenateを使ってmat_Abという拡大係数行列を作っていますが、引数としては係数行列と係数ベクトルを別々に渡せる方が便利でしょう。

###pivot選択をして桁落ちなどを防ぐ
　係数行列が非正則行列（$rank \boldsymbol{A} < N$）の場合、どこかの対角要素が0になります。そのため何の工夫もないと、1にしたい対角要素が0だった時にゼロ除算が発生し、エラーで止まってしまいます。また、小さい数で割り算をすると桁落ちが発生します。そうしたことを避けるため、pivot選択は必ずすべきです。今回のプログラムではpivot選択をした上で、pivotがほぼ0だったら計算を中止するようにしていますが、加えて検算をできる関数も作りました（行列の掛け算をしてるだけですが）。

###求解できない時は不定か不能か分かるようにする
　上の項目に関連しますが、係数行列が非正則行列だと解を求められません。しかし、係数行列が非正則の連立方程式は常に解を持たないかというと、そうとも限りません。
　例えば、次の連立方程式は解を持ちません（不能と言います）。

\begin{bmatrix}
  1& 1\\
  2& 2
\end{bmatrix}
\begin{Bmatrix}
  x_1\\ x_2
\end{Bmatrix}
 = \begin{Bmatrix}
  1\\ 1
\end{Bmatrix}

　実際、係数行列の行列式は$det \boldsymbol{A} = 0$となり、非正則であることを確認できます。より直感的な理解としては、上の2直線をグラフにするといいでしょう。$x_1 +x_2 = 1$と$2x_1 +2x_2 = 1$は、互いに平行な直線です。なので交点を持たず、解も存在しないことが理解できます。
　一方、次の例は係数ベクトル$\boldsymbol{b}$を変えただけですが、この連立方程式は一意でない解を持ちます（不定と言います）。

\begin{bmatrix}
  1& 1\\
  2& 2
\end{bmatrix}
\begin{Bmatrix}
  x_1\\ x_2
\end{Bmatrix}
 = \begin{Bmatrix}
  1\\ 2
\end{Bmatrix}

　こちらの行列式も$det \boldsymbol{A} = 0$で非正則ですが、式を見ると$x_1 +x_2 = 1$の両辺を2倍したものが$2x_1 +2x_2 = 2$になっていることが分かります。つまりこの2式は同じ式であり、グラフに描けばこの2直線は一致します。そのため解は、「$x_1 +x_2 = 1$上の全ての点」ということになります。
　上のプログラムの例では、mat_A2、vec_b2の組において計算できませんでしたが、これも不定解をもつ係数の組です。例えば$x_1=-1, x_2=1, x_3=1$を代入すると、連立方程式の解になっていることが分かります。今回のプログラムでは、係数ベクトルが0かどうかをみて、不定か不能か判定を出すようにしました。

###計算速度は非常に遅い（特にPython）
　今回のアルゴリズムは特に何の工夫もしていないため、計算速度は他の解法と比較しても一番遅いと思われます（$O(n^3)$くらい？）。また、今回はとくにfor文が遅いと言われるPythonを使っているため、余計に遅くなっています。線形連立方程式を解くとなると、数値計算の分野ではFortranが未だに現役です。Fortranでなければ、せめてCやJavaとかを使いたいところです（今ならJulia？）。

##自作プログラムは使えるのか
　今回に関しては「使えない」と言い切っていいと思います。場合にもよりますが、まともに連立1次方程式を解きたければ、まずはLAPACKとかのライブラリを使うべきでしょう。今回のプログラムでは計算速度が遅いし、性能も不十分です。
　連立1次方程式の数値解法（というか線形代数）は深い沼の世界で、その沼を潜るにはかなり抽象的な数学を必要とします。自分もその淵にだけ立ったことがありますが、ひとまずはライブラリを大人しく使うのが得策でしょう。ちなみに、PythonであればAnacondaを使う人も多いと思いますが、現状では内部にMKLが入っています。そしてMKLはBLASやLAPACKを利用しており、これらはFortran90で書かれています。NumPyの内部環境は以下のコードで調べられるので、一度確認しておくと良いと思います。

import numpy
print(numpy.show_config())

　最後に、SciPyでの計算例を載せておきます。これを実行してみると、mat_A2、vec_b2の組で「Ill-conditioned matrix detected. （条件の悪い行列が検出されました）」というメッセージの後、vex_x2 = [-0.28571429, -0.42857143, 1.71428571]と出ましたが、検算したら特殊解として合ってそうでした。
　Gauss-Jordan法は連立1次方程式の数値解法の中でも、「直接法」と呼ばれる解法に分類されます。一方でSciPyでは、主に計算速度の点から「反復法」に分類される手法を採用しているはずです。反復法であれば、不定解をもつ係数の組でも特殊解は出してくれる、ということなんでしょう。

import numpy as np  #数値計算用モジュール
import scipy.linalg  #SciPyの線形計算ライブラリ
import gauss_jordan


mat_A1 = np.array([
    [ 2.0, 1.0, 1.0],
    [ 1.0, 2.0, 1.0],
    [ 1.0, 1.0, 2.0]
])
vec_b1 = np.array(
    [ 7.0, 8.0, 9.0]
)
vec_x1 = scipy.linalg.solve(mat_A1,vec_b1)
print("vec_x1 = ", vec_x1)
gauss_jordan.solution_check(mat_A1, vec_b1, vec_x1)  #解の検算
print()


mat_A2 = np.array([
    [ 1.0, 2.0, 3.0],
    [ 5.0, 6.0, 7.0],
    [ 9.0,10.0,11.0]
])
vec_b2 = np.array(
    [ 4.0, 8.0,12.0]
)
vec_x2 = scipy.linalg.solve(mat_A2,vec_b2)
print("vec_x2 = ", vec_x2)
gauss_jordan.solution_check(mat_A2, vec_b2, vec_x2)  #解の検算
print()