LoginSignup
4
1

More than 5 years have passed since last update.

@MFFMMFFMM(MAKOTO FUKUDA)

逆行列を求めるコード

Posted at

『入門Python3』(Bill Lubanovic　著、斎藤 康毅　監訳、長尾 高弘　訳)を半分くらい読んだので試しに何かコードを書いてみようと思い, 逆行列をもとめるコードを書くことにしました.

はじめに

numpy使えよとか無駄に長いよとかそういう意見もあると思いますがこれはあくまで python の練習問題として自分で設定したものです. 初めて記事を書くので至らない点もあると思います. それを承知の上で付き合っていただけると嬉しいです. 本文は文語調で書きます.

環境

  • Windows10
  • Spyder 3.2.3

数学的な基礎

二通りの求め方が広く知られている.

行基本変形を使う方法

$n$次正方行列$A$に対してn次の単位行列Eを右に置いた$n×2n$行列
$$
\begin{bmatrix}
A & E
\end{bmatrix}
$$
を考え,この行列の区分け$A$が単位行列になるまで行基本変形し,つぎの行列を得る.
$$
\begin{bmatrix}
E & A'
\end{bmatrix}
$$
この時の上の行列の区分け$A'$がの行列$A$の逆行列$A^{-1}$である.

余因子行列を使う方法

$n$次正方行列$A$に逆行列が存在するとき, $A$の行列式を$|A|$, 余因子行列を$\tilde{A}$とするとAの逆行列は次のように表せる.
$$
A^{-1}=\frac{1}{|A|}\tilde{A}
$$
今回は余因子行列を用いた方法で逆行列を求める.

行列式を求める

余因子行列を求めるには行列式が必要だ.　行基本変形と余因子展開を組み合わせて行列式を求める.　$a_{11}\not=0$のとき
$$
{\left|
\begin{array}{cccc}
\begin{bmatrix}
a_{11} & 0 & \ldots & 0 \\
a_{21} & a_{22} & \ldots & a_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{m1} & a_{m2} & \ldots & a_{mn}
\end{bmatrix}
\end{array}
\right|}
=a_{11}×
{\left|
\begin{array}{cccc}
\begin{bmatrix}
a_{22} & \ldots & a_{2n} \\
\vdots & \ddots & \vdots \\
a_{m2} & \ldots & a_{mn}
\end{bmatrix}
\end{array}
\right|}
$$
という関係を用いて行列をどんどん小さくしていき, 最終的に変形すべき行列が$1×1$行列になったら今まで外に出してきた$(1,1)$成分をすべて掛け合わせると行列式が求まる.　$(1,1)$成分の値を取り出すのはとりあえず置いておいて, 与えられた$n$次正方行列の第1行が$(1,0,\ldots\ldots,0)$となるように行基本変形して新たな$n-1$次正方行列を返す関数を求める.引数は行列, つまりリストのリストだ.

例えば
$$A=\begin{bmatrix}4&3&2&1\\
1&4&3&2\\
2&1&4&3\\
3&2&1&4\end{bmatrix}$$
という行列に対して

\begin{bmatrix}
4&3&2&1\\
1&4&3&2\\
2&1&4&3\\
3&2&1&4
\end{bmatrix}
\rightarrow
\begin{bmatrix}
4&3-4\cdot\frac{3}{4}&2-4\cdot\frac{2}{4}&1-4\cdot\frac{1}{4}\\
1&4-1\cdot\frac{3}{4}&3-1\cdot\frac{2}{4}&2-1\cdot\frac{1}{4}\\
2&1-2\cdot\frac{3}{4}&4-2\cdot\frac{2}{4}&3-2\cdot\frac{1}{4}\\
3&2-3\cdot\frac{3}{4}&1-3\cdot\frac{2}{4}&4-3\cdot\frac{1}{4}
\end{bmatrix}
\rightarrow
\begin{bmatrix}
4&0&0&0\\
1&\frac{13}{4}&\frac{5}{2}&\frac{7}{4}\\
2&-\frac{1}{2}&3&\frac{5}{2}\\
3&-\frac{1}{4}&-\frac{1}{2}&\frac{13}{4}
\end{bmatrix}

