[파이썬과 선형대수] LU분해 그리고 직접 구현해보기

LU 분해(LU decomposition)

이전의 QR분해와 비슷하게 주어진 행렬을 두 행렬의 곱으로 분해하는 과정이다.
QR분해 때는 직교행렬인 Q와 행렬 R로 분해했었는데 LU분해의 경우 하삼각 행렬 L과 상삼각 행렬 U로 분해하는 것이 차이점이다.

$\begin{pmatrix}2&2&4 \\ 1&0&3\\2&1&2\end{pmatrix}=\begin{pmatrix} 1&0&0 \\ \frac{1}{2}&1&0\\1&1&1\end{pmatrix}\begin{pmatrix} 2&2&4\\0&-1&1 \\ 0&0&-3\end{pmatrix}$

보는 바와 같이 L은 lower triangular matrix로 대각의 위쪽 성분은 모두 0이다. 반대로 U는 upper triangular matrix로 대각 아래 성분이 모두 0인 것을 알 수 있다.

LU 분해 방법

LU 분해에 대해 알아보았다면 하는 방법을 알아보자.
LU분해를 하는 방법에는 기본 행렬을 이용하는 방법이 대표적이다.
어떤 방식을 사용하던 간에 LU분해는 가우스 소거법에서 비롯된다. 먼저 기본 행렬을 이용하는 방식에 대해 알아보자.

1. 기본 행렬을 이용하는 방법

가우스 소거법을 진행해보면 행렬은 사다리꼴 행렬(Row Echelon Matrix)만 남게 된다. 사다리꼴 행렬은 상삼각 행렬이므로 가우스 소거법을 진행한 후의 사다리꼴 행렬이 LU분해에서의 U라고 볼 수 있다.

가우스 소거법은 기본 행 연산에 의해 진행된다. 그리고 기본 행 연산을 표현하는 행렬이 기본 행렬이다. 하나의 기본 행렬에 하나의 기본 행 연산을 담기 때문에 사다리꼴 행렬이 되기까지 여러 개의 기본 행렬이 곱해질 것이다.
$$E_nE_{n-1}...E_2E_1A=U$$
$$A=(E_nE_{n-1}...E_2E_1)^{-1}U$$
$$L=(E_nE_{n-1}...E_2E_1)^{-1}\cdot\cdot\cdot(1)$$

예시

