LU分解を用いた連立方程式の解法と逆行列の導出の仕組みとPythonでのプログラミング

当記事ではLU分解を用いた「$1$次連立方程式の解法」と「逆行列」の導出の仕組みに関して詳しく確認し、確認した式を元にPythonでプログラムを作成します。特に逆行列は行列を扱う上で様々なところで用いられるので、仕組みの概要とプログラムの対応を理解しておくと良いと思います。

・参考
Wikipedia LU分解

仕組みの確認

LU分解の仕組みとPythonを用いたプログラミング

上記で詳しく取り扱った。下記で定義するLU_decompositionは当記事の実装で用いるので、抜粋を行なった。

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)-A1)

・実行結果

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

LU分解を用いた連立方程式の解法

$n$元$1$次連立方程式は$n$次正方行列の$A$と$n$次元ベクトル$x$と$b$を用いて$Ax=b$のように行列表記することができる。$Ax=b$は下記のように行列表記できる。
$$
\large
\begin{align}
Ax &= b \\
\left(\begin{array}{ccc} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{nn} \end{array} \right) \left(\begin{array}{c} x_{1} \\ \vdots \\ x_{n} \end{array} \right) &= \left(\begin{array}{c} b_{1} \\ \vdots \\ b_{n} \end{array} \right)
\end{align}
$$

ここで上記の式では$a_{11}, x_{1}, b_{1}$のように$3$種類の文字を用いたが、「連立方程式を解く」という視点で考えるときは$a_{ij}$と$b_{k}$には何らかの数字が入り、これらの値に対応する$x_{1},…,x_{n}$を導出するというのが一般的な連立方程式を解く流れである。

さて、連立方程式を解くにあたっては、$Ax=b$より$x=A^{-1}b$を用いるというのが一般的な流れであるが、逆行列の公式を用いることのできる$2 \times 2$行列以外は逆行列$A^{-1}$の計算が必要である。ここでは$A^{-1}$が得られていない場合のLU分解を用いた解法に関して確認を行う。

$A=LU$のようにLU分解を行うことができるとき、$Ax=b$は$LUx=b$と表すことができる。ここで$Ux=y$とおくと、$Ly=b$と$Ux=y$のように二段階で連立方程式を解き、$x$を求めることもできる。この計算は$L, U$が三角行列であることによりシンプルに計算を行うことができる。

以下、$Ly=b$と$Ux=y$に関してそれぞれ確認を行う。

$Ly=b$を$y$に関して解く

$$
\large
\begin{align}
b &= Ly \\
\left(\begin{array}{c} b_{1} \\ b_{2} \\ b_{3} \\ \vdots \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}{c} y_{1} \\ y_{2} \\ y_{3} \\ \vdots \end{array} \right) \\
&= \left(\begin{array}{c} y_{1} \\ l_{21}y_{1} + y_{2} \\ l_{31}y_{1} + l_{32}y_{2} + y_{3} \\ \vdots \end{array} \right)
\end{align}
$$

上記を元に$i=1$と$i>1$に場合分けし、$y_{i}$に関して方程式を解くと下記のように表せる。
・$i=1$
$$
\large
\begin{align}
b_{1} &= y_{1} \\
y_{1} &= b_{1}
\end{align}
$$

・$i>1$
$$
\large
\begin{align}
b_{i} &= y_{i} + \sum_{k=1}^{i-1} l_{ik}y_{k} \\
y_{i} &= b_{i} – \sum_{k=1}^{i-1} l_{ik}y_{k}
\end{align}
$$

上記を$i=1$から$n$まで昇順に計算することでベクトル$y$の各要素の値が得られる。

$Ux=y$を$x$に関して解く

$$
\large
\begin{align}
y &= Ux \\
\left(\begin{array}{c} \vdots \\ y_{n-2} \\ y_{n-1} \ y_{n} \end{array} \right) &= \left(\begin{array}{cccc} \ddots & \vdots & \vdots & \vdots \\ \cdots & u_{n-2,n-2} & u_{n-2,n-1} & u_{n-2,n} \\ \cdots & 0 & u_{n-1,n-1} & u_{n-1,n} \\ \cdots & 0 & 0 & u_{nn} \end{array} \right) \left(\begin{array}{c} \vdots \\ x_{n-2} \\ x_{n-1} \ x_{n} \end{array} \right) \\
&= \left(\begin{array}{c} \vdots \\ u_{n-2,n-2}x_{n-2} + u_{n-2,n-1}x_{n-1} + u_{n-2,n}x_{n} \\ u_{n-1,n-1}x_{n-1} + u_{n-1,n}x_{n} \\ u_{nn}x_{n} \end{array} \right)
\end{align}
$$

