[알고리즘 & 코딩테스트/코딩테스트] - 파이썬으로 행렬 구현하기 (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 이라고 합니다.
구현
지난 포스트에서 만든 행렬을 가지고 구현합니다. 구하는 과정은 위키를 참조하시길 바랍니다.
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 를 참고하시길 바랍니다.
'Theory > Algorithms' 카테고리의 다른 글
몬테카를로 시뮬레이션 (2) - 온라인 판매자 예제 (0) | 2022.02.26 |
---|---|
몬테카를로 시뮬레이션 (1) - 파이 계산하기 (0) | 2022.02.25 |
Trie (0) | 2021.07.15 |
Sparse Matrix (희소행렬) 구현하기 (0) | 2021.07.05 |
위상 정렬 (Topological Sort) (0) | 2021.03.05 |