はじめに
本記事はコンピュータで数値計算を行うときに、計算誤差により意図しない結果となることがあるよというお話です。
ある文系大学生（青色coder）に、PythonのNumPyが行列計算に便利だと力説しているときに、逆行列の話になりました。そして、試しに逆行列を求めてるとなぜか値がおかしい。原因がよくわかわない。あせる私。このままでは、理系大学院生としてどうなのかと思うので考察します。
問題点
適当に以下のコードで行列を作成して、逆行列を計算すると
import numpy as np
A = np.arange(9).reshape(3, 3)
print(A@np.linalg.inv(A))
#---> LinAlgError: Singular matrix
print(np.linalg.det(A))
#---> 0.0
となり、正則行列（逆行列が定義できる行列のこと）ではないので逆行列を計算できない。そこで、以下のコードを実行すると
A = np.arange(1,10).reshape(3, 3)
print(A@np.linalg.inv(A))
#--->
# [[ 0. 1. 0.]
# [ 0. 2. 0.]
# [-4. 3. 2.]]
なぜか単位行列とかけ離れた行列になります。
考察
行列式を計算すると
print(np.linalg.det(A))
#---> -9.51619735392994e-16
となり、正則行列ではありますが、とても小さな値となっています。
条件数の逆数rcondを計算すると、
print(1/np.linalg.cond(A))
#---> 1.9793046108741972e-17
となり、とても小さく特異性が高い行列だとわかります1。（特異性が高いというのは逆行列が存在しない特異行列に近いとということ）
特異性が高い行列だと、逆行列が正しく求められないみたいです。
対角成分に微小値を足すと近似的に求められます。
A = np.arange(1,10).reshape(3, 3)
A = A + 1e-8*np.eye(3)
print(A@np.linalg.inv(A))
#--->
#[[ 1.00000000e+00 7.10542736e-15 3.55271368e-15]
# [-7.10542736e-15 1.00000000e+00 7.10542736e-15]
# [-4.83857773e-15 -2.79601156e-14 1.00000000e+00]]
print(np.round(A@np.linalg.inv(A)))
#--->
#[[1.00000000e+00 0.00000000e+00 0.00000000e+00]
# [0.00000000e+00 1.00000000e+00 0.00000000e+00]
# [4.13602329e-09 1.88513379e-08 1.00000001e+00]]
となり、単位行列となります。
特異性が高い行列だと逆行列の精度は低くなるみたいです。気づかないとこわいですね。
ちなみに、MATLAB(数値解析ソフトウェアで、配列のインデックスがなぜか1から始まる謎言語1)だと、
「警告: 行列は、特異行列に近いか、正しくスケーリングされていません。結果は不正確になる可
能性があります。RCOND = 9.251859e-18。 」
と警告が出ます。助かります。
まとめ
逆行列を計算するときは、特異性に気をつけないといけません。
コンピュータでも計算誤差は生じます。
おまけ
予想外の結果が出るものを書いておきます。計算誤差ではなく、仕様のものもあります。
1.
True or False
print((0.1+0.2)==0.3)
2.
print(round(0.5))
3.
print(round(2.675, 2))
print(round(2.665, 2))
4.
import numpy as np
A = np.ones(1)
# --> [1]
B = A
A[0] = 2
print(A)
print(B)
結果と考察は時間の都合上後で追記します。
注釈
行列の特異性を見るにはdetではなくrcondを見る必要があります。
以下反例