LU분해 5분컷 이해: python

대부분의 행렬A는 행렬을 일단 L(하삼각행렬 + 대각행렬은 1)과 U(상삼각행렬)로 분해가 된다. 분해를 왜하냐고 묻는다면, 일단 분해하고나면, A의 특성을 구하는데, 여러모로 행렬연산이 편해진다. 예를 들어, LU을 분해하고나면 det(LU)=det(L)det(U)로 바로 계산할수도있고, 일차방정식, 역행렬도 구할 수 있다.

 

L(하삼각행렬): 다음의 꼴과 같이 하삼각행렬만 0이 아닌 값이 있고, 대각행렬은 모두 0이어야 한다. 유사하게 U(상삼각행렬)은 위의 삼각행렬을 제외하고 모두 0이어야한다. 다만 L의 대각행렬이 1이기 때문에, 상삼각행렬은 1일 필요는 없다. 대각행열까지를 포함한다. Python의 Numpy은 Tri+(l, u)을 통해 이를 구현하고있다. np.tril은 하삼각행렬을 만들어주는 함수이고, np.triu은 상삼각행렬을 만들어준다. 

import numpy as np

# 하삼각행렬
L = np.tril([[1,2,3],
            [4,5,6],
            [7,8,9]], k=-1)
I = np.identity(L.shape[0])
L = L + I

print(L)

[[1. 0. 0.]
 [4. 1. 0.]
 [7. 8. 1.]]

# 상삼각행렬
U = np.triu([[1,2,3],
            [4,5,6],
            [7,8,9]], k=0)
print(U)

[[1 2 3]
 [0 5 6]
 [0 0 9]]

 

LU 분해는 numpy로 직접구하거나 scipy에 지원하는 lu 분해를 이용한다. scipy은 p matrix을 주는데, 행렬을 순서를 변경하기위한 행렬이다. LU분해를 할 때, L의 대각행렬이 아무리 연산하려고해도 0이어서, 다음 하삼각행렬을 구할 수 없을때, 피봇팅(pivotting)해주게되는데, scipy.linalg.lu 은 이 피봇팅한것을 원복할 때 쓰는 P을 가지고있다. 따라서, $A=PLU$으로 연산할 수 있고, 내부적으로 연산한 L, U을 반환해주기도 한다 (permute_l=True: Perform the multiplication P*L). 보통은 이 P을 치환행렬(Pivoting 을 위한 행렬)이라고한다.

 

# scipy을 이용하는 방법

from scipy.linalg import lu as lu_decompose

mat = np.array([[1,2,3],
                [2,6,4],
                [8,9,1]])

L, U = lu_decompose(mat, permute_l=True)  # p matrix을 안받으려면 True로 설정

def lu_decompose(mat):
    rows, columns = mat.shape
    s = min(rows, columns)  # to determine s: (row by s) @ (s by colunns) 
    
    for k in range(s):
        x = 1 / mat[k,k]
        for i in range(k+1, rows):
            mat[i,k] = mat[i,k] * x
        for i in range(k+1, rows):
            for j in range(k+1, columns):
                mat[i,j] = mat[i,j] - mat[i,k] * mat[k,j]
            
    L = np.tril(mat, k=-1) + np.identity(rows)
    U = np.triu(mat, k=0)
    return L, U


mat = np.array([[1,2,3],
                [2,6,4],
                [8,9,1]])
L, U = lu_decompose(mat)

 

 

LU분해 사용법

---

1. det(LU)=det(L)det(U)이다.

det_A = np.linalg.det(mat)
detLU = np.linalg.det(L) * np.linalg.det(U)

print(det_A, detLU)
# -59.999999999999986 -59.999999999999986

 

2. 일차방정식을 LU로 풀 수 있다.

$Ax=y$ 이었으므로, $LUx=y$라고 쓸 수있다. 따라서, $x$에 $u$을 곱하고, $L$을 다시 곱하면, y가 된다 라는 말이다. 반대로, x을 찾기위해서는 y가 되는 Lz을 구하고, z가 되는 Ux의 x을 구하면된다.

 

$Ax=b$의 식을 구하는 것을 보면, x은 해를 구하는 방법이다. 한편, $Ax=I$으로 $b$을 $I$으로 바꾸면, A을 역행렬을 구하는 작업이다.

 

우선A의 LU을 구했다고치자. 다음의 전개가 가능하다.

 

$LUx=I$

 

$Ly=I, where Ux=y$

 

따라서, $Ly=I$을 먼저 구한 다음에 순차적으로 x까지구하면된다. 문제는 이러한 방법이 왜 편리하냐를 의미하는 것인데, 이는 전방대치(forward substitute)와 후방대치(backward substitute)을 쓸 수 있어서 그렇다. 아래와 같이 전방치(forward substitution)방법으로 구하는 방법은 하삼각행렬만 숫자가 0이 아니고, 나머지는 0이기에, 직관적으로 x1은 5/3임을 알 수 있다. 이를 구했기에, 2행..4행까지 바로 다 구할 수 있다. 그렇기에 편리하다는 의미이다.

 

Reference: https://www.essaytaste.com/solved-write-python-code-perform-forward-substitution-opposed-back-substitution-algorithm-solve-q37161323/