$$A=\begin{pmatrix}2&3&1 \\ 4&7&3\\6&18&5\end{pmatrix}$$
$$E_1=\begin{pmatrix} 1&0&0 \\ -2&1&0\\0&0&1\end{pmatrix},E_2=\begin{pmatrix} 1&0&0 \\ 0&1&0\\-3&0&1\end{pmatrix},A^{'}=\begin{pmatrix} 2&3&1\\0&1&1 \\ 0&9&2\end{pmatrix}$$
$$E_3=\begin{pmatrix} 1&0&0\\0&1&0 \\ 0&-9&1\end{pmatrix},A^{''}=U=\begin{pmatrix} 2&3&1\\0&1&1 \\ 0&0&-7\end{pmatrix}$$
위 식 (1)에 의해
$$L=\begin{pmatrix} 1&0&0\\2&1&0 \\ 3&9&1\end{pmatrix}$$

2. LU분해를 쉽게 하는 방법

기본행렬을 이용한 LU분해는 기본행렬을 하나하나 기록해놓은 다음에 한꺼번에 곱해주어야하기 때문에 빠르게 구하기 쉽지 않다.
정식적인 방법은 아니지만 꼼수 비슷하게 조금 더 쉽게 구하는 방법이 있다.
먼저 주어진 행렬을 똑같이 가우스 소거법을 해서 상삼각행렬을 만드는 과정을 진행한다.
그리고 규칙에 따라서 L을 작성하면 된다.
규칙 1. A의 대각 원소를 1로 바꾸기 위해 곱하는 수의 역수가 L의 대각 원소가 된다.
규칙 2. A의 대각 원소 아래에 위치한 원소를 0으로 만들기 위해 곱해주는 계수의 음수가 L의 같은 위치 원소가 된다.
규칙 2가 조금 헷갈릴 것이다. 위의 같은 행렬 A를 사용해 예를 들어보자.
$$A=\begin{pmatrix}2&3&1 \\ 4&7&3\\6&18&5\end{pmatrix}$$
먼저 A의 1행1열 원소 2를 1로 만들기 위해 1행에 $\frac{1}{2}$를 곱해야 하므로 L의 1행1열 원소는 그의 역수인 2다.

$$L_1=\begin{pmatrix}2&0&0 \\ 0&0&0\\0&0&0\end{pmatrix},A^{'}=\begin{pmatrix}1&\frac{3}{2}&\frac{1}{2} \\ 4&7&3\\6&18&5\end{pmatrix}$$
편의상 나머지 원소는 0으로 두자.
이후 1행에 -4를 곱해서 2행에 더한다. 이때, A의 2행1열 원소를 0으로 만들기 위함이므로 L의 2행1열 원소는 (-) 부호를 붙인 4가 된다.
추가로 1행에 -6를 곱해서 3행에 더한다. 똑같이 해보면 L의 3행1열 원소가 6이 된다.
$$L_2=\begin{pmatrix}2&0&0 \\ 4&0&0\\6&0&0\end{pmatrix}$$
이제 A는
$$A^{''}=\begin{pmatrix}1&\frac{3}{2}&\frac{1}{2} \\ 0&1&1\\0&9&2\end{pmatrix}$$
2행에 -9를 곱해서 3행에 더해야 한다. 그러므로 L의 3행2열 원소는 9다.
$$L_3=\begin{pmatrix}2&0&0 \\ 4&0&0\\6&9&0\end{pmatrix}$$
$$A^{''}=\begin{pmatrix}1&\frac{3}{2}&\frac{1}{2} \\ 0&1&1\\0&0&-7\end{pmatrix}$$
상삼각행렬이 완성되었다면, L의 대각은 1로 채워주어도 된다.
$$L=\begin{pmatrix}2&0&0 \\ 4&1&0\\6&9&1\end{pmatrix},U=\begin{pmatrix}1&\frac{3}{2}&\frac{1}{2} \\ 0&1&1\\0&0&-7\end{pmatrix}$$
혹은 U의 대각 원소를 모두 1로 만들어도 된다.
$$L=\begin{pmatrix}2&0&0 \\ 4&1&0\\6&9&-7\end{pmatrix},U=\begin{pmatrix}1&\frac{3}{2}&\frac{1}{2} \\ 0&1&1\\0&0&1\end{pmatrix}$$

처음 접해보는 방식이라 어려울 수 있지만, 기본행렬 여러 개를 구하거나 계산하지 않는다는 점에서 간편하다.

치환행렬의 필요성

치환행렬을 이전 글에서 알아보긴 했는데 왜 필요할까?

$$L=\begin{pmatrix}l_{11}&0&0 \\ l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\end{pmatrix},U=\begin{pmatrix}u_{11}&u_{12}&u_{13} \\ 0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}$$
LU분해는 위처럼 하삼각행렬과 상삼각행렬의 곱으로 나타난다. 이때, 분해하려는 행렬 A의 1행1열 원소가 0이라면 LU분해를 할 수 없다. 왜냐하면 $l_{11}\times u_{11}$ 값이 0이기 때문에 둘 중 하나의 값이 0이라는 뜻인데 이는 곧 L 또는 U가 특이행렬임을 의미하기 때문이다.

여기서 치환행렬의 필요성이 드러난다. 행이나 열을 바꿈으로써 대각 성분의 0을 없애면 LU분해를 할 수 있게 된다.
이렇게 치환행렬이 들어간 LU분해는 $PA=LU$로 나타낼 수 있고, 치환행렬 P는 직교행렬이기 때문에 $A=P^{-1}LU=P^TLU$로 표현할 수 있다.

LU분해의 고유성

가우스 소거법을 진행했을 때, 사다리꼴 행렬이 고유하지 않으므로, LU분해도 고유하지 않다. 다시 말해, 어떤 행렬을 하삼각행렬과 상삼각행렬로 나누는 경우는 무한하다.
다만, L의 대각 성분이 모두 1이라는 조건을 추가한다면, 최대계수 정방행렬 A의 LU분해는 고유하다고 말할 수 있다.

---

LU분해의 응용

선형방정식의 풀이

선형 방정식을 푸는 데에는 여러 방식이 있겠지만 LU분해를 이용해서도 해를 찾을 수 있다.
$$Ax=b$$
$$L(Ux)=b$$
이때, $Ux=y$로 치환하여 $Uy=b$를 먼저 계산해 $y$를 구한다.

$$Ux=y$$
다음으로 $Ux=y$를 계산해 $x$를 찾으면 된다.

역행렬 구하기

삼각 행렬의 역행렬은 일반적인 행렬에 비해 쉽게 구할 수 있다.
상삼각행렬의 역행렬은
$$A=\begin{pmatrix}a_{11}&a_{12}&a_{13} \\ 0&a_{22}&a_{23}\\0&0&a_{33}\end{pmatrix}, A^{-1}=\begin{pmatrix}\frac{1}{a_{11}}&-\frac{a_{12}}{a_{11}a_{22}}&\frac{a_{12}a_{23}-a_{13}a_{22}}{a_{11}a_{22}a_{33}} \\ 0&\frac{1}{a_{22}}&-\frac{a_{23}}{a_{22}a_{33}}\\0&0&\frac{1}{a_{33}} \end{pmatrix}$$
공식이 존재한다.
비록 $2 \times 2$행렬만큼 간단하진 않지만 그래도 공식이 존재하는 만큼 일반적인 행렬보다는 계산이 간단한 편이라고 볼 수 있다.

그래서 LU분해를 이용하면 역행렬도 조금 더 단순히 구할 수 있다.
$$A=P^TLU$$
$$A^{-1}=(P^TLU)^{-1}$$
$$A^{-1}=U^{-1}L^{-1}P$$

행렬식 값 구하기

$4 \times 4$이상의 행렬부터 행렬식을 구하기 쉽지 않다. 이때 LU분해를 사용한다면 조금 더 편하게 행렬식을 구할 수 있다.

LU분해는 주어진 행렬을 하삼각행렬와 상삼각행렬의 곱으로 분해하는 것이라고 했다. 삼각행렬의 행렬식은 대각 성분의 곱과 같다는 성질을 이용하여 원래 행렬의 행렬식을 쉽게 구할 수 있다.
$$det(A)=det(L) \times det(U)$$
만약 L의 대각 성분이 모두 1이라면,
$$det(A)=det(U)$$
로 표현할 수도 있다.

---

파이썬으로 LU분해 구현하기

LU분해를 하는 함수는 scipy.linalg 모듈에 이미 구현되어 있다.
해당 함수에 대한 자세한 설명은 아래 링크에 자세하게 작성해놓았다.
https://aaaaaaaaaaayowooji.tistory.com/49

[파이썬 라이브러리]numpy.scipy 함수 분석(1)(lu)

 

구현되어 있는 함수를 사용하는 것도 좋지만, 직접 구현도 해보도록 하자.

알고리즘

LU분해를 구현하는 방법에는 여러 가지가 있지만 조금 간단히 구현하기 위해서 Doolittle 알고리즘으로 구현해보도록 하자.

Doolittle 알고리즘

$L$ 행렬의 대각 성분을 1로 고정하고 계산하는 알고리즘이다.

$$A=\begin{pmatrix}a_{11}&a_{12}&a_{13} \\ a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{pmatrix},L=\begin{pmatrix}1&0&0 \\ l_{21}&1&0\\l_{31}&l_{32}&1\end{pmatrix},U=\begin{pmatrix}u_{11}&u_{12}&u_{13} \\ 0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}$$
직접 계산을 해보면,
$$u_{11}=a_{11},\;u_{12}=a_{12},\;u_{13}=a_{13}$$
$$u_{22}=a_{22}-l_{21}u_{13},\;u_{23}=a_{23}-l_{21}u_{12},\; u_{33}=a_{33}-l_{31}u_{13}-l_{32}u_{23}$$
위 식을 통해서 공식을 만들어보면
$$U[i][j]=A[i][j]-\sum\limits_{k=0}^{i-1}L[i][k] \times U[k][j]\;\;\;(j>i)$$

