invert_matrix.py

import numpy as np
import copy
from pprint import pprint

# 目的行の列が0だった場合に入れ替える
def swap_row(mtx, rev_mtx, row, col):
    num = mtx[row][col]
    if num == 0.0:
        org_vec = mtx[row]
        org_rev_vec = rev_mtx[row]
        for trg_row in range(row + 1, len(mtx)):
            trg_vec = mtx[trg_row]
            if trg_vec[col] == 0.0:
                continue
            temp_vec = mtx[trg_row]
            temp_rev_vec = rev_mtx[trg_row]
            mtx[row] = temp_vec
            mtx[trg_row] = org_vec
            rev_mtx[row] = temp_rev_vec
            rev_mtx[trg_row] = org_rev_vec
            break
        # 入れ替え先が全部0.0だったらNoneを返す
        else:
            return None, None
    return mtx, rev_mtx

# 掃き出し法
def calc_inv_matrix(mtx):
    row_count = len(mtx)
    org_mtx = copy.copy(mtx)
    # 単位行列で逆行列を初期化
    rev_mtx = []
    for i in range(row_count):
        rev_vec = [.0] * row_count
        rev_vec[i] = 1.0
        rev_mtx.append(rev_vec)
    # 掃き出し開始
    for i in range(row_count):  # 行ループ
        # 行入れ替え
        mtx, rev_mtx = swap_row(mtx, rev_mtx, i, i)
        # 入れ替え先も全部0.0だったらエラー
        if mtx is None:
            print('error :')
            return None
        # 現在の基準行をそれぞれ取得
        base_vec = mtx[i]
        base_rev_vec = rev_mtx[i]
        # 目的行の目的列が1.0になるように逆数を掛ける
        if base_vec[i] != 1.0:
            inv_num = 1.0 / base_vec[i]  # 逆数
            base_vec = list(map(lambda x: x * inv_num, base_vec))
            mtx[i] = base_vec
            base_rev_vec = list(map(lambda x: x * inv_num, rev_mtx[i]))
            rev_mtx[i] = base_rev_vec
        # 他の行の列を0にする
        for j in range(row_count):
            if i == j: # 現在の基準行はスキップ
                continue
            trg_vec = mtx[j]
            trg_rev_vec = rev_mtx[j]
            # すでに0になっていたらスキップ
            if trg_vec[i] == 0:
                continue
            # 0にするための比率を算出
            base_num = base_vec[i]
            trg_num = trg_vec[i]
            ratio = trg_num / base_num
            new_rev_vec = []
            new_vec = []
            # 他の列に対しては同様の比率で引き算する
            for k in range(row_count):
                # 元の行列を定数倍して引き算
                base_num = base_vec[k]
                trg_num = trg_vec[k]
                new_num = trg_num - (base_num * ratio)
                new_vec.append(new_num)
                # 逆行列も同様に定数倍して引き算
                base_rev_num = base_rev_vec[k]
                trg_rev_num = trg_rev_vec[k]
                new_rev_num = trg_rev_num - (base_rev_num * ratio)
                new_rev_vec.append(new_rev_num)
            # 補正後の行ベクトルに置き換え
            rev_mtx[j] = new_rev_vec
            mtx[j] = new_vec
    return rev_mtx

# 適当な行列
mtx = [
      [0, 7, 9, 5],
      [0, 5, 7, 10],
      [1, 3, 4, 8],
      [2, 3, 1, 4],
      ]

# 自前計算
inv_mtx = calc_inv_matrix(mtx[:])
# numpyで答えあわせ
np_inv_mtx = np.linalg.inv(np.array(mtx))

print('元の行列 :')
pprint(mtx, width=20)
print('自前計算の逆行列 :')
pprint(inv_mtx, width=100)

# numpyで直接計算した逆行列
print('numpyの逆行列答え :\n', np_inv_mtx)
# 自前計算の逆行列を元の行列に掛けて単位行列になるか検算
print('numpy検算 :\n', np.dot(inv_mtx, mtx))