Gaussian Eliminationアルゴリズムの実装

NumpyやScipyなどのモジュールを使用すると、RREFを取得することはできますが、REFを取得することができません。そのため、自力でREFを取得するアルゴリズムについて考えます。

配列の扱いが面倒なので、今回はPythonを使用します。

Helper functions

配列の先頭の０要素をカウントする関数

def get_zeros(row:np.ndarray) -> int:
    for i in range(len(row)):
        if row[i] != 0:
            return i
    return len(row)

行の基本変形(ERO)を行う関数

この関数がメインの関数になります。

def ERO(mat:np.ndarray) -> np.ndarray:
    A = [[0 for i in mat[0]] for j in mat]
    pivot = mat[0]
    print(f"pivot:[{','.join([str(i) for i in pivot])}]")
    A[0] = pivot
    for i in range(len(mat)):
        if i != 0:
            for j in range(len(mat[i])):
                toSub = pivot[j]*mat[i][0]/pivot[0]
                A[i][j] = mat[i][j]-toSub
    return np.array(A)

列の削減関数

getzerosの返り値に応じて行列を加工する関数。

def reduce_operation(mat:np.ndarray) -> np.ndarray:
    for n in range(1,len(mat)):
        if get_zeros(mat[n]) > 0:
            return np.array([[i for i in j[n:]] for j in mat[n:]])
    return mat

辞書のソート関数

スワップ回数を記録するための、辞書ソートを拡張した関数。

def sort_d(d:dict):
    items = list(d.items())
    n_swapped = 0
    swapped = True
    while swapped:
        swapped = False
        n_swapped+=1
        for i in range(len(items)-1):
            if items[i][1] > items[i+1][1]:
                tmp = items[i]
                items[i] = items[i+1]
                items[i+1] = tmp
                swapped = True
    return items , n_swapped

列のソート関数

def sort_array(mat:np.ndarray) -> np.ndarray:
    try:
        zeros = {i:get_zeros(n) for i,n in enumerate(mat)}
        return mat[[i[0] for i in sorted(zeros.items(), key=lambda i:i[1])]]
    except:
        return mat

Gaussian elimination

上記の関数を組み合わせて、gausiaan eliminationを行う関数を作成します。

例として以下の行列について考えます。

$A =\begin {pmatrix}1&-2&1&4&-5\\1&1&-2&3&-3\\2&-1&-1&2&2\\5&-1&0&5&5\\2&2&0&4&-1\end{pmatrix}$


A = np.array(
[[1,-2,1,4,-5],
[1,1,-2,3,-3],
[2,-1,-1,2,2],
[5,-1,0,5,5],
[2,2,0,4,-1]])

print(A)

def gaussian_eliminatioin(mat:np.ndarray) -> np.ndarray:
    
    A = ERO(mat)
    REF = [list(A[0])]
    for i in range(len(mat[0])-1):
        A = reduce_operation(A)
        A = ERO(A)
        A = sort_array(A)
        REF.append([0 for k in range(i+1)] + list(A[0]))
    return np.array(REF)

print(gaussian_eliminatioin(A))

>>>
[[ 1 -2  1  4 -5]
 [ 1  1 -2  3 -3]
 [ 2 -1 -1  2  2]
 [ 5 -1  0  5  5]
 [ 2  2  0  4 -1]]
pivot:[1,-2,1,4,-5]
pivot:[3.0,-3.0,-1.0,2.0]
pivot:[4.0,-12.0,24.0]
pivot:[10.0,-19.0]
pivot:[0.5]
[[  1.   -2.    1.    4.   -5. ]
 [  0.    3.   -3.   -1.    2. ]
 [  0.    0.    4.  -12.   24. ]
 [  0.    0.    0.   10.  -19. ]
 [  0.    0.    0.    0.    0.5]]

上記から、

$G(A)=\begin{pmatrix}1.0&-2.0&1.0&4.0&-5.0\\0.0&3.0&-3.0&-1.0&2.0\\0.0&0.0&4.0&-12.0&24.0\\0.0&0.0&0.0&10.0&-19.0\\0.0&0.0&0.0&0.0&0.5\end{pmatrix}$

小数から分数への変換

上記のndarrayではfloat型の値が並んでいるため、数学的には見栄えが良くありません。そこで、小数と整数を区別して自動的に行列を再構築する方法を考えます。

float→fracについては、以下の方法で変更することができそうです。

from decimal import Decimal
A = Decimal('0.5').as_integer_ratio()
>>>
(0,5)

上記を利用して、行列を再構築します。

from decimal import Decimal
from fractions import Fraction as frac 

A = [[str(frac(*Decimal(f'{i}').as_integer_ratio())) for i in j] for j in A]
[['1', '-2', '1', '4', '-5'], ['0', '3', '-3', '-1', '2'], ['0', '0', '4', '-12', '24'], ['0', '0', '0', '10', '-19'], ['0', '0', '0', '0', '1/2']]

$\begin {pmatrix}1&-2&1&4&-5\\0&3&-3&-1&2\\0&0&4&-12&24\\0&0&0&10&-19\\0&0&0&0&\frac{1}{2}\end{pmatrix}$