$U$행렬의 요소를 모두 구했다면 $L$의 요소를 계산해야 한다.
$$a_{21}=l_{21}u_{11},\;a_{31}=l_{31}u_{11},\;a_{32}=l_{31}u_{12}+l_{32}u_{22}$$
$$l_{21}=\frac{a_{21}}{u_{11}},\;l_{31}=\frac{a_{31}}{u_{11}},\;l_{32}=\frac{1}{u_{22}}(a_{32}-l_{31}u_{12}-l_{32}u_{22})$$
각 요소는 위와 같이 계산되므로 아래와 같은 공식을 만들어 볼 수 있다.
$$L[i][j]=\frac{1}{U[j][j]}(A[i][j]-\sum\limits_{k=0}^{j-1}L[i][k] \cdot U[k][j])$$

$U$의 대각 요소가 0이면 LU분해가 불가능하므로 0일 때 예외처리를 해야 한다.

코드

import numpy as np

def doolittle_lu(A):
    n = len(A)
    L = np.eye(n)  # 대각선이 1인 단위 행렬
    U = np.zeros_like(A, dtype=np.float64)

    for i in range(n):
        # U 계산
        for j in range(i, n):
            U[i][j] = A[i][j] - sum(L[i][k] * U[k][j] for k in range(i))

        # L 계산
        for j in range(i + 1, n):
            if U[i][i] == 0:
                raise ValueError("LU 분해가 불가능합니다 (0으로 나눌 수 없음).")
            L[j][i] = (A[j][i] - sum(L[j][k] * U[k][i] for k in range(i))) / U[i][i]

    return L, U

