LU分解

数万円のノート PC があれば、ガウス・ジョルダンの消去法 で 1000 変数の連立方程式を解くことができます。しかし、現代の複雑な科学計算では 10,000 あるいは 100,000 を超える変数の連立方程式を解くことを要求されることも珍しくありません。このような巨大連立方程式に備えて少しでも計算コストを下げようと、様々な高速アルゴリズムが考案されています。SciPy には PA=LU 分解とよばれるアルゴリズムで連立方程式を解く linalg.lu_solve() が実装されています。

LU分解

行列を下三角行列 (lower triangular matrix) と上三角行列 (upper triangular matrix) の積の形で表すことを LU分解 とよびます。
 \[\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}=\begin{bmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\end{bmatrix}\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}\tag{1}\]
サイズ $3$ の正方行列
 \[A=\begin{bmatrix}2&2&3\\2&3&4\\1&1&2\end{bmatrix}\tag{2}\]
を使って、LU分解の具体的な手順を見ていきます。最初に $A$ に操作をほどこして上三角行列 $U$ を作り出すことを目指します。分解しやすい行列を選んでいるので、この作業は $2$ 段階で終わります。

操作①　$2$ 行目から $1$ 行目を引く
 \[\begin{bmatrix}2&2&3\\0&1&1\\1&1&2\end{bmatrix}\tag{3}\]
操作②　$3$ 行目から $1$ 行目の $1/2$ 倍を引く
 \[\begin{bmatrix}2&2&3\\0&1&1\\0&0&\frac{1}{2}\\ \end{bmatrix}\equiv U\tag{4}\]
上三角行列が得られたので、これを $U$ とおいて操作を止めます。
① と ② は、それぞれ 基本行列 (下三角行列)
 \[E_{2,1,-1}=\begin{bmatrix}1&0&0\\-1&1&0\\0&0&1\end{bmatrix},\quad E_{3,1,-1/2}=\begin{bmatrix}1&0&0\\0&1&0\\-\frac{1}{2}&0&1\end{bmatrix}\tag{5}\]
を乗じることに等しい操作なので、
 \[E_{3,1,-1/2\ }E_{2,1,-1\ }A=U\tag{6}\]
と表せます。両辺に左側から $E_{3,1,-1/2\ }$ と $E_{2,1,-1\ }$ の逆行列を掛けると、
 \[A=E_{2,1,-1}^{-1\ }E_{3,1,-1/2}^{-1\ }U\tag{7}\]
となるので、$E_{2,1,-1}^{-1\ }E_{3,1,-1/2}^{-1}$ が下三角行列 $L$ となります (「下三角行列と上三角行列」の記事で説明したように、三角行列の逆行列や三角行列同士の積は三角行列です)。$E_{i,j,k}$ の逆行列は $E_{i,j,-k}$ なので、
 \[L=E_{2,1,1\ }E_{3,1,1/2}=\begin{bmatrix}1&0&0\\1&1&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&0&0\\0&1&0\\\frac{1}{2}&0&1\end{bmatrix}=\begin{bmatrix}1&0&0\\1&1&0\\\frac{1}{2}&0&1\end{bmatrix}\tag{8}\]
となります。すなわち、行列 $A$ は次のように分解できます。
 \[\begin{bmatrix}2&2&3\\2&3&4\\1&1&2\end{bmatrix}=\begin{bmatrix}1&0&0\\1&1&0\\\frac{1}{2}&0&1\end{bmatrix}\begin{bmatrix}2&2&3\\0&1&1\\0&0&\frac{1}{2}\end{bmatrix}\tag{9}\]
ただし、以上はあくまで人間が手計算で LU分解する手順を示したもので、数値計算ではより効率的なアルゴリズムを使います。

LDU分解

上の例で見たように、一般に LU分解においては $L$ の対角成分は $1$ に揃えられていますが、$U$ の対角成分は揃っていません。そこで、
 \[\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}=\begin{bmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\end{bmatrix}\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}\tag{10}\]
