前回

前回では行列を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行列は作成しない。