LU 분해(LU Decomposition) - (1) 분해 과정

LU 분해

---

  LU 분해는 연립방정식의 해를 푸는 방법 중 하나로, 일반적으로  $O(n^3)$의 시간복잡도를 가지는 가우스-조르당 소거법에 비해 적은 소요시간을 가지는 것으로 알려져 있다(실제로는 LU 분해를 사용해도 $O(n^3)$ 정도의 시간이 소요될 수 있다). LU 분해는 주어진 연립방정식 $Ax = b$에서의 $A$를  $L$, $U$로 표현되는 하삼각행렬(Lower Triangular Matrix)과 상삼각행렬(Upper Triangular Matrix)의 행렬곱으로 표현하는 과정과, 분해된 $L$,$U$ 행렬을 대입 연산(Back Substitution/Forward Substitution)을 통해 해를 구하는 두 과정으로 나뉘는데, 이번 포스트에서는 분해 과정에 대하여 설명하고, 해를 구하는 과정은 다음 포스트 참조

 

  LU 분해를 구현한 파이썬 알고리즘은 깃허브 링크에서 확인할 수 있다.

 

 

 

하삼각행렬과 상삼각행렬

---

  하삼각행렬과 상삼각행렬은, 각각 대각원소 아래/위 원소가 모두 0인, 다음과 같은 행렬을 말한다. 

 

하삼각행렬 : \begin{pmatrix} a_{1} & a_{2} & a_{3} \\ 0 & a_{4} & a_{5} \\ 0 & 0 & a_{6} \end{pmatrix}  상삼각행렬 :  \begin{pmatrix} a_{6} & 0 & 0 \\ a_{4} & a_{5} & 0 \\ a_{1} & a_{2} & a_{3} \end{pmatrix}

 

 

  이때, 대각원소를 제외한 모든 원소가 0인 대각행렬 또한 위 정의에 부합하므로, 대각행렬은 하삼각행렬이자 상삼각행렬이다.

 

 

  LU분해에서 $U$에 해당하는 상삼각 행렬은 $A$의 기본행 연산으로 생성된다. 이때, 다른 연산(다른 행과의 덧셈 등)을 외한 행교환 연산은 사용할 수 없다.  예외적으로, 첫 원소(1행x1열)가 0일 경우 행교환 연산을 1회(1열의 원소가 0이 아닌 다른 행과) 수행하게 되는데, 이때의 LU분해를 PLU 분해라고 한다.

 

 

  $A$를 기본행 연산을 $k$회 수행하여 상삼각행렬 $U$로 만드는 연산은 행렬 $A$의 앞에 기본행렬인 $E$들을 곱하는 것과 동일하므로 $ (E_{k}\cdots E_{1})\cdot A = U $의 형태로 나타낼 수 있다. 이때 상삼각행렬을 만드는 기본행 연산은 행렬 A의 대각선 아래 원소에 대해서만 연산을 수행하므로(해당 원소들을 0으로 만들기 위해), 기본행렬 E는 항상 하삼각행렬이다(물론 원소를 하나만 가질 수도 있음).

  행렬곱의 성질에 의해, 하삼각행렬끼리의 행렬곱은 항상 하삼각행렬이다(증명 생략). 따라서 기본행 연산을 통해 생성된 기본행렬들의 곱인 $(E_{k} \cdots E_{1})$을 전개하면 하삼각행렬이 되며, 이를 $E_{all}$ 이라고 정의했을 때 이 행렬의 역행렬인 $(E_{all})^{-1}$ 또한 하삼각행렬이다(증명 생략).

 

  즉, 상삼각행렬 $U$를 만들기 위해 $A$에 곱해진 기본행렬들을 모두 곱하고($E_{all}$), 이를 우변으로 넘겨 $U$와의 행렬곱 형태로 표현하면,

$E_{all}$의 역행렬인 $(E_{all})^{-1}$과 $U$의 곱이 각각 하삼각행렬과 상삼각행렬의 곱으로 표현되는 것이다. 

 

  수식으로 나타내면 다음과 같다.

 $(E_{k}\cdots E_{1})\cdot A = U $ #A를 기본행연산하여 U로 변환
$\Leftrightarrow E_{all}\cdot A = U $ #기본행연산을 위해 생성된 기본행렬을 모두 곱함
$\Leftrightarrow A = (E_{all})^{-1} \cdot U $ #기본행렬의 행렬곱의 역행렬을 취하여 U에 곱함
$\Leftrightarrow A = L \cdot U $

 

 

