対称行列の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}$ が可逆対称行列としよう。
**Step 1. $k = 1$
$d_{11} = a_{11}$ を代入して $\displaystyle l_{i1} = {{1} \over {d_{11}}} a_{i1}$ を計算する。
**Step 2. $k = 2, 3, \cdots , n-1$
Step 2-1. 次を計算する。 $$ d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2 $$
**Step 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\} $$
**Step 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$ のサイズがそれほど大きくなく、計算もそれなりにできるので、疑わしい場合は直接手計算しても良い。