のように分解されたとき、$U$ を対角行列 $D$ を使って
 
\[U=\begin{bmatrix}u_{11}&0&0\\0&u_{22}&0\\0&0&u_{33}\end{bmatrix}\begin{bmatrix}1&u_{12}/u_{11}&u_{13}/u_{11}\\0&1&u_{23}/u_{22}\\0&0&1\end{bmatrix}\tag{11}\]
のように分解すれば、対角成分を $1$ に揃えた上三角行列が現れます。この三角行列をあらためて $U$ とおいて、$A=LDU$ のように表すことを LDU 分解 とよびます。(9) 式で表された LU 分解
 \[\begin{bmatrix}2&2&3\\2&3&4\\1&1&2\end{bmatrix}=\begin{bmatrix}1&0&0\\1&1&0\\\frac{1}{2}&0&1\end{bmatrix}\begin{bmatrix}2&2&3\\0&1&1\\0&0&\frac{1}{2}\end{bmatrix}\tag{9}\]
においては、右側の上三角行列を (11) の方法で分解して、
 \[\begin{bmatrix}2&2&3\\2&3&4\\1&1&2\end{bmatrix}=\begin{bmatrix}1&0&0\\1&1&0\\\frac{1}{2}&0&1\end{bmatrix}\begin{bmatrix}2&0&0\\0&1&0\\0&0&\frac{1}{2}\end{bmatrix}\begin{bmatrix}1&1&\frac{3}{2}\\0&1&1\\0&0&1\end{bmatrix}\tag{12}\]
と表すことができます。LU分解と LDU分解の間に本質的な違いはありません。

PA=LU分解

任意の行列 $A$ が与えられたとき、必ずしも $E$ を乗じるだけで上三角行列を作り出せるとは限りません。しかし、$A$ が可逆であれば、$A$ に適切な行交換をほどこした行列 $PA$ については必ず LU 分解できることが知られています。このように行列 $PA$ に対して LU分解を実行することを PA=LU分解 とよびます。$PA=LU$ の左側に $P$ の逆行列を掛ければ $A=P^{-1}LU$ の形で表されます。$P^{-1}$ も置換行列なので、これをあらためて $P$ と置き直せば、$A=PLU$ の形に分解されたことになります。

たとえば、行列
 \[A=\begin{bmatrix}0&1&1\\1&3&3\\2&5&8\end{bmatrix}\tag{13}\]
が与えられたとします。1 行目と 2 行目を交換して (置換行列 $P=P_{1,2}$ を掛けて)、
 \[PA=\begin{bmatrix}1&3&3\\0&1&1\\2&5&8\end{bmatrix}\tag{14}\]
という行列を作ります。次に 3 行目から 1 行目の 2 倍を引きます ($E_{3,1,-2}$ を掛けます)。
 \[\begin{bmatrix}1&3&3\\0&1&1\\0&-1&2\end{bmatrix}\tag{15}\]
3 行目に 2 行目を加えると ($E_{3,2,1}$ を掛けると)、上三角行列 $U$ を得ます。
 \[\begin{bmatrix}1&3&3\\0&1&1\\0&0&3\end{bmatrix}\tag{16}\]
以上の操作を行列で書き表すと
 \[E_{3,2,1\ }E_{3,1,-2\ }PA=U\tag{17}\]
となるので、$E_{3,2,1\ }E_{3,1,-2}$ の逆行列を掛けて
 \[PA=E_{3,1,2\ }E_{3,2,-1\ }U\tag{18}\]
となります。したがって、下三角行列
 \[L=E_{3,1,2\ }E_{3,2,-1}=\begin{bmatrix}1&0&0\\0&1&0\\2&0&1\end{bmatrix}\begin{bmatrix}1&0&0\\0&1&0\\0&-1&1\end{bmatrix}=\begin{bmatrix}1&0&0\\0&1&0\\2&-1&1\end{bmatrix}\tag{19}\]
