可逆行列のLU分解
ビルドアップ
行列 $A \in \mathbb{R}^{m \times m}$ が左から掛けられた時 $(i, j)$ の成分を $0$ にする行列 $E_{ij}$ を $A$ の $ij$-消去オペレーターと定義しよう。
具体的に 正方行列 $(a_{ij}) \in \mathbb{R}^{m \times m}$ に対する $E_{ij}$ は、対角成分が $1$ であり、$(i,j)$-成分が $\displaystyle -m_{ij} = -{{a_{ij}} \over {a_{jj}}}$、残りの成分が $0$ で得られる。これは連立方程式の解法で同じ変数同士の係数を合わせて消去する操作を行列で表したものである。このような消去は連立方程式の解がただ一つだけ存在する時に意味があるので、我々が扱う $A$ は可逆行列とする。
例
簡単な例として、 $$ B := \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix} $$ とすると、$E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$ が得られる。そして $B = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$ なので $B$ に対する $(2,1)$-消去オペレーターは $$ E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ -m_{ij} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -{{1} \over {2}} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ で求めることができる。一方でこのような操作を何度も行った $$ U_{B} : = E_{32} E_{31} E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 2 \end{bmatrix} $$ は上三角行列であり、$E_{21} , E_{31} , E_{32}$ は全て下三角行列である。下三角行列同士の積は下三角行列なので $L_{B} : = E_{21} E_{31} E_{32}$ も下三角行列であり、$L^{-1} B = U$ で表せる。再び $B$ について整理すると $B = L_{B} U_{B}$ すなわち下lower三角行列と上upper三角行列の積に分解された形になる。
このような分解方法はあまりにも単純なため、$A \in \mathbb{R}^{n \times n}$ の場合についても特に問題なく一般化が可能だ。連立方程式 $A \mathbf{x} = \mathbf{b}$ の解法を考えると、$A$ が $A = LU$ に分解されるので $$ \begin{align*} LU \mathbf{x} =& \mathbf{b} \\ \mathbf{y} :=& U \mathbf{x} \end{align*} $$ と置くことができる。ここで $U$ が上三角行列なので、$\mathbf{y}$ は後退代入法を通じて容易に求めることができ、 $$ \begin{align*} L (U \mathbf{x} ) =& L \\ \mathbf{y} =& \mathbf{b} \end{align*} $$ では $L$ が下三角行列なので、$\mathbf{b}$ は前進代入法を通じて容易に求めることができる。従って、LU分解を実際に求めることができるということは、非常に多くの問題を非常に簡単かつ迅速に解決することができるという意味になる。
導出
今、直接的に $L$ と $U$ を探し出すために、$A := L U$ を満たす下三角行列 $L$ と上三角行列 $U$ を次のように定義しよう。 $$ \begin{align*} L : =& \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & 1 \end{bmatrix} \\ U : =& \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{bmatrix}\end{align*} $$ ここで、$L$ は消去オペレーターたちの積なので、対角成分が全て $1$ であることは容易に確認できる。この定義に従って、直接 $A = (a_{ij})$ の成分 $a_{ij}$ を計算してみると $$ a_{ij} = \sum_{s = 1}^{n} l_{is} u_{sj} = \sum_{s = 1}^{min(i , j)} l_{is} u_{sj} $$ となる。今、$a_{ij}$ を対角線上にある場合、対角線の上にある場合、対角線の下にある場合に分けて考えてみよう。
$i = j = k$ ならば、
$$ \begin{align*} a_{kk} =& \sum_{s = 1}^{\min(i, j)} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + l_{kk} u_{kk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + u_{kk} \end{align*} $$
これを $u_{kk}$ に対して整理すると、
$$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$
が得られる。同様に、$k < j$ なら $$ a_{kj} = \sum_{s = 1}^{k-1} l_{ks} u_{sj} + u_{kj} $$ を、$k < i$ なら $$ a_{ik} = \sum_{s = 1}^{k-1} l_{is} u_{sk} + l_{ik} u_{kk} $$ を得ることができる。これをそれぞれ $u_{kj}$ と $l_{ik}$ に対して整理すると次のようになる。 $$ \begin{align*} u_{kj} =& a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} \\ l_{ik} =& {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} \end{align*} $$
■
アルゴリズム
$(a_{ij}) \in \mathbb{R}^{n \times n}$ が可逆行列だとしよう。
ステップ 1. $k = 1$
$u_{1j} = a_{1j}$ を代入して、$\displaystyle l_{i1} = {{1} \over {u_{11}}} a_{i1}$ を計算する。
ステップ 2. $k = 2, 3, \cdots , n-1$
- ステップ 2-1. 次を計算する。 $$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$
- ステップ 2-2. $j = k+1 , k+2 , \cdots , n-1$ について次を計算する。 $$ u_{kj} = a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} $$
- ステップ 2-3. $i = k+1 , k+2 , \cdots , n-1$ について次を計算する。 $$ l_{ik} = {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} $$
ステップ 3. $k = n$ について次を計算する。 $$ u_{nn} = a_{nn} - \sum_{s = 1}^{n-1} l_{ns} u_{sn} $$
- アルゴリズムをよく見ると、$k$ が増加するにつれて、$L$ は左から右へ、$U$ は上から下へ計算されることが分かる。各ステップで必要な値は少なくともその前のステップで全て計算されているため、各ステップでの計算単位は、一行または一列になる。
コード
さて、上で例として挙げた行列 $B$ を実際に分解してみよう。
以下のコードは、前述のアルゴリズムを実際にRで実装したものである。
B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
LU<-function(A)
{
n=length(A[,1])
L<-diag(1,n)
U<-diag(0,n)
k=1
j=1
i=1
U[1,]<-A[1,]
L[,1]<-A[,1]/A[1,1]
for(k in 2:(n-1))
{
U[k,k] = A[k,k] - sum(L[k,1:(k-1)]*U[1:(k-1),k])
for(j in (k+1):n)
{
U[k,j] = A[k,j] - sum(L[k,1:(k-1)]*U[1:(k-1),j])
}
for(i in (k+1):n)
{
L[i,k] = ( A[i,k] - sum(L[i,1:(k-1)]*U[1:(k-1),k]) ) / U[k,k]
}
}
U[n,n] = A[n,n] - sum(L[n,1:(n-1)]*U[1:(n-1),n])
return(list(L=L,U=U))
}
LU(B)
require(Matrix)
expand(lu(B))
実行結果は以下の通りである。
R自体で実装されたLU分解と正確に一致することが確認できる。