物理のノート

逆行列

逆行列の数値計算($2\times2$正方行列)

逆行列$A^{-1}$は,行列式と余因子(Cofactor)を用いて,次のように定義される.
次の$2\times2$正方行列の場合, \begin{equation} A= \begin{pmatrix} a & b\\ c & d \end{pmatrix} , \end{equation} 逆行列$A^{-1}$は次のようになる.

\begin{equation} A^{-1} =\frac{1}{\lvert A \rvert } \begin{pmatrix} C_{1,1} & C_{2,1}\\ C_{1,2} & C_{2,2} \end{pmatrix} =\frac{1}{\lvert A \rvert } \begin{pmatrix} (-1)^{1+1}\Delta_{1,1} & (-1)^{2+1}\Delta_{2,1}\\ (-1)^{1+2}\Delta_{1,2} & (-1)^{2+2}\Delta_{2,2} \end{pmatrix} = \frac{1}{ad-bc } \begin{pmatrix} d & (-1)b\\ (-1)c & a \end{pmatrix}\\ = \frac{1}{ad-bc } \begin{pmatrix} d & -b\\ -c & a \end{pmatrix} \end{equation} この定義が示すように,逆行列が存在するには, 行列式$\lvert A \rvert \ne 0$が条件となる.

定義に従って,次の$2\times 2$正方行列の逆行列を求める. \begin{equation} A= \begin{pmatrix} 1 & 3\\ 2 & 4 \end{pmatrix} \end{equation} とすると, \begin{equation} A^{-1}=\frac{1}{1\cdot 4-3\cdot 2} \begin{pmatrix} 4 & (-1)3\\ (-1)2 & 1 \end{pmatrix}\\ =-\frac{1}{2} \begin{pmatrix} 4 & -3\\ -2 & 1 \end{pmatrix} = \begin{pmatrix} -2 & \frac{3}{2}\\ 1 & -\frac{1}{2} \end{pmatrix} \end{equation}
ところで,この$2\times 2$正方行列の逆行列は,Gauss-Jordan消去法(掃き出し法)で用いた数値計算を次のように改良して求めることができる.


=============
2 * 4 array
=============
array2 =
[[ 1.  3.  1.  0.]
 [ 2.  4.  0.  1.]]
=====================================
Gauss-Jordan Elimination(掃き出し法)
=====================================
array2 =
[[ 2.  2.  0.  1.]
 [ 1.  3.  1.  0.]]
array2 =
[[ 2.  2.  0.  1.]
 [ 1.  3.  1.  0.]]
array2 =
[[ 2.   2.   0.   0.5]
 [ 1.   3.   1.   0. ]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   3.   1.   0. ]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   1.   1.   0. ]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   1.   1.   0. ]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   1.   1.  -0.5]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   1.   1.  -0.5]]
array2 =
[[ 1.   2.   0.   0.5]
 [ 0.   1.   1.  -0.5]]
array2 =
[[ 1.   0.   0.   0.5]
 [ 0.   1.   1.  -0.5]]
array2 =
[[ 1.   0.  -2.   0.5]
 [ 0.   1.   1.  -0.5]]
array2 =
[[ 1.   0.  -2.   1.5]
 [ 0.   1.   1.  -0.5]]
=====================
determinant of matrix
=====================
det = 2.0

Pythonで実装したアルゴリズムは次の通りである.