を得るので、
 \[\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}0&1&1\\1&3&3\\2&5&8\end{bmatrix}=\begin{bmatrix}1&0&0\\0&1&0\\2&-1&1\end{bmatrix}\begin{bmatrix}1&3&3\\0&1&1\\0&0&3\end{bmatrix}\tag{20}\]
のように分解されます。

scipy.linalg.lu()

scipy.linalg.lu() は受け取った行列を A=PLU の形に分解します。適当な行列 A をつくって分解を試してみます。

# scipy_pa_ludecomposition

# In[1]

import numpy as np
from scipy.linalg import lu

np.random.seed(10)

# ランダムな要素をもつ行列を生成
A = np.random.randint(1, 10, (3, 3))

# PA=LU分解
P, L, U = lu(A)

print("A:\n{}\n".format(A))
print("P:\n{}\n".format(P))
print("L:\n{}\n".format(L))
print("U:\n{}".format(U))

'''
A:
[[5 1 2]
 [1 2 9]
 [1 9 7]]

P:
[[1. 0. 0.]
 [0. 0. 1.]
 [0. 1. 0.]]

L:
[[1.         0.         0.        ]
 [0.2        1.         0.        ]
 [0.2        0.20454545 1.        ]]

U:
[[5.   1.   2.  ]
 [0.   8.8  6.6 ]
 [0.   0.   7.25]]
'''

念のために PLU を計算して、A に一致することを確認しておきます。

# In[2]

# PLUを計算
print(P@L@U)

# [[5. 1. 2.]
#  [1. 2. 9.]
#  [1. 9. 7.]]

scipy.linalg.lu_factor()

scipy.linalg.lu_factor() は受け取った行列を PA=LU分解して、$L+U$ から $L$ の主対角成分を除いた行列 LU と置換行列 P を表現するピボットインデックスを返します。

# scipy_pa_ludecomposition_pivot

# In[1]

import numpy as np
from scipy.linalg import lu_factor

# 乱数シードを設定
np.random.seed(0)

# ランダムな要素をもつ行列を生成
A = np.random.randint(-4, 4, (3, 3))

# PA=LU分解
LU, piv = lu_factor(A)

print("A:\n{}".format(A))
print("LU:\n{}".format(LU))
print("piv:\n{}".format(piv))

'''
A:
[[ 0  3  1]
 [-4 -1 -1]
 [-1  3 -3]]

LU:
[[-4.         -1.         -1.        ]
 [ 0.25        3.25       -2.75      ]
 [-0.          0.92307692  3.53846154]]

piv:
[1 2 2]
'''

後述する scipy.lu_solve() に scipy.linalg.lu_factor() を渡せば連立方程式を解くことができます。

sympy.matrices.dense.MutableDenseMatrix.LUdecomposition()

SymPy の行列は LUdecomposition()メソッドを使って PA=LU 分解できます。
行列の要素に記号を設定すると、一般的な分解の形を確認できます。

# sympy_ludecomposition

# In[1]

from sympy import Symbol, var, init_printing
from sympy.matrices import Matrix

init_printing()

var("a:z")

# 2次正方行列
X = Matrix([[a, b],
            [c, d]])

# PA=LU分解
L, U, _ = X.LUdecomposition()

display(L)
display(U)

LU分解と連立方程式

記事の冒頭で述べたように、LU 分解は連立方程式 $A\boldsymbol{x}=\boldsymbol{b}$ を効率的に (ガウスジョルダンの消去法よりも速く) 解くために編み出されました。$A=LU$ とおけば、方程式は
 \[LU\boldsymbol{x}=\boldsymbol{b}\tag{21}\]
と表せます。さらに
 \[U\boldsymbol{x}=\boldsymbol{c}\tag{22}\]
