前回
前回では行列をLU分解した。
今回はLU分解した行列の逆行列を求める。
前回作成したlu_decomposition関数は引き続き使用する。
下三角行列Lの逆行列
ソースコード
def inverse_l(l):
length = len(l)
il = np.eye(length)
for i in range(1,length):
for j in range(i):
il[i,j] = -l[i,j]
for k in range(j+1,i):
il[i,j] -= l[i,k] * il[k,j]
return il
np.set_printoptions(linewidth=150)
a = np.random.rand(5,5)
l, u = lu_decomposition(a)
il = inverse_l(l)
print('L\n', l)
print('L^{-1}\n', il)
print('L * L^{-1}\n', l.dot(il))
print('L^{-1} * L\n', il.dot(l))
出力結果
L
[[ 1. 0. 0. 0. 0. ]
[ 1.14616816 1. 0. 0. 0. ]
[ 1.34499533 2.23780925 1. 0. 0. ]
[ 1.10881862 0.92510343 -0.19976126 1. 0. ]
[ 0.59746136 -0.5427907 0.222358 1.1189475 1. ]]
L^{-1}
[[ 1. 0. 0. 0. 0. ]
[-1.14616816 1. 0. 0. 0. ]
[ 1.21991038 -2.23780925 1. 0. 0. ]
[ 0.19519631 -1.37213103 0.19976126 1. 0. ]
[-1.70926204 2.57572808 -0.44588037 -1.1189475 1. ]]
L * L^{-1}
[[1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00]
[5.55111512e-17 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
L^{-1} * L
[[ 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -1.11022302e-16 0.00000000e+00 1.00000000e+00 0.00000000e+00]
[ 2.22044605e-16 2.22044605e-16 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
上三角行列Uの逆行列
ソースコード
def inverse_u(u):
length = len(u)
iu = np.zeros((length, length))
for i in range(length):
for j in range(i, length):
tmp = 0 if i != j else 1
for k in range(i,j):
tmp -= iu[i,k] * u[k,j]
iu[i,j] = tmp / u[j,j]
return iu
np.set_printoptions(linewidth=150)
a = np.random.rand(5,5)
l, u = lu_decomposition(a)
iu = inverse_u(u)
print('U\n', u)
print('U^{-1}\n', iu)
print('U * U^{-1}\n', u.dot(iu))
print('U^{-1} * U\n', iu.dot(u))
出力結果
U
[[ 0.31491714 0.52824061 0.35026388 0.5302676 0.15019698]
[ 0. -1.06948835 -0.4657366 -0.3557189 0.56946014]
[ 0. 0. -0.4410975 -0.38212418 -0.01388272]
[ 0. 0. 0. -1.04239703 -0.48815709]
[ 0. 0. 0. 0. 0.40882289]]
U^{-1}
[[ 3.17543845 1.56840936 0.86551327 0.76284342 -2.41103133]
[ 0. -0.93502655 0.98725584 -0.04283178 1.28480456]
[ 0. 0. -2.26707247 0.83106838 0.91535678]
[ 0. 0. 0. -0.95932737 -1.14548983]
[ 0. 0. 0. 0. 2.44604702]]
U * U^{-1}
[[ 1.00000000e+00 -3.48270747e-17 1.14348072e-17 -9.71222523e-17 -6.57160923e-18]
[ 0.00000000e+00 1.00000000e+00 2.84641731e-17 7.23787024e-17 -1.85570429e-16]
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00 7.35095055e-17 4.55219411e-17]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 -7.33637720e-17]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
U^{-1} * U
[[ 1.00000000e+00 1.80153811e-17 -7.71326778e-17 -2.31860133e-16 -7.43377334e-17]
[ 0.00000000e+00 1.00000000e+00 -2.31203199e-18 -6.05873099e-18 -1.64000591e-17]
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00 5.35325757e-17 -3.82652839e-18]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 6.88020856e-18]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
元の行列の逆行列
ソースコード
a = np.random.rand(5,5)
l, u = lu_decomposition(a)
il = inverse_l(l)
iu = inverse_u(u)
ia = iu.dot(il)
print('a\n', a)
print('a^{-1}\n', ia)
print('a * a^{-1}\n', a.dot(ia))
print('a^{-1} * a\n', ia.dot(a))
出力結果
a
[[0.77920088 0.1471498 0.12865165 0.85940836 0.08181253]
[0.16269642 0.9075763 0.99422554 0.54888386 0.23620424]
[0.56614793 0.69306587 0.77175149 0.91403137 0.21029721]
[0.66274032 0.77909366 0.37639368 0.57952536 0.88788831]
[0.07098181 0.87230157 0.22465406 0.82722447 0.38804975]]
a^{-1}
[[ 9.57427605 10.08736519 -14.45904346 0.21468055 -0.8140454 ]
[ 13.03504096 15.31052423 -21.92218184 -0.44649226 0.83434764]
[ -6.65020152 -6.45261992 10.92073055 0.18126003 -1.00330376]
[ -7.73321467 -9.66683796 13.65916377 -0.2849802 0.76425409]
[-10.71766426 -11.91899126 16.48371559 1.46697502 -0.19799977]]
a * a^{-1}
[[ 1.00000000e+00 -3.76312255e-16 3.30140460e-15 2.00929848e-16 7.86733982e-18]
[ 2.97110501e-15 1.00000000e+00 4.41381135e-15 -1.66752070e-16 2.64946104e-16]
[ 2.85341323e-15 -1.99495978e-15 1.00000000e+00 5.38596575e-18 1.03655022e-16]
[-8.49519864e-16 -7.10786885e-15 9.97693462e-15 1.00000000e+00 9.09196633e-17]
[ 1.00869260e-14 9.77628210e-15 -7.31079258e-15 -7.57730822e-16 1.00000000e+00]]
a^{-1} * a
[[ 1.00000000e+00 -2.20867121e-15 -3.25909098e-15 -1.77358142e-15 1.70542576e-16]
[ 2.43341874e-15 1.00000000e+00 -1.31749913e-15 1.89853211e-15 -5.99409485e-18]
[-9.42007514e-16 -3.74486862e-15 1.00000000e+00 -2.35845894e-15 -6.10110897e-16]
[ 4.81060848e-15 7.47562435e-15 7.94852745e-15 1.00000000e+00 1.29595009e-15]
[ 5.57567989e-16 1.46424490e-15 2.93137817e-15 5.09057328e-16 1.00000000e+00]]
まとめ
前回の記事と合わせて、任意の正方行列の逆行列を求めるソースコードを作成した。
逆行列が存在しない行列に対してはエラーを起こしてしまうが、私が使用する場面では逆行列が存在しないことはほとんど起こらないので、エラー処理だけしておいて、P行列は作成しない。