対称行列のLDU分解
概要
$L^{T}$ は上三角行列なので、$A = LU$ のコンテキストにおいて、$U$ を $U:= DL^{T}$ に変換すると考えても良い。一般的なLU 分解よりも条件が厳しくなった分、計算量はかなり減る。
定理
逆行列が存在する対称行列 $A$ に対し、$A = LDL^{T}$ を満たす下三角行列 $L$ と対角行列 $D$ が存在する。
証明
$A$ は逆行列が存在するから、$A = L U$ を満たす下三角行列 $L$ と上三角行列 $U$ が存在する。一方 $A = A^{T}$ が与えられると、
$$ L U = A = A^{T} = U^{T} L^{T} $$
LU 分解では $L$ は対角成分が全て $1$ である下三角行列として求まるから、$\det L = 1$、すなわち逆行列 $L^{-1}$ が存在する。$L^{-1}$ は下三角行列であり、 $$ L U = U^{T} L^{T} $$ に $U$ を整理すると、 $$ U = L^{-1} U^{T} L^{T} $$ 両辺に $(L^{T})^{-1}$ をかけると、 $$ U (L^{T})^{-1} = L^{-1} U^{T}$$ を得る。$U (L^{T})^{-1} = L^{-1} U^{T}$ の左辺は上三角行列で、右辺は下三角行列なので、 $$ D: = U (L^{T})^{-1} = L^{-1} U^{T} $$ は対角行列であることがわかる。これを $U$ に対して整理すると $U = DL^{T}$ が得られ、以下が得られる。 $$ A = LU = LDL^{T} $$
■
計算
以下のように定義して、具体的に $L$ と $D$ を計算してみよう。
$$ 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},\quad D : = \begin{bmatrix} d_{11} & 0 & 0 & \cdots & 0 \\ 0 & d_{22} & 0 & \cdots & 0 \\ 0 & 0 & d_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & d_{nn} \end{bmatrix} $$
$A = (a_{ij}) = LDL^{T}$ を直接計算してみよう。
$l_{ij} > 0$ であり、また $l_{ii} = 1$ であるから、
$$ \begin{align*} a_{ij} =& \sum_{\nu = 1}^{n} l_{i \nu} \sum_{\mu = 1}^{n} d_{\nu \mu} l_{\mu j}^{T} \\ =& \sum_{\nu = 1}^{n} l_{i \nu} \sum_{\mu = 1}^{n} \delta_{\nu \mu} d_{\nu \mu} l_{ j \mu} \\ =& \sum_{\nu = 1}^{n} l_{i \nu} d_{\nu \nu} l_{ j \nu} \\ =& \sum_{\nu = 1}^{ \min (i,j) } l_{i \nu} d_{\nu \nu} l_{ j \nu} \end{align*} $$
まずは $j \le i$ の場合を考えてみよう。
$$ a_{ij} = \sum_{\nu = 1}^{j} l_{i \nu} d_{\nu \nu} l_{j \nu} = \sum_{\nu = 1}^{j-1} l_{i \nu} d_{\nu \nu} l_{j \nu} + l_{ij} d_{jj} l_{jj} = \sum_{\nu = 1}^{j-1} l_{i \nu} d_{\nu \nu} l_{j \nu} + l_{ij} d_{jj} $$
ここで $i = j = k$ の場合、
$$ a_{kk} = \sum_{\nu = 1}^{k-1} l_{k \nu} d_{\nu \nu} l_{k \nu} + d_{kk} d_{kk} $$
について整理すると、
$$ d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2 $$
逆に、$j = k$ であり、$i = k+1, \cdots , n$ だとすると、
$$ a_{ik} = \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} + l_{ik} d_{kk} $$
$l_{ik}$ について整理すると、
$$ l_{ik} = {{1} \over {d_{kk}}} \left\{ a_{ik} - \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} \right\} $$
アルゴリズム
$(a_{ij}) \in \mathbb{R}^{n \times n}$ が逆行列を持つ対称行列だとする。
**ステップ1. $k = 1$
$d_{11} = a_{11}$ を代入して、$\displaystyle l_{i1} = {{1} \over {d_{11}}} a_{i1}$ を計算する。
**ステップ2. $k = 2, 3, \cdots , n-1$
ステップ2-1. 次を計算する。 $$ d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2 $$
ステップ2-2. $i = k+1 , k+2 , \cdots , n-1$ に関して次を計算する。 $$ l_{ik} = {{1} \over {d_{kk}}} \left\{ a_{ik} - \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} \right\} $$
ステップ3. $k = n$ に関して次を計算する。 $$ d_{nn} = a_{nn} - \sum_{\nu = 1}^{n-1} d_{\nu \nu} l_{n \nu}^2 $$
- 手計算では複雑に見えたが、アルゴリズムとして整理してみると、計算量が大幅に減少したことが分かる。いつでも$D$ の対角成分を先に決めた後、$L$ の値を列ごとに計算することを覚えておこう。
コード
今度は、LU 分解でやったことと同じ例を取り上げて、実際に LDU 分解を行ってみよう。
以下のコードは、R でそのアルゴリズムを実際に実装したものである。
B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
LDU<-function(A)
{
n=length(A[,1])
L<-diag(1,n)
D<-diag(0,n)
k=1
D[1,1] = A[1,1]
L[,1] = A[,1]/D[1,1]
for(k in 2:(n-1))
{
D[k,k] = A[k,k] - sum(diag(D)[1:(k-1)]*L[k,1:(k-1)]^2)
for(i in (k+1):n)
{
L[i,k] = ( A[i,k] - sum(L[i,1:(k-1)]*diag(D)[1:(k-1)]*L[k,1:(k-1)]) )/D[k,k]
}
}
D[n,n] = A[n,n] - sum(diag(D)[1:(n-1)]*L[n,1:(n-1)]^2)
return(list(L=L, D=D))
}
LDU(B)
LDU(B)$L%*%LDU(B)$D%*%t(LDU(B)$L)
실행 결과는 다음과 같다.
R 에는 따로 LDU 분해가 구현되어 있지 않아 검산만 해보았는데, 정확하게 분해가 되었음을 알 수 있다.
$B$ サイズがあまり大きくなく、計算もそこそこ管理可能なので、疑問があれば手計算をしてみてもいい。