#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('2 * 2 array')
print('=============')
array2 = np.array([[1.,3.,1.,0.],[2.,4.,0.,1.]])
arrayp = np.array([[0.,0.],[0.,0.]])
arrayq = np.array([[0.,0.],[0.,0.]])
tmp = np.array([[0.,0.,0.,0.],[0.,0.,0.,0.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination(掃き出し法)')
print('=====================================')
#n = 2  元の数
det = 1.0  
for k in range(0,2,1): 
    for i in range(k,2,1):
        if i !=k:
            if abs(array2[k,k]) < abs(array2[i,k]):
    # j列の範囲をk列目から3列目に変更する.                
                for j in range(k,4,1):
                    tmp[i,j] = array2[k,j]
                    array2[k,j] = array2[i,j]
                    array2[i,j] = tmp[i,j]           

    arrayp[k,k] = array2[k,k]
    det = det * arrayp[k,k]
    # j列の範囲をk列目から3列目に変更する.
    for j in range(k,3,1):
        array2[k,j+1] = array2[k,j+1] / arrayp[k,k]
        print('array2 =') 
        print(array2)     
    array2[k,k] = 1

    for i in range(0,2,1):
        arrayq[k,k] = array2[i,k]
        # j列の範囲をk列目から3列目に変更する.
        for j in range(k,4,1):
            if i != k:
                array2[i,j] = array2[i,j] - array2[k,j] * arrayq[k,k]
                print('array2 =') 
                print(array2)           
            else:
                array2[i,j] = array2[i,j] 
                
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det)) 

このアルゴリズムでは,行列$A$の列を右に拡張し,$2\times2$単位行列を加えている.
この状態で,掃き出し法(Gauss-Jordan消去法)を実行すると,実行結果には$2\times2$単位行列と逆行列$A^{-1}$が並んだ$2\times4$行列が前述のように得られる.

念のため,検算をしてみよう.
ある行列とその逆行列の積は単位行列になる. \[ AA^{-1}=I \] 行列同士の積の定義に沿ってコードを実装する.

#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('==========================')
print('2 * 2 array multiplication')
print('==========================')
a = np.array([[1.,3.],[2.,4.]])
b = np.array([[-2.,1.5],[1.,-0.5]])
c = np.array([[0.,0.],[0.,0.]])
#行列Aをaとする.
print('a =')
print(a)
#行列Aの逆行列をbとする.
print('b =')
print(b)

for i in range(0,2,1): 
    for j in range(0,2,1):
            for k in range(0,2,1):
                c[i,j] += a[i,k]*b[k,j]
        
#行列aと行列bの乗算結果を行列cとする.
print('c =')
print(c)

実行結果は次の通り.

==========================
2 * 2 array multiplication
==========================
a =
[[ 1.  3.]
 [ 2.  4.]]
b =
[[-2.   1.5]
 [ 1.  -0.5]]
c =
[[ 1.  0.]
 [ 0.  1.]]

単位行列が得られた.

逆行列の数値計算($3\times 3$正方行列)

$3\times3$正方行列についても同様に,逆行列を求めてみる.
行列式の数値計算で利用した次の行列を用いる. \begin{equation} \begin{pmatrix} 3 & 1 & 2 \\ 5 & 1 & 3 \\ 4 & 2 & 1 \\ \end{pmatrix} \end{equation} 結果は次のようになる.

[[ 1.    0.    0.   -1.25  0.75  0.25]
 [ 0.    1.    0.    1.75 -1.25  0.25]
 [ 0.    0.    1.    1.5  -0.5  -0.5 ]]

途中の計算プロセスと行列式を記す.

=============
3 * 4 array
=============
array2 =
[[ 3.  1.  2.  1.  0.  0.]
 [ 5.  1.  3.  0.  1.  0.]
 [ 4.  2.  1.  0.  0.  1.]]
=====================================
Gauss-Jordan Elimination(掃き出し法)
=====================================
array2 =
[[ 5.   0.2  3.   0.   1.   0. ]
 [ 3.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 5.   0.2  0.6  0.   1.   0. ]
 [ 3.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 5.   0.2  0.6  0.   1.   0. ]
 [ 3.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 5.   0.2  0.6  0.   0.2  0. ]
 [ 3.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 5.   0.2  0.6  0.   0.2  0. ]
 [ 3.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   1.   2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  2.   1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.   0.   0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 4.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   2.   1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   1.2  1.   0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   1.2 -1.4  0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   1.2 -1.4  0.   0.   1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   1.2 -1.4  0.  -0.8  1. ]]
array2 =
[[ 1.   0.2  0.6  0.   0.2  0. ]
 [ 0.   0.4  0.2  1.  -0.6  0. ]
 [ 0.   1.2 -1.4  0.  -0.8  1. ]]
array2 =
[[ 1.          0.2         0.6         0.          0.2         0.        ]
 [ 0.          1.2        -1.16666667  0.         -0.8         1.        ]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.2         0.6         0.          0.2         0.        ]
 [ 0.          1.2        -1.16666667  0.         -0.8         1.        ]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.2         0.6         0.          0.2         0.        ]
 [ 0.          1.2        -1.16666667  0.         -0.66666667  1.        ]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.2         0.6         0.          0.2         0.        ]
 [ 0.          1.2        -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.6         0.          0.2         0.        ]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.2         0.        ]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.2         0.        ]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333  0.        ]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.4         0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.2         1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.         -0.6         0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.         -0.33333333  0.        ]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.         -0.33333333 -0.33333333]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.5        -0.33333333 -0.33333333]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.5        -0.5        -0.33333333]]
