숄레스키 분해

Prerequisites

이번 포스팅의 내용을 더 잘 이해하기 위해선 아래의 내용에 대해 알고 오시는 것이 좋습니다.

LU 분해

LU 분해를 수행하는 또 다른 방법

LU 분해 편에서는 LU 분해란 Gaussian elimination을 수행하는 과정에서 사용하는 기본 행 연산을 이용해 얻게 되는 행렬 분해 방법이라고 소개한 바 있다.

그런데, 꼭 Gaussian elimination을 이용하지 않고 아래와 같이 행렬 $A$ 를 하삼각행렬과 상삼각행렬의 곱으로 분해된다고 가정하더라도 LU 분해의 결과를 그대로 얻을 수 있을 것이다.

임의의 $3\times 3$ 사이즈의 행렬 $A$ 에 대해 다음과 같이 행렬 분해가 가능하다고 생각해보자.

\begin{matrix} (1) & [\begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix}] = [\begin{matrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{matrix}] [\begin{matrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{matrix}] \end{matrix}

$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \\ \end{bmatrix}$

\begin{matrix} (2) & = [\begin{matrix} u_{11} & u_{12} & u_{13} \\ u_{11} l_{21} & u_{12} l_{21} + u_{22} & u_{13} l_{21} + u_{23} \\ u_{11} l_{31} & u_{12} l_{31} + u_{22} l_{32} & u_{13} l_{31} + u_{23} l_{32} + u_{33} \end{matrix}] \end{matrix}

$=\begin{bmatrix} u_{11} & u_{12} & u_{13} \\ u_{11}l_{21} & u_{12}l_{21} + u_{22} & u_{13}l_{21}+u_{23} \\ u_{11}l_{31} & u_{12}l_{31} + u_{22}l_{32} & u_{13}l_{31}+u_{23}l_{32}+u_{33} \\ \end{bmatrix}$

이 결과만을 놓고 보면 $u_{11}$ , $u_{12}$ , $u_{13}$ 이 $a_{11}$ , $a_{12}$ , $a_{13}$ 과 같은 것을 알 수 있고, 다음 행의 내부 값은 구해놓은 $u_{11}$ , $u_{12}$ , $u_{13}$ 로부터 구할 수 있음을 알 수 있다.

이런 방식으로 순차적으로 $L$ 과 $U$ 의 원소들을 구할 수 있다.

Symmetric 행렬의 LU 분해?

Symmetric 행렬의 경우¹ LU 분해를 다음과 같은 방식으로도 생각할 수 있다.

만약 $A$ 가 symmetric matrix라면 혹시 이런식으로도 분해될 수 있지는 않을까?

\begin{matrix} (3) & A = L L^{T} = L^{T} L \end{matrix}

$A=LL^T=L^TL$

왜냐면 symmetric matrix는 $A=A^T$ 이므로 $(LL^T)^T = LL^T$ 라고 쓸 수도 있을 것 같고, $L^T$ 는 상삼각행렬이므로 $LU$ 분해와 유사한 형태의 결론을 얻을 수도 있을 것 같기 때문이다.

그런데 $A$ 가 symmetric 행렬이라고 해서 항상 $A=LL^T=L^TL$ 로 분해할 수 있는 것은 아니다. 우리는 임의의 $L$ 이 어떤 특성을 가져야만 하는지 생각해보자.

행렬 $L$ 과 임의의 벡터 $x$ 의 곱 $Lx$ 를 생각해보자. 이 $Lx$ 벡터의 L2-norm 값은 항상 0보다 크거나 같다.

그리고 L2-norm을 계산하는 방법은 벡터의 내적으로도 가능한데, 이 과정을 다시 쓰면,

\begin{matrix} (4) & | L x |^{2} = (L x)^{T} (L x) \end{matrix}

$|Lx|^2 = (Lx)^T(Lx)$

이다.

그리고 전치 연산자의 성질에 의해 다음과 같이 정리해서 쓸 수 있다.

\begin{matrix} (5) & x^{T} L^{T} L x \end{matrix}

$x^TL^TLx$

그리고 괄호로 중간에 있는 $L^TL$ 을 묶어주면,

\begin{matrix} (6) & x^{T} (L^{T} L) x \end{matrix}

$x^T(L^TL)x$

이며 $(L^TL)$ 이 어떤 행렬 $A$ 라고 한다면,

\begin{matrix} (7) & x^{T} A x \end{matrix}

$x^TAx$

와 같이 쓸 수 있으며 이 계산은 임의의 벡터 $Lx$ 의 L2-norm을 계산하는 방법으로부터 왔으므로

\begin{matrix} (8) & x^{T} A x \geq 0 \end{matrix}

$x^TAx\geq 0$

이다. 우리는 $x^TAx\geq 0$ 과 같은 성질을 만족하는 행렬을 semi positive definite 행렬이라고 부른다. 다시 말해, $A=LL^T=L^TL$ 을 만족하는 행렬은 semi-positive definite 행렬이어야 한다.

정리하면,

행렬 A가 정방행렬이고
대칭행렬이면서
semi-positive definite

일 때 $A=LL^T=L^TL$ 과 같이 분해 가능하며 이 분해 방법을 Cholesky factorization (숄레스키 분해)라고 부른다.

Cholesky factorization 계산

Cholesky factorization은 앞선 LU 분해의 또 다른 계산 방법에서와 같은 맥락으로 계산할 수 있다.

행렬 $A$ 가 대칭행렬이면서 semi-positive definite이라고 가정하자.

그러면 행렬 $A$ 를 다음과 같이 분해할 수 있다고 가정할 수 있다.

\begin{matrix} (9) & A = L L^{T} = [\begin{matrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33} \end{matrix}] [\begin{matrix} L_{11} & L_{21} & L_{31} \\ 0 & L_{22} & L_{32} \\ 0 & 0 & L_{33} \end{matrix}] \end{matrix}

$A = LL^T= \begin{bmatrix} L_{11} & 0 & 0 \\ L_{21} & L_{22} & 0 \\ L_{31} & L_{32} & L_{33} \\ \end{bmatrix} \begin{bmatrix} L_{11} & L_{21} & L_{31} \\ 0 & L_{22} & L_{32} \\ 0 & 0 & L_{33} \\ \end{bmatrix}$

행렬 곱을 계산해보면 다음과 같은 결과를 얻을 수 있다.

\begin{matrix} (10) & \Rightarrow [\begin{matrix} L_{11}^{2} & (대칭) \\ L_{21} L_{11} & L_{21}^{2} + L_{22}^{2} \\ L_{31} L_{11} & L_{31} L_{21} + L_{32} L_{22} & L_{31}^{2} + L_{32}^{2} + L_{33}^{2} \end{matrix}] \end{matrix}

$\Rightarrow \begin{bmatrix} L_{11}^2 & & \text{(대칭)}\\\\ L_{21}L_{11} & L_{21}^2+L_{22}^2 & \\\\ L_{31}L_{11} & L_{31}L_{21}+L_{32}L_{22} & L_{31}^2+L_{32}^2+L_{33}^2 \end{bmatrix}$

위의 계산 결과와 행렬 $A$ 의 원소를 1:1로 비교하면 다음과 같이 정리할 수 있다는 것을 알 수 있다.

\begin{matrix} (11) & L = [\begin{matrix} \sqrt{a_{11}} & 0 & 0 \\ a_{21} / L_{11} & \sqrt{a_{22} - L_{21}^{2}} & 0 \\ a_{31} / L_{11} & (a_{32} - L_{31} L_{21}) / L_{22} & \sqrt{a_{33} - L_{31}^{2} - L_{32}^{2}} \end{matrix}] \end{matrix}

$L=\begin{bmatrix} \sqrt{a_{11}} & 0 & 0 \\\\ a_{21}/L_{11} & \sqrt{a_{22}-L_{21}^2} & 0 \\\\ a_{31}/L_{11} & (a_{32}-L_{31}L_{21})/L_{22} & \sqrt{a_{33}-L_{31}^2-L_{32}^2} \end{bmatrix}$

이것의 패턴을 일반화 하면 다음과 같다는 것 또한 생각해볼 수 있다.

\begin{matrix} (12) & L_{j j} = \sqrt{a_{j j} - \sum_{k = 1}^{j - 1} L_{j k}^{2}} \end{matrix}

$L_{jj}=\sqrt{a_{jj}-\sum_{k=1}^{j-1}L_{jk}^2}$

\begin{matrix} (13) & L_{i j} = \frac{1}{L_{j j}} (a_{i j} - \sum_{k = 1}^{j - 1} L_{i k} L_{j k}) for i > j \end{matrix}

$L_{ij}=\frac{1}{L_{jj}}\left(a_{ij}-\sum_{k=1}^{j-1}L_{ik}L_{jk}\right)\quad \text{for }i>j$

MATLAB 구현

Cholesky factorization의 계산 방법을 직접 구현한 알고리즘 중 하나로 Cholesky–Banachiewicz algorithm 이 있다.

이 알고리즘은 Cholesky factorization의 하삼각행렬을 매 행마다 계산해나가는 방식으로 구현되어 있다.

아래는 MATLAB으로 Cholesky-Banachiewicz 알고리즘을 구현한 것이다.

A = [4,12,-16;12,37,-43;-16,-43,98];
L = zeros(size(A));

for i = 1:size(A,1)
    for j = 1:i
        my_sum = 0;
        
        for k = 1:j-1
            my_sum = my_sum + L(i,k)*L(j,k);
        end
        
        if i == j
            L(i,j)=sqrt(A(i,i)-my_sum);
        else
            L(i,j)=(1/L(j,j)*(A(i,j)-my_sum));
        end
    end
end

L
transpose(chol(A))

위 결과에서 L과 transpose(chol(A))를 비교한 결과가 같다는 것을 확인할 수 있을 것이다.

참고문헌

Cholesky decomposition, 위키피디아

복소 행렬이라면 에르미트 행렬로 확장할 수 있다. 다만, 에르미트 행렬의 경우에도 대칭 행렬의 경우와 같은 방법으로 Cholesky factorization의 아이디어를 생각할 수 있으므로 이번 포스팅에서는 실수 성분만을 갖는 symmetric 행렬에 한정해 생각해보도록 하자. ↩

이전크래머 공식의 기하학적 의미

다음양의 정부호 행렬