上記を元に$i=n$と$i<n$に場合分けし、$y_{i}$に関して方程式を解くと下記のように表せる。
・$i=n$
$$
\large
\begin{align}
y_{n} &= u_{nn}x_{n} \\
x_{1} &= \frac{y_{n}}{u_{nn}}
\end{align}
$$

・$i<n$
$$
\large
\begin{align}
y_{i} &= u_{i,n}x_{i} + \sum_{k=1}^{n-i} u_{i,n-k+1}x_{n-k+1} \\
x_{i} &= \frac{\displaystyle y_{i} – \sum_{k=1}^{n-i} u_{i,n-k+1}x_{n-k+1}}{u_{i,n}}
\end{align}
$$

上記を$i=n$から$1$まで降順に計算することでベクトル$x$の各要素の値が得られる。

LU分解を用いた逆行列の導出

$n$次単位行列$E_{n}$を考え、その$i$列を$e_{i}$とおくとき、前項のようにLU分解などを用いることで$Ax_{i}=e_{i}$を解くことで$x_{i}$に関して解ける。ここで$e_{1},…,e_{n}$に関してそれぞれ$x_{1},…,x_{n}$を解くとこの結果は下記のようにまとめることができる。
$$
\large
\begin{align}
A \left(\begin{array}{ccc} x_{1} & \cdots & x_{n} \end{array} \right) &= \left(\begin{array}{ccc} e_{1} & \cdots & e_{n} \end{array} \right) \\
AX &= E_{n}
\end{align}
$$

ここで上記の両辺の左から$A^{-1}$をかけることで$A^{-1}=X$が得られる。よって、$b = e_{1},…,e_{n}$に対して前項の連立方程式の解法を用いることで逆行列が得られる。

Pythonを用いたプログラムの作成

連立方程式の解の計算

前節の「LU分解を用いた連立方程式の解法」の内容に基づいて、下記のように連立方程式を解くことができる。

def solve_equation(A,b):
    n = A.shape[0]
    L, U = LU_decomposition(A)
    x, y = np.zeros(n), np.zeros(n)
    for i in range(n):
        if i==0:
            y[i] = b[i]
        if i>0:
            y[i] = b[i] - np.dot(L[i,:i],y[:i])
    for i in range(n):
        if i==0:
            x[n-1-i] = y[n-1-i]/U[n-1-i,n-1-i]
        if i>0:
            x[n-1-i] = (y[n-1-i] - np.dot(U[n-1-i,(n-i):],x[(n-i):]))/U[n-1-i,n-1-i]
    return x, y

A1 = np.array([[2., 1.], [1., 2.]])
b1 = np.array([3., 3.])

x, _ = solve_equation(A1,b1)
print(x)

・実行結果

> print(x)
[ 1.  1.]

$2x+y=3, x+2y=3$の連立方程式の解が$x=1,y=1$であることから計算結果は妥当であると考えることができる。

逆行列の計算

前節の「LU分解を用いた逆行列の導出」の内容に基づいて、下記のように逆行列を計算することができる。

def calc_inv_mat(A):
    n = A.shape[0]
    E = np.eye(n)
    A_inv = np.zeros([n,n])
    for i in range(n):
        A_inv[:,i], _ = solve_equation(A,E[:,i])
    return A_inv

A1 = np.array([[2., 1.], [1., 2.]])
A1_inv = calc_inv_mat(A1)
print(A1_inv)

・実行結果

> print(A1_inv)
[[ 0.66666667 -0.33333333]
 [-0.33333333  0.66666667]]

$2$次正方行列の逆行列の公式を用いることで上記の計算結果が妥当であることが確かめられる。下記では同様に$3$次正方行列の逆行列の計算結果の確認も行なった。

A2 = np.array([[2., 1., 5.], [1., 6., 2.], [7., 2., 1.]])
A2_inv = calc_inv_mat(A2)

print(A2_inv)
print(np.dot(A2,A2_inv))
print("error: {:.3f}".format(np.linalg.norm(np.dot(A2,A2_inv)-np.eye(3))))

・実行結果

> print(A2_inv)
[[-0.01092896 -0.04918033  0.15300546]
 [-0.07103825  0.18032787 -0.00546448]
 [ 0.21857923 -0.01639344 -0.06010929]]
> print(np.dot(A2,A2_inv))
[[  1.00000000e+00  -6.93889390e-18   0.00000000e+00]
 [  1.11022302e-16   1.00000000e+00   0.00000000e+00]
 [  6.93889390e-16  -9.02056208e-17   1.00000000e+00]]
> print("error: {:.3f}".format(np.linalg.norm(np.dot(A2,A2_inv)-np.eye(3))))
error: 0.000

上記の誤差の計算にはフロベニウスノルムを用いたが、$0$に近い値が得られることから結果が概ね妥当であることが確認できる。