というように変形して

\begin{bmatrix}
\frac{13}{4}&\frac{5}{2}&\frac{7}{4}\\
-\frac{1}{2}&3&\frac{5}{2}\\
-\frac{1}{4}&-\frac{1}{2}&\frac{13}{4}
\end{bmatrix}

という行列を取り出す.

def reduct(matrix):
    l = len(matrix)
    w = [[0 for j in range(len(matrix[0])-1)] for k in range(len(matrix)-1)]
    for row in range(1,len(matrix[0])):
        for column in range(1,len(matrix)):
            s = matrix[column][row]-(matrix[column][0]/matrix[0][0])*matrix[0][row]
            w[column-1][row-1] = s
    return w

>>>A=[[4,1,2,3],[3,4,1,2],[2,3,4,1],[1,2,3,4]]
>>>print(reduct(A))
[[3.25, -0.5, -0.25], [2.5, 3.0, -0.5], [1.75, 2.5, 3.25]]

これでやっと行列式を求める準備が整った.　以下のコードでは u という, 与えられた行列の列数に長さが等しいリストを用意し, reduct()するたびに$(1,1)$成分を保存している.　最後に u の各成分を掛け合わせれば行列式が求められる.

def determinant(matrix):
    l = len(matrix)  
    u = [0 for r in range(0,l)]
    d = 1
    def reduct(matrix):
        w = [[0 for j in range(len(matrix[0])-1)] for k in range(len(matrix)-1)]
        for row in range(1,len(matrix[0])):
            for column in range(1,len(matrix)):
                s = matrix[column][row]-(matrix[column][0]/matrix[0][0])*matrix[0][row]
                w[column-1][row-1] = s
        return w


    def loop(matrix):
        if len(matrix) > 1:
            u[-len(matrix)] = matrix[0][0]
            matrix = reduct(matrix)
        elif len(matrix) == 1:
            u[-1]= matrix[0][0]
        return matrix


    while l > 0:
        matrix = loop(matrix)
        l = l-1


    for ele in range(0,len(u)):
        d = d*u[ele]

    return d

>>> A=[[1,2,3,4],[4,1,2,3],[3,4,1,2],[2,3,4,1]]
>>> print("|A|=",determinant(A))
|A|= -160.0

行列Aに逆行列が存在しない場合, 0で割る処理が含まれてエラーになる.

>>> A=[[1,2,3,4],[1,2,3,4],[3,4,1,2],[2,3,4,1]]
>>> print("|A|=",determinant(A))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 26, in determinant
  File "<stdin>", line 17, in loop
  File "<stdin>", line 9, in reduct
ZeroDivisionError: float division by zero

余因子行列を求める

$A$の余因子行列の第$(j,i)$成分は$A$の第$(i,j)$余因子である.　$A$の第$(i,j)$余因子は$A$の第$(i,j)$小行列の行列式に$(-1)^{i+j}$をかけたものである. 以下に示すのは与えられた行列の第$(i,j)$小行列を返す関数である. elseはなくてもよい.

def coef(matrix,row,column):
    u = [[0 for k in range(len(matrix[0])-1)] for l in range(len(matrix)-1)]
    for j in range(0,len(matrix[0])):
        for i in range(0,len(matrix)):
            if   (j < (row-1)) and (i < (column-1)):
                u[i][j] = matrix[i][j]
            elif (j < (row-1)) and (i > (column-1)):
                u[i-1][j] = matrix[i][j]
            elif (j > (row-1)) and (i < (column-1)):
                u[i][j-1] = matrix[i][j]
            elif (j > (row-1)) and (i > (column-1)):
                u[i-1][j-1] = matrix[i][j]
            else :
                pass
    return u

