LU分解

LU分解#

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

\begin{array}{r} A = L U = [\begin{array}{cccc} 1 \\ l_{21} & 1 \\ ⋮ & ⋱ & ⋱ \\ l_{n 1} & \dots & l_{n, n - 1} & 1 \end{array}] [\begin{array}{cccc} u_{11} & u_{12} & \dots & u_{1 n} \\ u_{22} & \dots & u_{2 n} \\ ⋱ & ⋮ \\ u_{n n} \end{array}] \end{array}

をLU分解と呼ぶ。

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

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

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

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

様々な $b$ というのは、たとえば $B = (b_{1}, b_{2}, \dots, b_{n})$ として行列の連立一次方程式

A X = B

を解くことができる。 $B = I$ （ $b_{1} = e_{1}, b_{2} = e_{2}, \dots, b_{n} = e_{n}$ ）とすれば $A X = I$ を解いて $X = A^{- 1}$ を求めることができる。

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

例：行列式#

\begin{array}{r} \begin{aligned} | A | & = | L U | \\ = | L | | U | (∵ 行列式の性質 | A B | = | A | | B | より) \\ = 1 \times | U | (∵ L は単位下三角行列のため) \\ = | U | \end{aligned} \end{array}

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

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

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

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

$L U$ の積を計算すると

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

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

\begin{array}{r} \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} \end{array}

これらを解いていく

一般化#

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

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

\begin{array}{r} L = (l_{1}, \dots, l_{s}), U = (\begin{array}{c} u_{1}^{T} \\ ⋮ \\ u_{s}^{T} \end{array}) \end{array}

とすると

A = l_{1} u_{1}^{T} + \dots + l_{s} u_{s}^{T}

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

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

\begin{array}{r} l_{1} = \frac{1}{a_{11} (1)} (\begin{array}{c} a_{11} (1) \\ ⋮ \\ a_{m 1} (1) \end{array}), u_{1}^{T} = (a_{11} (1), \dots, a_{1 n} (1)) \end{array}

と置くと、

\begin{array}{r} A (1) - l_{1} u_{1}^{T} = (\begin{array}{cccc} 0 & 0 & \dots & 0 \\ 0 \\ ⋮ & A (2) \\ 0 \end{array}), すなわち A = l_{1} u_{1}^{T} + (\begin{array}{cccc} 0 & 0 & \dots & 0 \\ 0 \\ ⋮ & A (2) \\ 0 \end{array}) \end{array}

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

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

\begin{array}{r} l (2) = \frac{1}{a_{11} (2)} (\begin{array}{c} a_{11} (2) \\ ⋮ \\ a_{m^{'} 1} (2) \end{array}), u (2)^{T} = (a_{11} (2), \dots, a_{1 n^{'}} (2)) \end{array}

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

\begin{array}{r} A (2) - l (2) u (2)^{T} = (\begin{array}{cccc} 0 & 0 & \dots & 0 \\ 0 \\ ⋮ & A (3) \\ 0 \end{array}) \end{array}

となる。ここで

\begin{array}{r} l_{2} = (\begin{array}{c} 0 \\ l (2) \end{array}), u_{2}^{T} = (0, u (2)^{T}) \end{array}

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

\begin{array}{r} (\begin{array}{ccccc} 0 & 0 & \dots & \dots & 0 \\ 0 \\ ⋮ & A (2) \\ ⋮ \\ 0 \end{array}) = l_{2} u_{2}^{T} + (\begin{array}{ccccc} 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \dots & 0 \\ 0 & 0 \\ ⋮ & ⋮ & A (3) \\ 0 & 0 \end{array}) \end{array}

から $ $A = l_{1} u_{1}^{T} + l_{2} u_{2}^{T} + (\begin{array}{ccccc} 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \dots & 0 \\ 0 & 0 \\ ⋮ & ⋮ & A (3) \\ 0 & 0 \end{array})$ $

となる。

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

A = l_{1} u_{1}^{T} + \dots + l_{s} u_{s}^{T} + O

そうしたら

\begin{array}{r} L = (l_{1}, \dots, l_{s}), U = (\begin{array}{c} u_{1}^{T} \\ ⋮ \\ u_{s}^{T} \end{array}) \end{array}

のように行列 $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分解で解く#

連立一次方程式

A x = b

に対し、LU分解すると

A x = L U x = b

であり、変数 $y$ を

U x = y

とおき代入すると

L y = b

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

$L y = b$ （前進代入）は

\begin{array}{r} \begin{aligned} y_{1} := b_{1} \\ y_{2} := b_{2} - l_{21} y_{1} \\ ⋮ \\ 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} \end{array}

$U x = y$ （後退代入）は

\begin{array}{r} \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}} \\ ⋮ \\ 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} \end{array}

という計算

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}

あるいは $A x = 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]])

$A x = 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$ つきの分解は $P A = L U$ となる。

$A x = b$ は2つの方程式 $L y = P b, U x = y$ と同値。

L (U x) = P^{- 1} b

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

$L y = b$ （前進代入）は

y_{i} = b_{i} - \sum_{j = 1}^{i - 1} l_{i j} y_{j}

この $y$ は枢軸選択つきの前進消去で得られる ${(b_{1}^{(1)}, \dots, b_{n}^{(n)})}^{⊤}$ に等しい

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

x_{i} = \frac{1}{u_{i i}} (b_{i}^{(i)} - \sum_{j = i + 1}^{n} u_{i j} x_{j}) (i = n, n - 1, \dots, 1)

で解ける

# データの作成
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}, \dots, e_{n})$ として

A X = I

を解く

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
一般化した計算手順：『プログラミングのための線形代数』