LU분해의 컴퓨터 알고리즘

---

 

앞서 설명했듯, 기본행 연산을 위한 기본행렬들의 곱에 대해 역행렬을 취한 다음 우변으로 넘기는 것이 기본적인 계산법이지만, 프로그래밍 언어로 이를 구현할 때는 기본행 연산에 맞는 기본행의 역행렬을 우변에 바로 곱하고, 이를 계산하는 것이 더 용이하다. 즉, 다음과 같이 연산이 진행되는 것이다.

 

$A = ((E_{1})^{-1}\cdots (E_{k})^{-1}) \cdot U $
$\Leftrightarrow A = (E_{all})^{-1} \cdot U $
$\Leftrightarrow A = L \cdot U $

 

$E$ 대신 $E^{-1}$를 바로 구한다는 점만 제외하면 큰 차이는 없으며, 이를 고려하여 LU 분해를 파이썬 코드로 다음과 같이 구현할 수 있다.

def lu_decomposition(matrix):
    
    if len(matrix) == len(matrix[0]):
        dim = len(matrix)
        L = np.zeros((dim,dim))
    else : 
        print("정방행렬이 아닙니다.")
 
 
    #첫 원소가 0일 경우 행교환이 필요
    while matrix[0][0] == 0:
        temp = matrix[0].copy()
        matrix[0] = matrix[1].copy()
        matrix[1] = temp.copy()
    
    
    U=matrix.copy()
    #입력 행렬의 1차 검증 : 값이 모두 0인 열이 있는지 -> 변수의 값을 구할 수 없음
    idx = 0
    for c in range(dim): #열
        for r in range(dim): #행
            idx += U[r][c]
 
        if idx == 0: #첫 열이 모두 0
            print(c+1,"열 의 값이 모두 0입니다.")
 
            
    for i in range(dim):
        ero_product = U[i][i].copy() # Elementary Row Operation(기본행연산)
        if ero_product == 0:
            pass
        else :
            L[i][i] = ero_product
            U[i] /= ero_product
            for j in range(i+1,dim):
                ero_sum = U[j][i].copy() 
                L[j][i] = ero_sum
                U[j] -= ero_sum*U[i].copy()
    
    #검증(L,U의 곱이 입력한 행렬과 같은지)
    if matrix.all() == np.dot(L,U).all():
        print("분해 완료")
    else : 
        print("분해 실패")
        return 
    
 
    return L,U

 

  lu_decomposition 함수는 단일 행렬을 입력값으로 받으며, 함수의 앞 부분에는 정방행렬인지 검증하는 코드와, 입력 행렬의 첫 원소가 0일 경우, 1행을 다른 행과 교체하여 PLU분해를 수행하는 알고리즘을 삽입하였다.

 

  LU분해의 경우 중첩 for 문을 통해 구현하였는데, 각 $k$열에 대해 $k+1$행부터 마지막 행까지(3x3 행렬의 경우 1열의 2,3행 / 2열의 3행) 연산을 반복하게 되면 정방행렬의 대각원소 아래의 원소만 소거되여 결과인 행렬 $U$가 상삼각행렬로 나타난다. 이때 행렬 $L$에 각각의 기본행 연산의 역연산(덧셈 <-> 뺄셈 / 곱셈 <-> 나눗셈)을 기록하면, 기본행 연산이 모두 끝났을 때 행렬 $L$은  $(E_{1})^{-1}\cdots (E_{k})^{-1}$

의 결과인 하삼각행렬로 표현된다.

 

 

행렬 A를$ A = \begin{pmatrix} 2 & 2 & 1 \\ 2 & 1 & 3 \\ 1 & 2 & 3 \end{pmatrix}$로 정의하고, LU분해한 결과는 다음과 같다.

 

$A = L \cdot U = \begin{pmatrix} 2 & 0 & 0 \\ 2 & -1 & 0 \\ 1 & 1 & 9/2 \end{pmatrix} \cdot \begin{pmatrix} 1 & 1 & 1/2 \\ 0 & 1 & -2 \\ 0 & 0 & 1 \end{pmatrix}$

 

 

 

 

참고자료

1. 프로그래머를 위한 선형대수(히라오카 카즈유키 저, 이창신 옮김) - 길벗(2017)

2. [Youtube] 쑤튜브 : 10분 선형대수 강의 모음(youtu.be/x-Ewz1ukXEA)