LU分解 - データサイエンス関連+αのメモ

与えられた行列 $A$ を単位下三角行列 $L$ と上三角行列 $U$ との積で表すこと

A = LU = \left[\begin{array}{cccc} 1 & & & \\ l_{21} & 1 & & \\ \vdots & \ddots & \ddots & \\ l_{n 1} & \cdots & l_{n, n-1} & 1 \end{array}\right]\left[\begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1 n} \\ & u_{22} & \cdots & u_{2 n} \\ & & \ddots & \vdots \\ & & & u_{n n} \end{array}\right]

(1)

をLU分解と呼ぶ。

分解すると何が嬉しいのか？¶

行列式や一次方程式を少ない計算量で求められる材料として $L,U$ が手に入り、全体としての計算を効率化できる可能性がある。

例：連立1次方程式・逆行列¶

例えば様々な $b$ に対して $Ax=b$ を解く場合は $A=LU$ と一旦分解して、後で $L,U$ を使い回せばよい

様々な $b$ というのは、たとえば $B=(b_1, b_2, \cdots, b_n)$ として行列の連立一次方程式

AX = B

(2)

を解くことができる。 $B=I$ （ $b_1 = e_1, b_2 = e_2, \cdots, b_n = e_n$ ）とすれば $AX=I$ を解いて $X=A^{-1}$ を求めることができる。

計算量はLU分解に $O(n^3/3)$ かかる。そこから逆行列を計算する場合は $O(2n^3/3)$ かかり、全体で $O(n^3/3 + 2n^3/3)=O(n^3)$ となる。これは掃き出し法と同程度だが、もとの行列が疎行列だとLUも疎行列になるので、入力が疎行列の場合も含めた計算量でいうとLU分解のほうが効率的らしい。

例：行列式¶

\begin{align} |A| &= |LU|\\ &= |L| |U| \quad (\because \text{ 行列式の性質 } |AB|=|A||B| \text{ より })\\ &= 1\times |U| \quad (\because L \text{ は単位下三角行列のため })\\ &= |U| \end{align}

(3)

三角行列の行列式は対角行列の積になるため、 $|L|=1$ であり $|A|=|U|$ になるし $|U|$ も $U$ の対角成分の積をとるだけでいい

計算手順（ドゥーリトル法）¶

$L$ の対角成分は1に固定し、次のような分解を仮定する

A = \left[\begin{array}{cc} 1 & 0 \\ l_{21} & 1 \end{array}\right] \left[\begin{array}{cc} u_{11} & u_{12} \\ 0 & u_{22} \end{array}\right]

(4)

$LU$ の積を計算すると

A=\left[\begin{array}{ccc} u_{11} & u_{12}\\ l_{21} u_{11} & l_{21} u_{12} + u_{22} \end{array}\right]

(5)

となるため、 $A$ の各要素との対応がとれる

\begin{aligned} & a_{11}=u_{11} \\ & a_{12}=u_{12} \\ & a_{21}=l_{21} u_{11} \\ & a_{22}=l_{21} u_{12}+u_{22} \end{aligned}

(6)

これらを解いていく

一般化¶

もっと一般化すると次のようになる

$m\times n$ 行列 $A=(a_{ij})$ に対して、 $s=\min(m,n)$ とおき、 $m\times s$ の $L$ と $s\times n$ の $U$ で分解することを考える。

L = \left(\boldsymbol{l}_1, \ldots, \boldsymbol{l}_s\right), \quad U = \left(\begin{array}{c} \boldsymbol{u}_1^T \\ \vdots \\ \boldsymbol{u}_s^T \end{array}\right)

(7)

とすると

A = l_1 \boldsymbol{u}_1^T+\cdots+\boldsymbol{l}_s \boldsymbol{u}_s^T

(8)

とできるので、この表記を使っていく。

反復計算するので $A = A(1) = (a_{ij}(1))$ と表記する

