대칭행렬의 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$ 의 사이즈가 그렇게 크지 않고 계산도 그럭저럭 할만한 편이므로 정 의심된다면 직접 손계산을 해봐도 좋다.