概要
有限要素法などでよく使用されるスカイライン法ですが、Pythonで実装している例があまり見つからなかったので投稿しました。また、スカイライン法に関する資料も現時点(2021/9)ではインターネット上に少ないと感じたので、情報のまとめも兼ねて投稿します。
スカイライン法とは
有限要素法において計算する行列は一般的には大規模な疎行列(ほとんどが0成分)で対称行列になります。スカイライン法はそのような大規模な疎行列に対して少ないメモリ容量で高速に解くためのアルゴリズムです。
目的は下記の式を解くことです。
Kx = f \\
K : 既知の疎行列\\
x : 未知のベクトル\\
f : 既知のベクトル\\
疎行列Kの例
K =
\begin{pmatrix}
k_{11} & k_{12} & 0 & 0 & 0 & 0\\
k_{21} & k_{22} & k_{23} & 0 & 0 & 0\\
0 & k_{32} & k_{33} & k_{34} & 0 & 0\\
0 & 0 & k_{43} & k_{44} & k_{45} & 0\\
0 & 0 & 0 & k_{54} & k_{55} & k_{56}\\
0 & 0 & 0 & 0 & k_{65} & k_{66}
\end{pmatrix}
有限要素法では、Kを剛性行列, xを変位ベクトル, fを荷重ベクトルとして解きます。
スカイライン法ではKをLDU分解し、下三角行列Lを使って前進代入、上三角行列Uを使って後退代入を行うことで未知のベクトルxを解きます。
前提知識
LDU分解(修正コレスキー分解)
行列Kが係数行列であれば、下記のように行列を分解することができます。
K = LDU \\
L : 下三角行列\\
D : 対角行列\\
U : 上三角行列\\
L=
\begin{pmatrix}
1 & 0 & 0 & 0 & \cdots\\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1 \\
\vdots & & & & \ddots \\
\end{pmatrix}
\\
D=
\begin{pmatrix}
d_{11} & 0 & 0 & 0 & \cdots\\
0 & d_{22} & 0 & 0 \\
0 & 0 & d_{33} & 0 \\
0 & 0 & 0 & d_{44} \\
\vdots & & & & \ddots
\end{pmatrix}
\\
U =
\begin{pmatrix}
1 & u_{12} & u_{13} & u_{14} & \cdots\\
0 & 1 & u_{23} & u_{24} \\
0 & 0 & 1 & u_{34} \\
0 & 0 & 0 & 1 \\
\vdots & & & & \ddots
\end{pmatrix}
ここで、Kの全ての主座小行列の行列式が0以外ならばLDU分解を行うことが可能となります。
また、Kが対称行列のときは下記の式が成り立ちます。
K = LDL^T
有限要素法における剛性行列Kは一般的に上記を満たします。
また、LDU分解を行うことで未知のベクトルxは下記のように解くことが可能となります。
LDL^Tx = f\\
y = DL^Tx \\
Ly = f
ここで、yは下三角行列Lが左にかかっているので、下記のように表せます。
Ly=
\begin{pmatrix}
1 & 0 & 0 & 0 & \cdots\\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1 \\
\vdots & & & & \ddots \\
\end{pmatrix}
\begin{pmatrix}
y_1 \\
y_2 \\
y_3 \\
y_4 \\
\vdots
\end{pmatrix}
=
\begin{pmatrix}
f_1 \\
f_2 \\
f_3 \\
f_4 \\
\vdots
\end{pmatrix}
\\
y_1 = f_1 \\
y_2 = f_2 - l_{21}y_1 \\
y_3 = f_3 - l_{31}y_1 - l_{32}y_2 \\
\vdots\\
y_n = f_n - \sum_{k=1}^{n-1}l_{nk}f_k \\
y_i = f_i - \sum_{k=1}^{i-1}l_{ik}f_k \\
上式のように$y_1$, $y_2$・・・$y_n$の順に求める方法を前進代入法とといい、容易にyを求めることが可能です。
また、xは下記のように表せるので、
y = DL^Tx \\
L^Tx = D^{-1}y
求めたyより下記の式を作ることができます。
L^Tx=
\begin{pmatrix}
1 & l_{21} & l_{31} & l_{41} & \cdots & &l_{n1}\\
0 & 1 & l_{32} & l_{42} & & & l_{n2} \\
0 & 0 & 1 & l_{43} & & & l_{n3} \\
0 & 0 & 0 & 1 & & & l_{n4}\\
\vdots & & & & \ddots & & \vdots\\
0 & \cdots & & & 0 & 1 & l_{n,n-1}\\
0 & \cdots & & & & 0 & 1
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
x_4 \\
\vdots\\
x_{n-1}\\
x_n
\end{pmatrix}
=
\begin{pmatrix}
\frac{y_1}{d_1} \\
\frac{y_2}{d_2} \\
\frac{y_3}{d_3} \\
\frac{y_4}{d_4} \\
\vdots \\
\frac{y_{n-1}}{d_{n-1}} \\
\frac{y_n}{d_n}
\end{pmatrix}
\\
x_n = \frac{y_n}{d_n}\\
x_{n-1} = \frac{y_{n-1}}{d_{n-1}} - l_{n,n-1} x_n\\
x_{n-2} = \frac{y_{n-2}}{d_{n-2}} - l_{n,n-2} x_n - l_{n-1,n-2} x_{n - 1} \\
\vdots\\
x_1 = \frac{y_1}{d_1} - \sum_{k=2}^{n} l_{k1}x_k \\
x_i = \frac{y_i}{d_i} - \sum_{k=i + 1}^{n} l_{ki}x_k
上式のように$x_n$, $x_{n-1}$・・・$x_1$の順に求める方法を後退代入法とといい、こちらも容易にxを求めることが可能です。
つまり、行列KをLDU分解することができれば、求めた行列L,Dとベクトルfによって容易に未知のベクトルxを求めることが可能となります。
また、分解後のLにはKのバンド性を保持するという特徴があります。下記に例を示します。
K =
\begin{pmatrix}
k_{11} \\
k_{21} & k_{22}\\
k_{31} & k_{32} & k_{33} & & symm\\
0 & k_{42} & k_{43} & k_{44} \\
0 & k_{52} & k_{53} & k_{54} & k_{55} \\
0 & 0 & 0 & k_{64} & k_{65} & k_{66}
\end{pmatrix}
\\
L =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0\\
l_{21} & 1 & 0 & 0 & 0 & 0\\
l_{31} & l_{32} & 1 & 0 & 0 & 0\\
0 & l_{42} & l_{43} & 1 & 0 & 0\\
0 & l_{52} & l_{53} & l_{54} & 1 & 0\\
0 & 0 & 0 & l_{64} & l_{65} & 1
\end{pmatrix}
スカイライン法のアルゴリズム
スカイライン行列の作成
スカイライン法では行列Kを1次元の配列として扱います。行列Kを1次元の配列として扱う際には行列を格納する順番が資料によって違ったりしますが、本稿では下記の順番に1次元配列に格納することにします。
\begin{pmatrix}
① & ② & ④ & ⑦ \\
② & ③ & ⑤ & ⑧ \\
④ & ⑤ & ⑥ & ⑨ \\
⑦ & ⑧ & ⑨ & ⑩ \\
& & & & \ddots
\end{pmatrix}
\\
ここで、行列Kは対称行列なので実際に1次元配列に格納する際は上三角か下三角どちらか片方だけで十分になります。本稿では下三角で格納することにします。また、格納している最中に行の初めに0が存在する場合は0の値を飛ばして格納します。下記に例を示します。
K =
\begin{pmatrix}
k_{11} & \\
k_{21} & k_{22} & & & symm \\
k_{31} & k_{32} & k_{33}\\
0 & k_{42} & k_{43} & k_{44} \\
0 & 0 & k_{53} & k_{54} & k_{55}\\
\end{pmatrix}
\\
スカイライン行列 :
\begin{bmatrix}
k_{11} & k_{21} & k_{22} & k_{31} & k_{32} & k_{33} & k_{42} & k_{43} & k_{44} &
k_{53} & k_{54} & k_{55}
\end{bmatrix}
\\
nd =
\begin{bmatrix}
0 & 1 & 3 & 6 & 9 & 12
\end{bmatrix}
ここで、ndはスカイライン行列の対角項が現れる位置を示したもので、スカイライン行列が何行何列にあるのかを把握するために必要な情報となります。ndの最初はダミーとして必ず0にしておきます。
ここで、ndの作り方について簡単に解説します。
まず初めに行列の各行において最初の0ではない要素から対角項までの要素数を並べた列を作成します。
K =
\begin{pmatrix}
k_{11} & \\
k_{21} & k_{22} & & & symm \\
k_{31} & k_{32} & k_{33}\\
0 & k_{42} & k_{43} & k_{44} \\
0 & 0 & k_{53} & k_{54} & k_{55}\\
\end{pmatrix}
\\
\begin{bmatrix}
0 & 1 & 2 & 3 & 3 & 3
\end{bmatrix}
次に左から順に総和をとった数字に置き換えることでndを作成することができます。
\begin{bmatrix}
0 & 1 & 2 & 3 & 3 & 3
\end{bmatrix}
\\
nd =
\begin{bmatrix}
0 & 1 & 3 & 6 & 9 & 12
\end{bmatrix}
Kのi行j列の要素にアクセスするためには、ndを使用して下記の番号のスカイライン行列にアクセスします。
nd_{i+1} + j -i
# スカイライン行列, スカイライン表を作成する
def MakeSkyline(self):
skylineMat = []
tmp = [0]
skylineTable = []
for i in range(self.matK.shape[0]):
skyFlg = False
cnt = 0
for j in range(i + 1):
if skyFlg == True or not self.matK[i, j] == 0:
skyFlg = True
skylineMat.append(self.matK[i, j])
cnt += 1
tmp.append(cnt)
for i in range(len(tmp)):
sum = 0
for j in range(i + 1):
sum += tmp[j]
skylineTable.append(sum)
#print(skylineTable)
#print(skylineMat)
return skylineMat, skylineTable
ndはコード中ではスカイライン表として定義しています。
LDU分解
スカイライン行列を作成したら、次にLDU分解を行います。
ここで、LDU分解を行うと、行列Kの要素は下記の数式で表すことができます。
K = LDL^T\\
k_{ij} = \sum_{k=1}^{j-1}l_{ik}d_{kk}l_{jk} + l_{ij}d_{jj} \quad (i > j)\\
k_{ii} = \sum_{k=1}^{i-1}l_{ik}^2d_{kk} + d_{ii} \quad (i = j)
上式を変形すると下記の式となります。
d_{11} = k_{11} \\
l_{ij} = \frac{1}{d_{jj}}(k_{ij} - \sum_{k=1}^{j-1}l_{ik}d_{kk}l_{jk}) \quad (i > j)\\
d_{ii} = k_{ii} - \sum_{k=1}^{i-1}l_{ik}^2d_{kk} \quad (i = j)
上式より、行列L, Dの各要素の値を求めることができます。また、上式はさらに下記の式に置き換えることが可能です。
\bar{l}_{ij} = k_{ij} - \sum_{k=1}^{j-1}l_{ik}d_{kk}l_{jk} \\
l_{ij} = \frac{\bar{l}_{ij}}{d_{jj}} \quad (i > j)\\
l_{ik} = \frac{\bar{l}_{ik}}{d_{kk}}\\
d_{ii} = k_{ii} - \sum_{k=1}^{i-1}l_{ik}^2d_{kk} = k_{ii} - \sum_{k=1}^{i-1}l_{ik}\bar{l}_{ik} \quad (i = j)
まとめると
\bar{l}_{11} = k_{11}\\
\bar{l}_{ij} = k_{ij} - \sum_{k=1}^{j-1}l_{ik}d_{kk}l_{jk} \\
l_{ij} = \frac{\bar{l}_{ij}}{d_{jj}} \quad (i > j)\\
d_{ii} = k_{ii} - \sum_{k=1}^{i-1}l_{ik}\bar{l}_{ik} \quad (i = j)\\
となり、$\bar{l}_{ij}$を求めることで、三重積から二重積に減らすことができます。
# LDU分解を行う
def LDUdecomp(self, skylineMat, skylineTabel):
# 初期のL, D行列を生成する(行数 : len(skylineTabel) - 1)
matL = np.matrix(np.eye(len(skylineTabel) - 1))
matD = np.matrix(np.zeros((len(skylineTabel) - 1, len(skylineTabel) - 1)))
# LDU分解を行う
matD[0, 0] = skylineMat[0]
matL_bar = np.matrix(np.zeros((len(skylineTabel) - 1, len(skylineTabel) - 1)))
matL_bar[0, 0] = skylineMat[0]
for i in range(1, len(skylineTabel) - 1):
for j in range(i + 1):
skyNo = skylineTabel[i + 1] + j - i # スカイライン行列の番号を計算する
# matK_ijが0のとき
if skyNo <= skylineTabel[i]:
continue
# matL_barを計算する
llbarSum = 0.0
if not j==0:
vecL_bar = matL_bar[i : i + 1, 0 : j]
vecL = matL[j : j + 1, 0 : j]
llbarSum = np.dot(vecL_bar, vecL.T)
matK_ij = skylineMat[skyNo - 1] # skyNoはスカイライン行列の始めを1として計算しているので-1する
matL_bar[i, j] = matK_ij - llbarSum
# matLを計算する
if not matD[j, j] == 0.0 and not i == j:
matL[i, j] = matL_bar[i, j] / matD[j, j]
# matDを計算する
matD[i ,i] = matL_bar[i, i]
#print("matL\n", matL)
#print("matL_bar\n", matL_bar)
#print("matD\n", matD)
#print("matL * matD * matL.T\n", matL * matD * matL.T)
return matL, matD
LDU分解を行う際はバンド性が保持されることを利用して、$k_{ij} = L_{ij} = 0$を利用しています。
前進代入、後退代入
最後に前進代入、後退代入を行いアルゴリズムを終了します。計算方法は上述を参照ください。
# 前進代入を行う
vecy = np.zeros(len(skylineTabel) - 1)
for i in range(len(vecy)):
lySum = 0
for k in range(i):
lySum += matL[i, k] * vecy[k]
vecy[i] = self.vecf[i] - lySum
#print(vecy)
# 後退代入を行う
vecx = np.zeros(len(skylineTabel) - 1)
for i in reversed(range(len(vecx))):
lxSum = 0
for k in range(i + 1, len(vecx)):
lxSum += matL[k, i] * vecx[k]
vecx[i] = vecy[i] / matD[i, i] - lxSum
Pythonコード全体
import numpy as np
import numpy.linalg as LA
import time
class Skyline:
def __init__(self, matK, vecf):
# 対称行列かチェックする
if not np.array_equal(matK, matK.T):
print("対称行列を入力してください")
exit()
# 主座小行列の行列式が0でないことをチェックする
# for i in range(matK.shape[0]):
# matA = matK[0:i + 1, 0:i + 1]
# if LA.det(matA) == 0:
# print("この行列はLDU分解することができません")
# exit()
# ベクトルの次元が行列の次元と一致しているかチェックする
if not matK.shape[0] == vecf.shape[0]:
print("行列とベクトルの次元が一致していません")
exit()
self.matK = matK
self.vecf = vecf
# スカイライン法で解く
def Solve(self):
# スカイライン行列を作成する
skylineMat, skylineTabel = self.MakeSkyline()
# LDU分解を行う
matL, matD = self.LDUdecomp(skylineMat, skylineTabel)
# 前進代入を行う
vecy = np.zeros(len(skylineTabel) - 1)
for i in range(len(vecy)):
lySum = 0
for k in range(i):
lySum += matL[i, k] * vecy[k]
vecy[i] = self.vecf[i] - lySum
#print(vecy)
# 後退代入を行う
vecx = np.zeros(len(skylineTabel) - 1)
for i in reversed(range(len(vecx))):
lxSum = 0
for k in range(i + 1, len(vecx)):
lxSum += matL[k, i] * vecx[k]
vecx[i] = vecy[i] / matD[i, i] - lxSum
return vecx
# スカイライン行列, スカイライン表を作成する
def MakeSkyline(self):
skylineMat = []
tmp = [0]
skylineTable = []
for i in range(self.matK.shape[0]):
skyFlg = False
cnt = 0
for j in range(i + 1):
if skyFlg == True or not self.matK[i, j] == 0:
skyFlg = True
skylineMat.append(self.matK[i, j])
cnt += 1
tmp.append(cnt)
for i in range(len(tmp)):
sum = 0
for j in range(i + 1):
sum += tmp[j]
skylineTable.append(sum)
#print(skylineTable)
#print(skylineMat)
return skylineMat, skylineTable
# LDU分解を行う
def LDUdecomp(self, skylineMat, skylineTabel):
# 初期のL, D行列を生成する(行数 : len(skylineTabel) - 1)
matL = np.matrix(np.eye(len(skylineTabel) - 1))
matD = np.matrix(np.zeros((len(skylineTabel) - 1, len(skylineTabel) - 1)))
# LDU分解を行う
matD[0, 0] = skylineMat[0]
matL_bar = np.matrix(np.zeros((len(skylineTabel) - 1, len(skylineTabel) - 1)))
matL_bar[0, 0] = skylineMat[0]
for i in range(1, len(skylineTabel) - 1):
for j in range(i + 1):
skyNo = skylineTabel[i + 1] + j - i # スカイライン行列の番号を計算する
# matK_ijが0のとき
if skyNo <= skylineTabel[i]:
continue
# matL_barを計算する
llbarSum = 0.0
if not j==0:
vecL_bar = matL_bar[i : i + 1, 0 : j]
vecL = matL[j : j + 1, 0 : j]
llbarSum = np.dot(vecL_bar, vecL.T)
matK_ij = skylineMat[skyNo - 1] # skyNoはスカイライン行列の始めを1として計算しているので-1する
matL_bar[i, j] = matK_ij - llbarSum
# matLを計算する
if not matD[j, j] == 0.0 and not i == j:
matL[i, j] = matL_bar[i, j] / matD[j, j]
# matDを計算する
matD[i ,i] = matL_bar[i, i]
#print("matL\n", matL)
#print("matL_bar\n", matL_bar)
#print("matD\n", matD)
#print("matL * matD * matL.T\n", matL * matD * matL.T)
return matL, matD
def test1():
matK1 = np.matrix([[1, 2, 3],
[2, 6, 0],
[3, 0, 9]])
vecf1 = np.array([1, 2, 3])
print("matK1\n", matK1)
print("vecf1\n", vecf1)
print("vecx1\n", LA.inv(matK1) @ vecf1)
skyline = Skyline(matK1, vecf1)
print(skyline.Solve())
print("\n")
matK2 = np.matrix([[4, -3, 1, 0],
[-3, 5, -3, 1],
[1, -3, 5, -3],
[0, 1, -3, 4]])
vecf2 = np.array([1, 2, 3, 4])
print("matK2\n", matK2)
print("vecf2\n", vecf2)
print("vecx2\n", LA.inv(matK2) @ vecf2)
skyline = Skyline(matK2, vecf2)
print(skyline.Solve())
print("\n")
matK3 = np.matrix([[1, 2, 0, 0],
[2, 3, 4, 6],
[0, 4, 5, 7],
[0, 6, 7, 8]])
vecf3 = np.array([1, 2, 3, 4])
print("matK3\n", matK3)
print("vecf3\n", vecf3)
print("vecx3\n", LA.inv(matK3) @ vecf3)
skyline = Skyline(matK3, vecf3)
print(skyline.Solve())
print("\n")
matK4 = np.matrix([[1, 2, 0, 0],
[2, 0, 4, 6],
[0, 4, 5, 7],
[0, 6, 7, 0]])
vecf4 = np.array([1, 2, 3, 4])
print("matK4\n", matK4)
print("vecf4\n", vecf4)
print("vecx4\n", LA.inv(matK4) @ vecf4)
skyline = Skyline(matK4, vecf4)
print(skyline.Solve())
print("\n")
def test2(length):
# ランダムな対称行列を作成する
flg = True
while flg:
matR = np.random.random_integers(0,2,size=(length, length))
matK = np.matrix(np.zeros((length, length)))
for i in range(length):
for j in range(i + 1):
matK[i, j] = matR[i, j]
matK[j, i] = matK[i, j]
matK[0, 0] = 1 # この値が0になると失敗するため、必ず1にする
flg = False
for i in range(matK.shape[0]):
matA = matK[0:i + 1, 0:i + 1]
if LA.det(matA) == 0:
flg = True
# ランダムなベクトルを作成する
vecf = np.random.random_integers(0,100,size=(length))
print(matK)
print(vecf)
# Numpyで計算した時の処理時間を計測する
st = time.time()
vecx = LA.inv(matK) @ vecf
print(vecx)
end = time.time()
print("Numpyの処理時間 : %f[s]" % (end - st))
# スカイライン法の処理時間を計測する
st = time.time()
skyline = Skyline(matK, vecf)
skyline.Solve()
end = time.time()
print("スカイライン法の処理時間 : %f[s]" % (end - st))
if __name__ == '__main__':
# 値が正しいかチェックする
test1()
# 処理時間を比較する
test2(100)
テスト関数として入力チェックのためのtest1関数、numpyとの処理速度を比較するためのtest2関数を用意しています。
実行してみると分かりますが、残念ながらnumpyには処理速度で負けています。numpyはネイティブコード上で計算するの対してこちらはインタプリタ上なので、速度で上回ることはやはり不可能かなと思ってます。なので、スカイライン法を実装する際はC言語などのネイティブコード上で計算可能な言語を使用した方が良いかもしれません。
#参考文献
コレスキー分解 - PukiWiki for PBCG Lab
Microsoft Word - マトリックス法第5章.doc
大規模疎行列に対する連立一次方程式のソルバー：Skyline法 | 人力飛行機設計日記＠mtk_birdman