array2 =
[[ 1.          0.          0.83333333  0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          0.66666667  1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.          0.          0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.33333333 -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.75       -0.16666667]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.75        0.25      ]
 [ 0.          1.         -1.16666667  0.         -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.75        0.25      ]
 [ 0.          1.          0.          0.         -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.75        0.25      ]
 [ 0.          1.          0.          1.75       -0.66666667  0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.          0.          0.         -1.25        0.75        0.25      ]
 [ 0.          1.          0.          1.75       -1.25        0.83333333]
 [ 0.          0.          1.          1.5        -0.5        -0.5       ]]
array2 =
[[ 1.    0.    0.   -1.25  0.75  0.25]
 [ 0.    1.    0.    1.75 -1.25  0.25]
 [ 0.    0.    1.    1.5  -0.5  -0.5 ]]
=====================
determinant of matrix
=====================
det = 4.0

Pythonコードを以下に記す.

#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('=============')
print('3 * 6 array')
print('=============')
array2 = np.array([[3.,1.,2.,1.,0.,0.],[5.,1.,3.,0.,1.,0.],[4.,2.,1.,0.,0.,1.]])
arrayp = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
arrayq = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
tmp = np.array([[0.,0.,0.,0.,0.,0.],[0.,0.,0.,0.,0.,0.],[0.,0.,0.,0.,0.,0.]])
print('array2 =')
print(array2)
print('=====================================')
print('Gauss-Jordan Elimination(掃き出し法)')
print('=====================================')
#n = 3  ##元の数はこの場合,3
det = 1.0  
for k in range(0,3,1): 
    for i in range(k,3,1):
        if i !=k:
            if abs(array2[k,k]) < abs(array2[i,k]):
    # j列の範囲をk列目から5列目に変更する.                
                for j in range(k,6,1):
                    tmp[i,j] = array2[k,j]
                    array2[k,j] = array2[i,j]
                    array2[i,j] = tmp[i,j]           

    arrayp[k,k] = array2[k,k]
    det = det * arrayp[k,k]
    # j列の範囲をk列目から5列目に変更する.
    for j in range(k,5,1):
        array2[k,j+1] = array2[k,j+1] / arrayp[k,k]
        print('array2 =') 
        print(array2)     
    array2[k,k] = 1

    for i in range(0,3,1):
        arrayq[k,k] = array2[i,k]
        # j列の範囲をk列目から5列目に変更する.
        for j in range(k,6,1):
            if i != k:
                array2[i,j] = array2[i,j] - array2[k,j] * arrayq[k,k]
                print('array2 =') 
                print(array2)           
            else:
                array2[i,j] = array2[i,j] 
                
##行列式の計算結果
print('=====================')
print('determinant of matrix')
print('=====================')
print('det = {}'.format(det))  

行列とその逆行列の積が単位行列になるか確かめてみる.

==========================
3 * 3 array multiplication
==========================
a =
[[ 3.  1.  2.]
 [ 5.  1.  3.]
 [ 4.  2.  1.]]
b =
[[-1.25  0.75  0.25]
 [ 1.75 -1.25  0.25]
 [ 1.5  -0.5  -0.5 ]]
c =
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

単位行列が得られることが確かめられた.

上記,$3\times3$行列同士を乗算するコードを以下に記す.

#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
import numpy as np
print('==========================')
print('3 * 3 array multiplication')
print('==========================')
a = np.array([[3.,1.,2.],[5.,1.,3.],[4.,2.,1.]])
b = np.array([[-1.25,0.75,0.25],[1.75,-1.25,0.25],[1.5,-0.5,-0.5]])
c = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
print('a =')
print(a)
print('b =')
print(b)

for i in range(0,3,1): 
    for j in range(0,3,1):
            for k in range(0,3,1):
                c[i,j] += a[i,k]*b[k,j]
        
print('c =')
print(c)

index.htmlに戻る