逆行列を求める

第$(i,j)$小行列が分かれば余因子はすぐわかる. あとはこれを新しい$n$次正方行列 inv の第$(j,i)$成分に代入していけばよい. ここでは逆行列そのものが欲しいので各成分を$|A|$で割った値を代入する. 以上のことをまとめると以下のようなコードになる.

def inverse(matrix):
    def determinant(matrix):
        l = len(matrix)
        u = [0 for r in range(0,l)]  
        d = 1
        def reduct(matrix):
            w = [[0 for j in range(len(matrix[0])-1)] for k in range(len(matrix)-1)]
            for row in range(1,len(matrix[0])):
                for column in range(1,len(matrix)):
                    s = matrix[column][row]-(matrix[column][0]/matrix[0][0])*matrix[0][row]
                    w[column-1][row-1] = s
            return w

        def loop(matrix):
            if len(matrix) > 1:
                u[-len(matrix)] = matrix[0][0]
                matrix = reduct(matrix)
            elif len(matrix) == 1:
                u[-1]= matrix[0][0]
            else:
                print("error")
            return matrix

        while l > 0:
            matrix = loop(matrix)
            l = l-1

        for ele in range(0,len(u)):
            d = d*u[ele]
        return d

    def coef(matrix,row,column):
        u = [[0 for k in range(len(matrix[0])-1)] for l in range(len(matrix)-1)]
        for j in range(0,len(matrix[0])):
            for i in range(0,len(matrix)):
                if   (j < (row-1)) and (i < (column-1)):
                    u[i][j] = matrix[i][j]
                elif (j < (row-1)) and (i > (column-1)):
                    u[i-1][j] = matrix[i][j]
                elif (j > (row-1)) and (i < (column-1)):
                    u[i][j-1] = matrix[i][j]
                elif (j > (row-1)) and (i > (column-1)):
                    u[i-1][j-1] = matrix[i][j]
                else :
                    pass
        return u

    inv=[[0 for i in range(0,len(matrix[0]))]for j in range(0,len(matrix))]   

    for column in range(1,len(inv[0])+1):
        for row in range(1,len(inv)+1):
            inv[row-1][column-1] = ((-1)**(row+column))*determinant(coef(matrix,row,column))/determinant(matrix)

    return inv

試しに何か入れてみよう.

>>>A=[[1,2,3,4],[4,1,2,3],[3,4,1,2],[2,3,4,1]]
>>>print("A^(-1)=",inverse(A))
A^(-1)= [[-0.225, 0.275, 0.025, 0.025], [0.025, -0.225, 0.275, 0.025000000000000022], [0.025, 0.025, -0.225, 0.275], [0.275, 0.025, 0.025000000000000012, -0.225]]

これが$A$の逆行列かどうか確かめるために新たに行列の積を返す関数を定義する.

def matrix_product(matrix_1,matrix_2):
    t = 0
    w = [[0 for i in range(len(matrix_1[0]))] for j in range(len(matrix_2))]
    for column in range(0,len(matrix_2)):
        for row in range(0,len(matrix_1[0])):
            for i in range(0,len(matrix_2[column])):
                s = matrix_1[i][row]*matrix_2[column][i]
                t += s
            w[column][row] = t
            t = 0 
    return w

>>>print("E=",matrix_product(A,inverse(A)))
E= [[1.0000000000000002, 2.7755575615628914e-17, 0.0, 4.85722573273506e-17], [1.1102230246251565e-16, 1.0, 1.1102230246251565e-16, 0.0], [0.0, 1.1102230246251565e-16, 1.0, 0.0], [5.551115123125783e-17, 1.1102230246251565e-16, 1.1102230246251565e-16, 1.0]]

これで検算もできた.

おわりに

初めて長めのコードを書いたのでコードを書く上でのマナーは知りません.タグとかあってるかよくわからないです. 批判があれば喜んで受け入れるので遠慮なく書いてほしいです.

4
1
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
1