-Python- Numpy 逆行列,行列式,階数の計算

Numpyには,行列計算の機能(関数)が備わっています.
以下に,逆行列計算,行列式計算,階数計算のプログラム例を示します.

逆行列計算
逆行列を計算するプログラムです.元の行列(list_a)に対して np.linalg.inv 関数を使って逆行列list_b)を計算しています.mat_c = mat_a * mat_b によって,元の行列と計算された逆行列を掛け合わせると単位行列が得られています.

import numpy as np
list_a = (0, 4)
list_b = (1, 5)
mat_a = np.matrix([list_a, list_b])
mat_b = np.linalg.inv(mat_a)
 
print(mat_a)
print(mat_b)
 
mat_c = mat_a * mat_b

 

print(mat_c)
 
[[0 4]
 [1 5]]
[[-1.25  1.  ]
 [ 0.25  0.  ]]
[[1. 0.]
 [0. 1.]]
 

ただし,逆行列の計算には注意が必要です.
下記のプログラム例のような場合は,元の行列と計算した逆行列との積が単位行列になりません.この原因は,情報落ち,桁落ちが原因だと思われます(正確には,何が原因なのかは良くわかりません.

import numpy as np
 
list_a = (1., 2.)
list_b = (3., 4.)
 
mat_a = np.matrix([list_a,list_b])
mat_b = np.linalg.inv(mat_a)
 
print(mat_a)
print(mat_b)
 
mat_c = mat_a * mat_b
 
print(mat_c)
 
実行結果は以下のようになります.
[[1. 2.]
 [3. 4.]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[[1.00000000e+00 1.11022302e-16]
 [0.00000000e+00 1.00000000e+00]]
mat_b = np.linalg.inv(mat_a)
の部分を,以下のようにsolve() を使うように変更してみても結果は変わりません.
mat_b = np.linalg.solve(mat_a, np.eye(2))
 
実行結果
[[1. 2.]
 [3. 4.]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[[1.00000000e+00 1.11022302e-16]
 [0.00000000e+00 1.00000000e+00]]

逆行列の計算には,用心しなければならないという事が明らかになりました.

現段階では,この解決法はわかりませんが,計算誤差と割り切って,以下のように丸めてしまえば,一応正しい結果が得られます.
 
import numpy as np
 
list_a = (1., 2.)
list_b = (3., 4.)
 
mat_a = np.matrix([list_a,list_b])
mat_b = np.linalg.inv(mat_a)
 
print(mat_a)
print(mat_b)
 
mat_c = mat_a * mat_b
mat_c2 = np.round(mat_c, decimals=5)

print(mat_c2)

実行結果
[[1. 2.]
 [3. 4.]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[[1. 0.]
 [0. 1.]]
 
行列式の計算
行列式の計算を行うプログラム例を以下に示します.np.linalg.det は,行列式を求める関数です.
 
import numpy as np
list_a = (-1, 4, 1)
list_b = (1, 1, 2)
list_c = (0, 9, 1)
 
mat_a = np.matrix([list_a, list_b, list_c])
det_a = np.linalg.det(mat_a)
 
print(mat_a)
print(det_a)
 
実行結果
[[-1  4  1]
 [ 1  1  2]
 [ 0  9  1]]
22.000000000000004
 
正しい答えは22.なのですが,これも末尾に'4'がついています.
先の逆行列と同様に,とても小さな誤差なので,良いとみなすこともできなくはないですが,これまた注意が必要ということはわかりました....

これも,先と同様に det_a の計算後に以下の行を追加して...
 
det_a2 = np.round(det_a, decimals=5)
...
print(det_a2)

として det_a2 を表示させると以下のように正しい結果が得られます(これが良いかどうかは別として...)

実行結果
22.0

・階数(rank)の計算
 
import numpy as np
list_a = (1, 1, -3)
list_b = (-1, 0, 5)
list_c = (0, 3, 6)
 
mat_a = np.matrix([list_a, list_b, list_c])
rank_a = np.linalg.matrix_rank(mat_a)
 
print(mat_a)
print(rank_a)
 
実行結果は以下のようになります.
[[ 1  1 -3]
 [-1  0  5]
 [ 0  3  6]]
2