\boldsymbol{l}_1=\frac{1}{a_{11}(1)}\left(\begin{array}{c} a_{11}(1) \\ \vdots \\ a_{m 1}(1) \end{array}\right), \quad \boldsymbol{u}_1^T=\left(a_{11}(1), \ldots, a_{1 n}(1)\right)

(9)

と置くと、

A(1)-\boldsymbol{l}_1 \boldsymbol{u}_1^T=\left(\begin{array}{c|ccc} 0 & 0 & \cdots & 0 \\ \hline 0 & & & \\ \vdots & & A(2) & \\ 0 & & & \end{array}\right), \quad \text { すなわち } A=l_1 \boldsymbol{u}_1^T+\left(\begin{array}{c|ccc} 0 & 0 & \cdots & 0 \\ \hline 0 & & & \\ \vdots & & A(2) & \\ 0 & & \end{array}\right)

(10)

となる。残ったブロック行列を $A(2)$ としている。

ゼロじゃない要素を $\boldsymbol{l}(2), \boldsymbol{u}(2)$ とおき、同様に

\boldsymbol{l}(2)=\frac{1}{a_{11}(2)}\left(\begin{array}{c} a_{11}(2) \\ \vdots \\ a_{m^{\prime} 1}(2) \end{array}\right), \quad \boldsymbol{u}(2)^T=\left(a_{11}(2), \ldots, a_{1 n^{\prime}}(2)\right)

(11)

とすると（ここで $m'=m-1, n'=n-1$ ）、同様に

A(2)-l(2) \boldsymbol{u}(2)^T=\left(\begin{array}{c|ccc} 0 & 0 & \cdots & 0 \\ \hline 0 & & & \\ \vdots & & A(3) & \\ 0 & & & \end{array}\right)

(12)

となる。ここで

\boldsymbol{l}_2=\left(\begin{array}{c} 0 \\ \boldsymbol{l}(2) \end{array}\right), \quad \boldsymbol{u}_2^T=\left(0, \boldsymbol{u}(2)^T\right)

(13)

のように頭に0をつけてサイズを戻せば

\left(\begin{array}{c|cccc} 0 & 0 & \cdots & \cdots & 0 \\ \hline 0 & & & & \\ \vdots & & A(2) & & \\ \vdots & & & & \\ 0 & & & & \end{array}\right)=l_2 \boldsymbol{u}_2^T+\left(\begin{array}{cc|ccc} 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & \cdots & 0 \\ \hline 0 & 0 & & & \\ \vdots & \vdots & & A(3) & \\ 0 & 0 & & \end{array}\right)

(14)

から

A=l_1 \boldsymbol{u}_1^T+\boldsymbol{l}_2 \boldsymbol{u}_2^T+\left(\begin{array}{cc|ccc} 0 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 0 & \cdots & 0 \\ \hline 0 & 0 & & & \\ \vdots & \vdots & & A(3) & \\ 0 & 0 & & & \end{array}\right)

(15)

となる。

これを何度も繰り返していくと、最後は残る行列がゼロ行列になる

A=l_1 \boldsymbol{u}_1^T+\cdots+\boldsymbol{l}_s \boldsymbol{u}_s^T+O

(16)

そうしたら

L = \left(\boldsymbol{l}_1, \ldots, \boldsymbol{l}_s\right), \quad U = \left(\begin{array}{c} \boldsymbol{u}_1^T \\ \vdots \\ \boldsymbol{u}_s^T \end{array}\right)

(17)

のように行列 $L,U$ を作ればよい

pythonの例¶

import numpy as np

# 対象の行列
A = np.array([
    [1, 2],
    [3, 4]
])

# iter 1
l1 = (1 / A[0,0]) * A[:, 0]
u1 = A[0, :]
print(f"{l1=}, {u1=}")

l1=array([1., 3.]), u1=array([1, 2])

# ベクトルの外積
np.outer(l1, u1)

