LU分解#
与えられた行列\(A\)を単位下三角行列\(L\)と上三角行列\(U\)との積で表すこと
をLU分解と呼ぶ。
分解すると何が嬉しいのか？#
行列式や一次方程式を少ない計算量で求められる材料として\(L,U\)が手に入り、全体としての計算を効率化できる可能性がある。
例：連立1次方程式・逆行列#
例えば様々な\(b\)に対して\(Ax=b\)を解く場合は\(A=LU\)と一旦分解して、後で\(L,U\)を使い回せばよい
様々な\(b\)というのは、たとえば\(B=(b_1, b_2, \cdots, b_n)\)として行列の連立一次方程式
を解くことができる。\(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分解のほうが効率的らしい。
例：行列式#
三角行列の行列式は対角行列の積になるため、\(|L|=1\)であり\(|A|=|U|\)になるし\(|U|\)も\(U\)の対角成分の積をとるだけでいい
計算手順（ドゥーリトル法）#
\(L\)の対角成分は1に固定し、次のような分解を仮定する
\(LU\)の積を計算すると
となるため、\(A\)の各要素との対応がとれる
これらを解いていく
一般化#
もっと一般化すると次のようになる
\(m\times n\)行列\(A=(a_{ij})\)に対して、\(s=\min(m,n)\)とおき、\(m\times s\)の\(L\)と\(s\times n\)の\(U\)で分解することを考える。
とすると
とできるので、この表記を使っていく。
反復計算するので\(A = A(1) = (a_{ij}(1))\)と表記する
と置くと、
となる。残ったブロック行列を\(A(2)\)としている。
ゼロじゃない要素を\(\boldsymbol{l}(2), \boldsymbol{u}(2)\)とおき、同様に
とすると（ここで\(m'=m-1, n'=n-1\)）、同様に
となる。ここで
のように頭に0をつけてサイズを戻せば
から $\( 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) \)$
となる。
これを何度も繰り返していくと、最後は残る行列がゼロ行列になる
そうしたら
のように行列\(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分解で解く#
連立一次方程式
に対し、LU分解すると
であり、変数\(y\)を
とおき代入すると
になることから変数\(y\)を求める。
\(Ly = b\)（前進代入）は
\(Ux = y\)（後退代入）は
という計算
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.])
逆行列#
あるいは\(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)
P
array([[0., 1.],
[1., 0.]])
L
array([[1. , 0. ],
[0.33333333, 1. ]])
U
array([[3. , 4. ],
[0. , 0.66666667]])
P @ L @ U
array([[1., 2.],
[3., 4.]])
置換行列\(P\)つきの分解は\(PA=LU\)となる。
\(Ax = b\)は2つの方程式\(Ly=Pb, Ux=y\)と同値。
となり、並べ替えをもとに戻した\(b = P^{-1} b\)を使ってやればいい
\(Ly = b\)（前進代入）は
この\(y\)は枢軸選択つきの前進消去で得られる\(\left(b_1^{(1)}, \cdots, b_n^{(n)}\right)^{\top}\)に等しい
\(Ux = y\)（後退代入）は消去法と同じく
で解ける
# データの作成
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)\)として
を解く
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
一般化した計算手順：『プログラミングのための線形代数』