とおけば、
 \[L\boldsymbol{c}=\boldsymbol{b}\tag{23}\]
となります。最初に (23) を解いて $\boldsymbol{c}$ を得て、それを (22) に代入して $\boldsymbol{x}$ を求めます。方程式を解く手順を二段階に分けたことで、一見するとかえって煩わしくなったと考えるかもしれません。しかし数値計算的には非常に巧妙な仕掛けがほどこされています。三角行列をベクトルに乗じた形の方程式は簡単に解けます。方程式 (23) を具体的に書き表せば
 \[\begin{bmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\end{bmatrix}\begin{bmatrix}c_1\\c_2\\c_3\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\b_3\end{bmatrix}\tag{24}\]
のような形になります。すなわち、
 \[\begin{align*}l_{11}c_1=b_1\tag{25}\\[6pt]l_{21}c_1+l_{22}c_2=b_2\tag{26}\\[6pt]l_{31}c_1+l_{32}c_2+l_{33}c_3=b_3\tag{27}\end{align*}\]
という連立方程式なので、(25) から $c_1$ を得て、それを (26) に代入すると $c_2$ を得ます。さらに $c_1$ と $c_2$ を (27) に代入して $c_3$ が求められます。$U\boldsymbol{x}=\boldsymbol{c}$ についても同様の手順で解けます。これはとても単純なアルゴリズムなので実装も容易であり、実質的には LU分解が済んだ時点で方程式は解けているといえます。

$A$ のサイズが $n\times n$ であるとき、ガウス・ジョルダン消去法による連立方程式の処理にかかる計算量は $O(n^3)$ です。これに対して、LU分解を用いた場合の計算量は $O(n^2)$ なので、LU分解のほうが圧倒的に高速です。
 
LU分解を用いるもう１つの利点は、分解を一度だけ行ってメモリに記憶させておけば、異なる $b$ に対して $A\boldsymbol{x}=\boldsymbol{b}$ を素早く解くことができるということです。先にも述べたように、面倒な処理は分解過程に集中していて、三角行列を含む方程式の処理は簡単だからです。

scipy.linalg.lu_solve()

scipy.linalg.lu_solve(LU, b) は PA=LU分解を使って連立方程式 LUx = b を解きます。LU には scipy.linalg.lu_factor() の戻り値を渡します。10 変数の連立方程式を解いてみましょう。

# scipy_pa_lu_solve

# In[1]

import numpy as np
from scipy.linalg import lu_factor, lu_solve

np.random.seed(0)

# サイズ10の正方行列を生成
A = np.random.randint(-6, 6, (10, 10))

# 10要素のベクトルを生成
b = np.random.randint(0, 10, 10)

# Ax=bを解く
x = lu_solve(lu_factor(A), b)

print("A:\n{}".format(A))
print("b:\n{}".format(b))
print("x:\n{}".format(x))

'''
A:
[[-1 -6 -3  5 -3  1  3 -3 -1 -4]
 [-2  1  0  2  2  4 -5  0  1  1]
 [ 2 -5 -1  3  2  3 -2 -3 -6 -3]
 [-1 -6 -4 -3  2 -5 -3 -3 -3  1]
 [-6 -5  3  3 -6  4 -2  1 -3  5]
 [-4  1 -4 -6 -6 -2 -1 -1  0  2]
 [-2 -5 -2  3  4  4  2 -5 -5  1]
 [ 3  3 -3  0  1  5 -4  5 -6 -3]
 [-1  3  4 -2  5 -2  0 -2 -2 -3]
 [-2 -2  2 -2 -3  4  1 -1 -1 -6]]

b:
[1 5 9 3 0 5 0 1 2 4]

x:
[ 0.9223642   0.59032185  0.30489475 -0.01901043 -1.05452909
0.22507008 -1.59292475 -2.22283744 -0.21934587 -0.2440726 ]
'''