array([[1., 2.],
       [3., 6.]])

A2 = A - np.outer(l1, u1)
A2

array([[ 0.,  0.],
       [ 0., -2.]])

A_ = A2[1:, 1:]

l2_ = (1 / A_[0,0]) * A_[:, 0]
u2_ = A_[0, :]
print(f"{l2_=}, {u2_=}")

l2_=array([1.]), u2_=array([-2.])

# 残りがゼロに
A3 = A_ - np.outer(l2_, u2_)
A3

array([[0.]])

# 先頭にゼロを足して次元を合わせる
l2 = np.concatenate([np.zeros(shape=(1,)), l2_])
u2 = np.concatenate([np.zeros(shape=(1,)), u2_])

# それぞれのiterationの外積の和が元の行列になる
np.outer(l1, u1) + np.outer(l2, u2)

array([[1., 2.],
       [3., 4.]])

L = np.array([l1, l2]).T
L

array([[1., 0.],
       [3., 1.]])

U = np.array([u1, u2])
U

array([[ 1.,  2.],
       [ 0., -2.]])

L @ U

array([[1., 2.],
       [3., 4.]])

行列式をLU分解で求める¶

print(f"""
|A| = {np.linalg.det(A):.1f}
|L||U| = {np.linalg.det(L) * np.linalg.det(U):.1f}
Uの対角成分の積 = {np.diag(U)[0] * np.diag(U)[1]}
""")


|A| = -2.0
|L||U| = -2.0
Uの対角成分の積 = -2.0

連立一次方程式をLU分解で解く¶

連立一次方程式

Ax = b

(18)

に対し、LU分解すると

Ax = LUx = b

(19)

であり、変数 $y$ を

Ux=y

(20)

とおき代入すると

Ly = b

(21)

になることから変数 $y$ を求める。

$Ly = b$ （前進代入）は

\begin{aligned} & y_1:=b_1 \\ & y_2:=b_2-l_{21} y_1 \\ & \vdots \\ & y_{n-1}:=b_{n-1}-\sum_{j=1}^{n-2} l_{n-1, j} y_j \\ & y_n:=b_n-\sum_{j=1}^{n-1} l_{n j} y_j \end{aligned}

(22)

$Ux = y$ （後退代入）は

\begin{aligned} & x_n:=\frac{y_n}{u_{n n}} \\ & x_{n-1}:=\frac{y_{n-1}-u_{n-1, n} x_n}{u_{n-1, n-1}}\\ & \vdots\\ & x_2:=\frac{y_2-\sum_{j=3}^n u_{2 j} x_j}{u_{22}} \\ & x_1:=\frac{y_1-\sum_{j=2}^n u_{1 j} x_j}{u_{11}} \end{aligned}

(23)

という計算

def forward_substitution(L, b):
    # 前進代入
    n = L.shape[0]
    y_ = np.zeros_like(b, dtype=np.double)
    for i in range(n):
        sumj = 0
        for j in range(i):
            sumj += L[i, j] * y_[j]
        y_[i] = b[i] - sumj
    return y_

def backword_substitution(U, y):
    # 後退代入
    n = U.shape[0]
    x = np.zeros_like(y, dtype=np.double)
    for i in range(n-1, -1, -1):
        sumj = 0
        for j in range(n-1, i-1, -1):
            sumj += U[i, j] * x[j]
        x[i] = (y[i] - sumj) / U[i, i]
    return x

xGETRS¶

前進代入・後退代入を行うのがLAPACKのxGETRSサブルーチン

b = np.array([1, 2])

from scipy.linalg.lapack import dgetrf, dgetrs
LU, P, info = dgetrf(A)
dgetrs(LU, P, b)

(array([0. , 0.5]), 0)

y = np.linalg.inv(L) @ b
y

array([ 1., -1.])

x = np.linalg.inv(U) @ y
x

array([0. , 0.5])