# 테스트
A = np.array([[2,3,1], 
              [4,7,3],
              [6,18,5]])
L, U = doolittle_lu(A)

print("L 행렬:")
print(L)
print("U 행렬:")
print(U)

실행결과

L 행렬:
[[1. 0. 0.]
 [2. 1. 0.]
 [3. 9. 1.]]
U 행렬:
[[ 2.  3.  1.]
 [ 0.  1.  1.]
 [ 0.  0. -7.]]

위에서 했던 같은 행렬 A로 LU분해를 해보았다.
결과가 위에서 한 것과 같음을 알 수 있다.

scipy.linalg.lu와의 차이점

Doolittle 알고리즘을 사용해 직접 구현한 LU분해는 대각선에 0이 있으면 에러를 반환해버렸다. 하지만 위에서 알아보았듯이 치환행렬을 이용한다면 LU분해를 할 수 있는 경우가 존재한다. 이런 경우를 해결하기 위해서 scipy.linalg.lu는 치환행렬 P까지 고려하는 알고리즘을 사용한다. 그래서 scipy의 lu 함수는 P,L,U 세 개를 반환한다.

출처

서적

개발자를 위한 실전 선형대수학/한빛미디어/마이크,코헨 지음

Website

LU 분해 위키백과
LU분해에 대한 타 블로그 글(1)
LU분해에 대한 타 블로그 글(2)