Python을 이용한 Sholesky 분해 구현

Sholesky 분해(Cholesky Decomposition)는 대칭 행렬을 하삼각 행렬의 곱으로 분해하는 방법으로, 주로 선형 대수 문제의 해결 및 확률론, 통계학에서 활용된다. 여기에서는 Python을 사용하여 Sholesky 분해를 구현하는 방법을 단계별로 설명하겠다.

필요 라이브러리 임포트

import numpy as np

Sholesky 분해 개요

Sholesky 분해는 $\mathbf{A}$ 가 대칭이고 양의 준정부호(Positive Definite)인 $\mathbf{A} = \mathbf{L}\mathbf{L}^T$ 를 만족하는 하삼각 행렬 $\mathbf{L}$ 을 찾는 분해 기법이다. 여기서 $\mathbf{L}$ 는 하삼각 행렬이고, $\mathbf{L}^T$ 는 $\mathbf{L}$ 의 전치 행렬이다.

알고리즘 설명

Sholesky 분해를 구하는 일반적인 알고리즘은 다음과 같다:

$\mathbf{L}$ 를 영행렬로 초기화한다. ( $\mathbf{L}_{ij} = 0$ )
반복을 통해 각 행과 열의 요소를 아래의 식을 사용하여 계산한다:

$L_{ii} = \sqrt{A_{ii} - \sum_{k=1}^{i-1} L_{ik}^2}$

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

Python으로 구현

아래는 Python을 이용하여 Sholesky 분해를 구현하는 코드이다:

def cholesky_decomposition(A):
    """
    Perform a Cholesky decomposition of matrix A.
    A must be a symmetric and positive-definite matrix.

    :param A: A symmetric, positive-definite matrix
    :return: Lower triangular matrix L such that A = L L^T
    """
    n = A.shape[0]
    L = np.zeros_like(A)

    for i in range(n):
        for j in range(i+1):
            if i == j:  # Diagonal elements
                sum_sq = sum(L[i][k]**2 for k in range(j))
                L[i][j] = np.sqrt(A[i][i] - sum_sq)
            else:
                sum_prod = sum(L[i][k] * L[j][k] for k in range(j))
                L[i][j] = (A[i][j] - sum_prod) / L[j][j]

    return L

사용 예시

위의 함수를 테스트하기 위해 예제 행렬을 사용해 보자. 다음은 $\mathbf{A}$ 가 대칭이며 양의 준정부호인 예제이다:

A = np.array([[4, 12, -16],
              [12, 37, -43],
              [-16, -43, 98]])

L = cholesky_decomposition(A)
print("Lower triangular matrix L:")
print(L)

print("Verify A = L * L.T:")
print(np.dot(L, L.T))

이 코드를 실행하면 $\mathbf{L}$ 과 $\mathbf{L}^T$ 의 곱이 $\mathbf{A}$ 를 재구성하는 것을 볼 수 있다.

np.dot() 함수를 사용해 $\mathbf{L}$ 과 $\mathbf{L}^T$ 를 곱하여 원본 행렬 $\mathbf{A}$ 와 동일한지 검증한다.

코드 설명

cholesky_decomposition 함수는 입력으로 행렬 $\mathbf{A}$ 를 받으며, $\mathbf{L}$ 을 반환한다.
입력 행렬 $\mathbf{A}$ 는 대칭이고 양의 준정부호여야 한다.
n 은 행렬 $\mathbf{A}$ 의 크기를 나타낸다.
$\mathbf{L}$ 은 $\mathbf{A}$ 와 동일한 크기의 영행렬로 초기화된다.
이중 반복문을 통해 $\mathbf{L}$ 의 각 요소를 계싼한다.

완성된 코드 실행 결과

위에서 설명한 코드를 실행하면 다음과 같은 결과가 나올 것이다:

Lower triangular matrix L:
[[ 2.          0.          0.        ]
 [ 6.          1.          0.        ]
 [-8.         -5.          3.        ]]

Verify A = L * L.T:
[[  4.  12. -16.]
 [ 12.  37. -43.]
 [-16. -43.  98.]]

이 결과를 통해 우리가 구한 하삼각 행렬 $\mathbf{L}$ 이 원본 행렬 $\mathbf{A}$ 를 재구성함을 확인할 수 있다. 이는 코드가 올바르게 작동함을 의미한다.

Python 내장 함수와 비교

원하는 경우, Python의 라이브러리인 numpy에서 제공하는 내장 함수 numpy.linalg.cholesky를 사용하여 결과를 검증할 수도 있다:

L_numpy = np.linalg.cholesky(A)
print("Lower triangular matrix using numpy.linalg.cholesky:")
print(L_numpy)

이 결과는 우리가 수작업으로 만든 Cholesky 분해와 일치해야 한다. numpy를 활용하면 더욱 간편하게 작업을 수행할 수 있지만, 직접 구현하면서 내부 알고리즘을 깊이 이해할 수 있다.

에러 처리

실제 응용 환경에서는 입력 행렬이 항상 조건에 부합하지 않을 수 있다. 따라서, 에러 처리를 통해 이러한 상황을 미리 예측하고 적절한 조치를 취하는 것이 중요하다.

def cholesky_decomposition_safe(A):
    """
    Perform a Cholesky decomposition of matrix A with error handling.
    A must be a symmetric and positive-definite matrix.

    :param A: A symmetric, positive-definite matrix
    :return: Lower triangular matrix L such that A = L L^T
    """
    if A.shape[0] != A.shape[1]:
        raise ValueError("Matrix A must be square")

    if not np.allclose(A, A.T):
        raise ValueError("Matrix A must be symmetric")

    try:
        return cholesky_decomposition(A)
    except ValueError as e:
        raise ValueError(f"Matrix A is not positive definite: {e}")

try:
    L_safe = cholesky_decomposition_safe(A)
    print("Lower triangular matrix L with safe checks:")
    print(L_safe)
except ValueError as e:
    print(e)

Cholesky 분해는 선형 대수와 통계학에서 매우 유용한 기법으로, Python을 이용하면 이를 쉽게 구현하고 사용할 수 있다. 이를 직접 구현하면서 이론적인 이해와 함께 코드의 실체를 파악하는 좋은 기회가 되었다. 에러 처리와 같은 실제 응용 상황에서의 문제 해결 방법도 포함하여 보다 견고한 소프트웨어를 개발할 수 있다.