logo

대칭행렬의 LDU 분해 📂행렬대수

대칭행렬의 LDU 분해

개요

LTL^{T} 는 상삼각행렬이므로 A=LUA = LU 에서의 UUU:=DLTU:= DL^{T} 으로 바꾼다고 보면 된다. 일반적인 LU 분해보다 조건이 까다로워진만큼 계산량은 많이 줄어든다.

정리

가역대칭행렬 AA 에 대해 A=LDLTA = LDL^{T} 를 만족하는 하삼각행렬 LL대각행렬 DD 가 존재한다.

증명

AA가역행렬이므로 A=LUA = L U 를 만족하는 하삼각행렬 LL상삼각행렬 UU 가 존재한다. 한편 A=ATA = A^{T} 이므로

LU=A=AT=UTLT L U = A = A^{T} = U^{T} L^{T}

LU 분해에서 LL 은 대각선 성분이 모두 11 인 하삼각행렬로 구해지므로 detL=1\det L = 1, 즉 역행렬 L1L^{-1} 이 존재한다. L1L^{-1} 는 하삼각행렬이고 LU=UTLT L U = U^{T} L^{T} UU 에 대해 정리하면 U=L1UTLT U = L^{-1} U^{T} L^{T} 양변의 오른쪽에 (LT)1(L^{T})^{-1} 을 곱하면 U(LT)1=L1UT U (L^{T})^{-1} = L^{-1} U^{T} 을 얻는다. U(LT)1=L1UTU (L^{T})^{-1} = L^{-1} U^{T} 의 좌변은 상삼각행렬이고 우변은 하삼각행렬이므로 D:=U(LT)1=L1UT D: = U (L^{T})^{-1} = L^{-1} U^{T} 대각행렬임을 알 수 있다. 이를 UU 에 대해서 정리하면 U=DLTU = DL^{T} 이고 다음을 얻는다. A=LU=LDLT A = LU = LDL^{T}

계산

이제 구체적으로 LLDD 를 계산하기 위해 아래와 같이 정의해보자.

L:=[1000l21100l31l3210ln1ln2ln31],D:=[d110000d220000d330000dnn] 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=(aij)=LDLTA = (a_{ij}) = LDL^{T} 를 직접 계산해보자.

lij>0l_{ij} > 0 이고 lii=1l_{ii} = 1 이므로

aij=ν=1nliνμ=1ndνμlμjT=ν=1nliνμ=1nδνμdνμljμ=ν=1nliνdννljν=ν=1min(i,j)liνdννljν \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*}

우선 jij \le i 인 경우로 가정해보자.

aij=ν=1jliνdννljν=ν=1j1liνdννljν+lijdjjljj=ν=1j1liνdννljν+lijdjj 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=ki = j = k 인 경우

akk=ν=1k1lkνdννlkν+dkkdkk a_{kk} = \sum_{\nu = 1}^{k-1} l_{k \nu} d_{\nu \nu} l_{k \nu} + d_{kk} d_{kk}

에 대해 정리하면

dkk=akkν=1k1dννlkν2 d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2

반대로 j=kj = k 이고 i=k+1,,ni = k+1, \cdots , n 이라고 하면

aik=ν=1k1liνdννlkν+likdkk a_{ik} = \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} + l_{ik} d_{kk}

likl_{ik} 에 대해 정리하면

lik=1dkk{aikν=1k1liνdννlkν} l_{ik} = {{1} \over {d_{kk}}} \left\{ a_{ik} - \sum_{\nu = 1}^{k-1} l_{i \nu} d_{\nu \nu} l_{k \nu} \right\}

알고리즘

(aij)Rn×n(a_{ij}) \in \mathbb{R}^{n \times n} 가 가역대칭행렬이라고 하자.

**Step 1. k=1k = 1

d11=a11d_{11} = a_{11} 을 대입하고 li1=1d11ai1\displaystyle l_{i1} = {{1} \over {d_{11}}} a_{i1} 을 계산한다.


**Step 2. k=2,3,,n1k = 2, 3, \cdots , n-1

  • Step 2-1. 다음을 계산한다. dkk=akkν=1k1dννlkν2 d_{kk} = a_{kk} - \sum_{\nu = 1}^{k-1} d_{\nu \nu} l_{k \nu}^2

  • **Step 2-2. i=k+1,k+2,,n1i = k+1 , k+2 , \cdots , n-1 에 대해 다음을 계산한다. lik=1dkk{aikν=1k1liνdννlkν} 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=nk = n 에 대해 다음을 계산한다. dnn=annν=1n1dννlnν2 d_{nn} = a_{nn} - \sum_{\nu = 1}^{n-1} d_{\nu \nu} l_{n \nu}^2


  • 손으로 계산할 땐 복잡해보였지만 알고리즘으로 정리해서 보면 계산량이 확 줄어들었음을 알 수 있다.언제나 DD 의 대각성분을 먼저 구한 후 LL 의 값을 한 열 단위로 구해나간다는 것을 기억하자.

코드

이제 LU 분해에서 했던 것과 같은 예제를 가져와서 LDU 분해를 실제로 해보자.

5A1AC6BB2.png

아래의 코드는 위의 알고리즘을 실제로 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)

실행 결과는 다음과 같다. 20171128\_143203.png

R 에는 따로 LDU 분해가 구현되어 있지 않아 검산만 해보았는데, 정확하게 분해가 되었음을 알 수 있다.

BB 의 사이즈가 그렇게 크지 않고 계산도 그럭저럭 할만한 편이므로 정 의심된다면 직접 손계산을 해봐도 좋다.