# estimate b
A @ x

array([1., 2.])

逆行列¶

A^{-1} = U^{-1} L^{-1}

(24)

あるいは $Ax=b$ の $x$ を求める操作を $n$ 回繰り返す

np.linalg.inv(A)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

np.linalg.inv(U) @ np.linalg.inv(L)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

$Ax=b$ の $x$ を求める操作を $n$ 回繰り返す

y = L @ b
b = np.array([1, 2])

scipy.linalg.lu¶

scipy.linalg.lu — SciPy v1.12.0 Manual

from scipy.linalg import lu

P, L, U = lu(A)

array([[0., 1.],
       [1., 0.]])

array([[1.        , 0.        ],
       [0.33333333, 1.        ]])

array([[3.        , 4.        ],
       [0.        , 0.66666667]])

P @ L @ U

array([[1., 2.],
       [3., 4.]])

置換行列 $P$ つきの分解は $PA=LU$ となる。

$Ax = b$ は2つの方程式 $Ly=Pb, Ux=y$ と同値。

L(Ux) = P^{-1} b

(25)

となり、並べ替えをもとに戻した $b = P^{-1} b$ を使ってやればいい

$Ly = b$ （前進代入）は

y_i = b_i - \sum^{i-1}_{j=1} l_{ij} y_j

(26)

この $y$ は枢軸選択つきの前進消去で得られる $\left(b_1^{(1)}, \cdots, b_n^{(n)}\right)^{\top}$ に等しい

$Ux = y$ （後退代入）は消去法と同じく

x_i=\frac{1}{u_{i i}}\left(b_i^{(i)}-\sum_{j=i+1}^n u_{i j} x_j\right) \quad(i=n, n-1, \cdots, 1)

(27)

で解ける

# データの作成
import numpy as np

A = np.array([
    [1,2,3],
    [4,5,6],
    [7,8,0],
])
x_true = np.array([1, 2, 3])
b = A @ x_true
b

array([14, 32, 23])

from scipy.linalg import lu
P, L, U = lu(A)

# b := P^{-1} b
b_pi = np.linalg.inv(P) @ b
b_pi

array([23., 14., 32.])

def forward_substitution(L, b):
    # 前進代入
    n = L.shape[0]
    y_ = np.zeros_like(b, dtype=np.double)
    for i in range(n):
        sumj = 0
        for j in range(i):
            sumj += L[i, j] * y_[j]
        y_[i] = b[i] - sumj
    return y_

def backword_substitution(U, y):
    # 後退代入
    n = U.shape[0]
    x = np.zeros_like(y, dtype=np.double)
    for i in range(n-1, -1, -1):
        sumj = 0
        for j in range(n-1, i-1, -1):
            sumj += U[i, j] * x[j]
        x[i] = (y[i] - sumj) / U[i, i]
    return x

y_ = forward_substitution(L, b_pi)
y_

array([23. , 10.71428571, 13.5 ])

x = backword_substitution(U, y_)
x

array([1., 2., 3.])

逆行列を解く¶

$I=(e_1, \cdots, e_n)$ として

AX=I

(28)

を解く

n = A.shape[0]
I = np.eye(n)
I

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

X = np.zeros_like(I)
for i in range(n):
    b = I[:, i]
    b_pi = np.linalg.inv(P) @ b
    y_ = forward_substitution(L, b_pi)
    x = backword_substitution(U, y_)
    X[:, i] = x

X.round(3)

array([[-1.778,  0.889, -0.111],
       [ 1.556, -0.778,  0.222],
       [-0.111,  0.222, -0.111]])

np.linalg.inv(A).round(3)

array([[-1.778,  0.889, -0.111],
       [ 1.556, -0.778,  0.222],
       [-0.111,  0.222, -0.111]])

参考文献¶

ドゥーリトル法：LU分解 - Wikipedia
一般化した計算手順：『プログラミングのための線形代数』