Cholesky Decomposition (촐레스키 분해) - 파이썬 구현

[알고리즘 & 코딩테스트/코딩테스트] - 파이썬으로 행렬 구현하기 (4) - 마무리

---

지난 포스트에서 numpy 없이 구현한 행렬을 이용해 cholesky decomposition (촐레스키 분해) 를 해보도록 하겠습니다. 먼저 정방행렬 $A$가 있을 때, $x\neq 0$인 모든 $x$에 대해 $x^tAx>0$ 이면 이 행렬 $A$를 positive definite 하다고 합니다. ($x^tAx\geq 0$인 경우 semi-positive definite 하다고 합니다) 이때 행렬 $A$가 대칭이고 positive definite 하다면 $A=LL^t$를 만족하는, 대각 위의 성분이 모두 0인 하삼각행렬 (lower triangular matrix) $L$이 존재하게 되는데, $A$로부터 $L$을 구하는 것을 cholesky decomposition 이라고 합니다. 

 

구현

지난 포스트에서 만든 행렬을 가지고 구현합니다. 구하는 과정은 위키를 참조하시길 바랍니다.

https://en.wikipedia.org/wiki/Cholesky_decomposition

https://en.wikipedia.org/wiki/Cholesky_decomposition

def Cholesky(A):
    import copy, math
    if not is_almost_symmetric(A):  # 대칭
        raise ArithmeticError('not symmetric')

    L = copy.deepcopy(A)
    for k in range(L.ncols):
        if L[k,k] <= 0:  # 대각 원소는 제곱이 될 것이므로 0보다 커야함
            raise ArithmeticError('not positive definite')
        p = L[k,k] = math.sqrt(L[k,k])
        for i in range(k+1, L.nrows):  
            L[i,k] /= p
        for j in range(k+1, L.nrows):
            p = float(L[j,k])
            for i in range(k+1, L.nrows):
                L[i,j] -= p*L[i,k]
    ## Upper triangle
    for i in range(L.nrows):
        for j in range(i+1, L.ncols):
            L[i,j] = 0
    
    return L

 

행렬의 원소가 0인지 확인하는 is_almost_zero 함수를 만들어 $A-LL^t$의 값이 모두 0인지를 확인합니다.

def is_almost_zero(A, ap=1e-6, rp=1e-4):
    for r in range(A.nrows):
        for c in range(A.ncols):
            delta = abs(A[r,c])
            if delta > ap and delta > max(abs(A[r,c], abs[r,c]))*rp:
                return False

    return True
    
    
A = Matrix([[4,2,1],[2,9,3],[1,3,16]])
L = Cholesky(A)
print(is_almost_zero(A - L*L.T))
>>> True

---

입력 행렬이 대칭이 아니거나 positive definite이 아니라면 cholesky decomposition이 실패한다는 것으로 대칭 행렬이 positive definite 인지 확인하기 위해서도 사용할 수 있습니다.

def is_positive_definite(A):
    if not is_almost_symmetric(A):
        return False

    try:
        Cholesky(A)
        return True
    except RuntimeError:
        return False

 

적용

Cholesky decomposition의 간단한 적용으로는 다음과 같은 확률 분포를 따르는 벡터 $x$를 만드는 데 사용할 수 있습니다. 대칭이고 positive definite인 행렬 $A$에 대해 $p(x) \sim exp(-\frac{1}{2} x'A^{-1}x)$ 를 따르는 확률 분포가 있다고 가정합니다. Cholesky decomposition에 의해 $A=LL^t$이고 이를 대입하고 $u=L^{-1}x$로 바꾸면, 

$p(x) \sim exp(-\frac{1}{2} (L^{-1}x)'(L^{-1}x))$

$p(x) \sim exp(-\frac{1}{2} u'u)$

이고 벡터 $u=(u_0,u_1,...)$라 정의하면 다음을 얻을 수가 있는데, $u_i$는 평균이 0이고 표준편차가 1인 가우스 확률 변수가 됩니다.

$p(x) \sim exp(-\frac{1}{2} u_0^2)exp(-\frac{1}{2} u_1^2)exp(-\frac{1}{2} u_2^2) ...$

요약하면 공분산 행렬 $A$가 주어졌을 때, 같은 공분산을 지닌 확률 벡터 $x$ 또는 난수를 만들 수 있다는 것이죠.

import random
def RandomList(A):
    L = Cholesky(A)
    while True:
        u = Matrix([[random.gauss(0,1)] for c in range(L.nrows)])
        yield (L *u).flatten()

• $u_i$는 평균 0, 표준편차가 1인 가우스 확률 변수이므로 "random.gauss(0,1)"로 생성합니다.
• $u=L^{-1}x$로 정의했으므로 $x=Lu$로 구합니다.
• RandomList는 제너레이터 함수로 생성자이며 반복적으로 호출할 수 있습니다. 자세한 사항은 [개발 잡학/Python] - 제너레이터와 yield 를 참고하시길 바랍니다.