LU分解(LU decomposition)の原理の理解とPythonを用いた実装

当記事ではLU分解(LU decomposition)の原理の理解とPythonを用いた実装に関して取り扱います。LU分解は方程式を順に解く手法であり、手計算を行うと大変なのでPythonなどで計算を行うとスムーズですが、$L$や$U$に関して解くところが複雑であるのでなるべくわかりやすく取りまとめました。

・参考
Wikipedia LU分解

LU分解の原理

行列$L$と$U$の定義

以下、$n$次正方行列の$A$を下記のような下三角行列$L$と上三角行列$U$を用いて$A=LU$のように分解を行うことを考える。$L$はLower、$U$はUpperにそれぞれ対応すると考えれば良い。
$$
\large
\begin{align}
L &= \left(\begin{array}{cccc} 1 & 0 & 0 & \cdots \\ l_{21} & 1 & 0 & \cdots \\ l_{31} & l_{32} & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) \\
U &= \left(\begin{array}{cccc} u_{11} & u_{12} & u_{13} & \cdots \\ 0 & u_{22} & u_{23} & \cdots \\ 0 & 0 & u_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right)
\end{align}
$$

ここで上記の$L$の対角成分を$1$、$U$の対角成分を$u_{ii}$とおいたことに注意が必要である。また、行列$A$は下記のように表す。
$$
\large
\begin{align}
A &= \left(\begin{array}{cccc} a_{11} & a_{12} & a_{13} & \cdots \\ a_{21} & a_{22} & a_{23} & \cdots \\ a_{31} & a_{32} & a_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right)
\end{align}
$$

下三角行列の$L$と上三角行列の要素の導出

前項で行なった$A, L, U$の定義を元に$A=LU$を要素表記で計算を行う。
$$
\large
\begin{align}
& A = LU \\
& \left(\begin{array}{cccc} a_{11} & a_{12} & a_{13} & \cdots \\ a_{21} & a_{22} & a_{23} & \cdots \\ a_{31} & a_{32} & a_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) = \left(\begin{array}{cccc} 1 & 0 & 0 & \cdots \\ l_{21} & 1 & 0 & \cdots \\ l_{31} & l_{32} & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) \left(\begin{array}{cccc} u_{11} & u_{12} & u_{13} & \cdots \\ 0 & u_{22} & u_{23} & \cdots \\ 0 & 0 & u_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) \\
&= \left(\begin{array}{cccc} u_{11} & u_{12} & u_{13} & \cdots \\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} & \cdots \\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right)
\end{align}
$$

上記の結果を元に、$a_{ij}=(LU)_{ij}$や$l_{ij}, u_{ij}$は下記のように$i$と$j$に関して場合分けして考えることができる。
・$i=1, j=1$
$$
\large
\begin{align}
a_{11} &= u_{11} \\
l_{11} &= 1 \\
u_{11} &= a_{11}
\end{align}
$$

・$i=1, j>1$
$$
\large
\begin{align}
a_{1j} &= u_{1j} \\
l_{1j} &= 0 \\
u_{1j} &= a_{1j}
\end{align}
$$

・$i>1, j=1$
$$
\large
\begin{align}
a_{i1} &= l_{i1} u_{11} \\
l_{i1} &= \frac{a_{i1}}{u_{11}} \\
u_{i1} &= 0
\end{align}
$$

・$i=j>1$
$$
\large
\begin{align}
a_{ii} &= \sum_{k=1}^{i-1} (l_{ik}u_{ki}) + u_{ii} \\
l_{ii} &= 1 \\
u_{ii} &= a_{ii} – \sum_{k=1}^{i-1} (l_{ik}u_{ki})
\end{align}
$$

