ライブラリが豊富でない言語でLU分解したかったため、コードを作成することになった。
個人的に見やすく、書きやすいPythonでコードを作成してから、同じコードを別の言語で作成することにした。
以下がその作成したPythonコード。

ソースコード

import numpy as np

def lu_decomposition(a):
  length = len(a)
  l = np.eye(length)
  u = np.zeros((length,length))

  for i in range(length):
    for j in range(i):
      l[i,j] = a[i,j]
      for k in range(j):
        l[i,j] = l[i,j] - l[i,k] * u[k,j]
      l[i,j] = l[i,j] / u[j,j]

    for j in range(i, length):
      u[i,j] = a[i,j]
      for k in range(i):
        u[i,j] = u[i,j] - l[i,k] * u[k,j]
  return l,u

a = np.random.rand(5,5)
l, u = lu_decomposition(a)
print(a)
print(l)
print(u)
print(l.dot(u))
print(l.dot(u) == a)

出力結果

[[0.81391896 0.1476644  0.14297368 0.29154456 0.20082125]
 [0.58814889 0.47336317 0.85443161 0.50634996 0.73427392]
 [0.3060707  0.74831109 0.80533407 0.42147219 0.98555506]
 [0.78881249 0.7541937  0.67006895 0.85316669 0.52989888]
 [0.05798169 0.22727099 0.64835561 0.50293417 0.91297261]]
[[ 1.          0.          0.          0.          0.        ]
 [ 0.72261358  1.          0.          0.          0.        ]
 [ 0.37604566  1.88944704  1.          0.          0.        ]
 [ 0.9691536   1.66662873  1.07893841  1.          0.        ]
 [ 0.07123767  0.59115362 -0.29079903  0.68459484  1.        ]]
[[ 0.81391896  0.1476644   0.14297368  0.29154456  0.20082125]
 [ 0.          0.36665888  0.75111688  0.2956759   0.58915775]
 [ 0.          0.         -0.66762614 -0.24682583 -0.20314528]
 [ 0.          0.          0.          0.34414316 -0.42745376]
 [ 0.          0.          0.          0.          0.78394202]]
[[0.81391896 0.1476644  0.14297368 0.29154456 0.20082125]
 [0.58814889 0.47336317 0.85443161 0.50634996 0.73427392]
 [0.3060707  0.74831109 0.80533407 0.42147219 0.98555506]
 [0.78881249 0.7541937  0.67006895 0.85316669 0.52989888]
 [0.05798169 0.22727099 0.64835561 0.50293417 0.91297261]]
[[ True  True  True  True  True]
 [ True  True  True  True  True]
 [ True False  True False False]
 [ True  True  True  True False]
 [ True False  True  True  True]]

誤差が多少出てしまうのはご愛嬌。

LU分解では使用するメモリを省略するために、L行列とU行列を一つの行列にまとめることがあるが、今回は動作確認がしたかったため、別の行列にしている。

参考

次回