行列の簡約化アルゴリズム

簡約な階段行列を機械的に求めるコードを考えました。 処理は以下の通りです。 すべて、行基本変形の操作に対応しています。

1 先頭が0から始まる行の入れ替え

今後の処理において都合が悪いので、行先頭にある0の個数順に行列をソートします。

2 主成分を探す

その行で最初に現れる0でない数を主成分とします。

3 主成分を1にする

主成分が1になるように、行全体を主成分で割ります。

4 主成分の下を0だけにする

y+1列目以降について、matrix[y][pivot_pos]の下がすべて0になるように、y+n行目から、y行目をmatrix[y+n][pivot_pos]倍したものを引きます。

2~4の操作は、上から順に全ての行に適用してください。

5 先頭が0から始まる行の入れ替え

もう一度、行先頭にある0の個数順に行列をソートします。 これによって、行列を階段状にします。

6 主成分を探す

その行で最初に現れる0でない数を主成分とします。

7 主成分を1にする

主成分が1になるように、行全体を主成分で割ります。

8 主成分の上も0だけにする

簡約な階段行列の条件を満たすように、主成分の列の上には何もないようにします。 (下については、手順4で解決しています) y-n行目以前について、y行目をmatrix[y-n][pivot_pos]倍したものを引いてください。

6~8の操作も、上から順に全ての行に適用してください。

ここまでの操作を行うと、元の行列が簡約化されているはずです。 確認してみてください。

サンプルコードは以下の通りです。

def sort_key(array):
    num = 0
    for elm in array:
        if elm == 0:
            num += 1
        else:
            return num
    return num


matrix = [[0, 1, 2], [3, 4, 5], [6, 7, 8]]

y_size = len(matrix)
x_size = len(matrix[0])
pivot_pos = -1

matrix.sort(key=sort_key)

for y in range(y_size):
    for x in range(pivot_pos+1, x_size):
        if matrix[y][x] != 0:
            pivot_pos = x
            break

    pivot_elm = matrix[y][pivot_pos]
    if pivot_elm == 0:
        break

    for x in range(x_size):
        matrix[y][x] /= pivot_elm

    for y_lower in range(y+1, y_size):
        c = matrix[y_lower][pivot_pos]
        for x in range(x_size):
            matrix[y_lower][x] -= c * matrix[y][x]

matrix.sort(key=sort_key)

pivot_pos = -1
for y in range(y_size):
    for x in range(pivot_pos+1, x_size):
        if matrix[y][x] != 0:
            pivot_pos = x
            break

    pivot_elm = matrix[y][pivot_pos]
    if pivot_elm == 0:
        break

    for x in range(x_size):
        matrix[y][x] /= pivot_elm
    
    for y_upper in range(y):
        c = matrix[y_upper][pivot_pos]
        for x in range(x_size):
            matrix[y_upper][x] -= c * matrix[y][x]

print(matrix)