・$i>j>1$
$$
\large
\begin{align}
a_{ij} &= \sum_{k=1}^{j-1} (l_{ik}u_{kj}) + l_{ij}u_{jj} \\
l_{ij} &= \frac{\displaystyle a_{ij} – \sum_{k=1}^{j-1} l_{ik}u_{kj}}{u_{jj}} \\
u_{ij} &= 0
\end{align}
$$

・$j>i>1$
$$
\large
\begin{align}
a_{ij} &= \sum_{k=1}^{i-1} (l_{ik}u_{kj}) + u_{ij} \\
l_{11} &= 0 \\
u_{11} &= a_{ij} – \sum_{k=1}^{i-1} l_{ik}u_{kj}
\end{align}
$$

上記の$l_{ij}, u_{ij}$の計算にあたって、$L, U$の要素が用いられるが、どれも「$i$行より上の行」か「$i$行かつ$j$列より左の列」の要素が対応するので、$i$に関する繰り返し処理の中に$j$に関する繰り返し処理を行うことで$L, U$の要素を用いることができる。

ここで$L, U$の要素に関して式から消去を行うことも本来的には可能だが、式変形が複雑であるので、次節ではここで取り扱ったように$L, U$の要素を徐々に確定する計算処理を用いる。

LU分解の実装

LU分解の実装は下記のように行うことができる。

import numpy as np

def LU_decomposition(X):
    n = X.shape[0]
    L,U = np.zeros([n,n]), np.zeros([n,n])
    for i in range(n):
        for j in range(n):
            if i==0:
                if j==0:
                    L[i,j] = 1
                U[i,j] = X[i,j]
            elif j==0:
                L[i,j] = X[i,j]/U[0,0]
            elif i==j:
                L[i,j] = 1
                U[i,j] = X[i,j] - np.dot(L[i,:j],U[:i,j])
            elif i>j:
                L[i,j] = (X[i,j] - np.dot(L[i,:j],U[:j,j]))/U[j,j]
            elif j>i:
                U[i,j] = X[i,j] - np.dot(L[i,:i],U[:i,j])
    return L, U

A1 = np.array([[2., 1.], [1., 2.]])
L1, U1 = LU_decomposition(A1)

print(L1)
print(U1)
print(np.dot(L1,U1))

・実行結果

> print(L1)
[[ 1.   0. ]
 [ 0.5  1. ]]
> print(U1)
[[ 2.   1. ]
 [ 0.   1.5]]
> print(np.dot(L1,U1))
[[ 2.  1.]
 [ 1.  2.]]

上記では$2$次正方行列に関して取り扱ったが、より複雑な例に関しても確認するにあたって、下記では$5$次正方行列に関してLU分解を行う。

A2 = np.array([[2., 1., 2., 4., 6.], [1., 2., 3., 2., 1.], [5., 1., 2, 3., 4.], [6., 1., 2., 3, 5.], [3., 1., 2., 6., 3.]])
L2, U2 = LU_decomposition(A2)

print(L2)
print(U2)
print(np.dot(L2,U2)-A2)

・実行結果

> print(L2)
[[ 1.          0.          0.          0.          0.        ]
 [ 0.5         1.          0.          0.          0.        ]
 [ 2.5        -1.          1.          0.          0.        ]
 [ 3.         -1.33333333  1.33333333  1.          0.        ]
 [ 1.5        -0.33333333  0.33333333  7.          1.        ]]
> print(U2)
[[  2.           1.           2.           4.           6.        ]
 [  0.           1.5          2.           0.          -2.        ]
 [  0.           0.          -1.          -7.         -13.        ]
 [  0.           0.           0.           0.33333333   1.66666667]
 [  0.           0.           0.           0.         -14.        ]]
> print(np.dot(L2,U2)-A2)
[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]]

上記のLは下三角行列、Uは上三角行列であることがそれぞれ確認できる。また、print(np.dot(L2,U2)-A2)では$LU$によって$A$が復元できたかについて確認を行なった。全ての要素が0であることが確認できることから、適切に行列の分解が行えたことが確認できる。