선형시스템 [특수 주제] LU분해

김록기·2023년 7월 27일

선형대수학 수학 파이썬

기초선형대수

목록 보기

5/13

다음과 같은 선형시스템을 생각해봅시다:

$x + 2y + 3z = 9$

$0x - y + 2z = 4$

$0x + 0y + z = 1$

3번째 행에서 $z = 1$ 을 구하고, 이를 2번째 행에 대입하여 $-y + 2(1) = 4$ 를 풀면 $y = -6$ 을 얻고, 이 결과들을 첫번째 행에 대입하는 방식으로 $x=18$ 을 얻을 수 있습니다. 이와 같은 방식을 역대입(Backward Substitution)이라고 하며, 역대입은 매우 빠르기 때문에, 삼각행렬로 표현된 선형시스템의 해는 빠르게 구할 수 있습니다.

이 선형시스템을 행렬과 벡터로 표현하면 다음과 같습니다:

$A = \begin{bmatrix} 1 & 2 & 3 \\ 0 & -1 & 2 \\ 0 & 0 & 1 \end{bmatrix}, \quad \bold{b} = \begin{bmatrix} 9 \\ 4 \\ 1 \end{bmatrix}$

이렇게, 행렬 $A$ 와 같은 형태의 행렬을 삼각행렬이라 부릅니다. 삼각행렬은 하삼각행렬과 상삼각행렬로 나눌 수 있습니다. 하삼각행렬은 주대각선 아래의 모든 원소가 0인 행렬이며, 상삼각행렬은 주대각선 위의 모든 원소가 0인 행렬입니다.

LU분해와 선형시스템의 풀이

LU 분해는 정사각행렬 A를 하삼각행렬(Lower Triangular Matrix) L과 상삼각행렬(Upper Triangular Matrix) U로 분해하는 기법을 말합니다. 즉, A = LU로 분해하는 것을 의미합니다.

LU 분해를 수행하면 선형시스템 $A\bold x = \bold b$ 가 매우 간단하게 해결됩니다.

Forward Substitution: 분해된 하삼각행렬 L과 원래의 벡터 $\boldsymbol{b}$ 를 이용하여 $L\boldsymbol{c} = \boldsymbol{b}$ 를 푸는 Forward Substitution 과정을 진행합니다. Forward Substitution은 L이 하삼각행렬이므로, 첫번째 행에서부터 순차적으로 해를 구해야합니다.

Backward Substitution: 다음으로, Forward Substitution 구한 벡터 $\boldsymbol{c}$ 를 이용하여 상삼각행렬 U와 해벡터 $\boldsymbol{x}$ 사이의 선형시스템 $U\boldsymbol{x} = \boldsymbol{c}$ 를 Backward Substitution으로 푸는 과정을 진행합니다. Backward Substitution은 U가 상삼각행렬이므로, 마지막 행에서부터 순차적으로 해를 구해야합니다.

LU분해의 방법

LU 분해의 방법은 주로 가우스 소거법(Gaussian Elimination)을 기반으로 합니다. 기본행연산들의 예시를 다시 보여주기 위해서, 여기서는 가우스 소거법을 기반으로 하는 방법을 소개합니다.

목적 : 정사각행렬 A를 하삼각행렬(Lower Triangular Matrix) L과 상삼각행렬(Upper Triangular Matrix) U로 분해합니다.
조건 : $A$ 의 모든 principal minor가 가역행렬이어야합니다. (물론, $A$ 가 대칭행렬일 필요는 없습니다.)
Step 1: 초기에는 L을 단위행렬로, U를 A로 설정합니다.
Step 2: i=1부터 $n-1$ 까지 다음을 반복합니다.
행렬 U의 i번째 행 $\bold u_i=\begin{bmatrix}U_{i1} & \ldots & U_{in}\end{bmatrix}$ 에 대해서, 다음을 반복합니다.
- $U$ 의 $j$ 번째 행 $\bold u_j$ 를 $\bold u_j -\dfrac{U_{ji}}{U_{ii}}\bold u_i$ 로 바꾸어줍니다. for $j >i$
  수학적으로 위의 과정은 다음과 같은 의미를 가집니다. 행렬 U에 대해서, 기본행연산 $E_{ij}$ 를 연속적으로 적용하여 $i$ 번째 행의 선행선분 $U_{ii}$ 아래의 모든 성분을 0으로 변환합니다.
  $U=E_{ij}U$ for $j >i$ .
- $L_{ij}$ 에 $\dfrac{U_{ji}}{U_{ii}}$ 을 저장합니다. for $j >i$
  수학적으로 위의 과정은 다음과 같은 의미를 가집니다. 행렬 L에 $U$ 에 적용한 기본행연산 $E_{ij}$ 의 역연산을 반대로 적용합니다. $L=L E_{ij}^{-1}$

모든 행렬이 LU분해가 가능한 것은 아닙니다. 그러나, 행렬 $A$ 가 LU분해 가능하지 않더라도, ( $A$ 가 가역이기만 하면) LU분해를 이용하여, 선형시스템을 풀수 있습니다. LU분해와 유사하지만, 순열행렬 P 추가로 포함한 형태로, 행렬 A를 다음과 같이 분해하여, 문제를 해결합니다.

$PA = LDU$ 또는 $PA = LU$

여기서 P는 순열행렬로서 행을 재배열하는 기능을 합니다. 그리고 행렬 $D$ 는 $L,U$ 의 대각성분을 1로 두기 위해서, 사용하는 대각행렬입니다.

LU분해의 수학적 설명

기본행렬의 예시들

기본행연산을 첨가행렬 $[A |\bold b]$ 에 적용하는 것은, 기본행렬을 $[A |\bold b]$ 의 왼쪽에 곱하는 것으로 볼 수 있습니다. 여기서, 기본행렬이란, 단위행렬에 기본행연산을 적용해서 얻어진 행렬입니다.

행의 교환

기본행렬 : $\begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix}, \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\0 & 1 & 0 \end{bmatrix}$ .
행에 0이 아닌 상수 $k$ 를 곱하기

기본행렬 : $\begin{bmatrix} k & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & k & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & k \end{bmatrix}$
한 행에 다른 행의 상수배( $k$ )를 더하기
기본행렬 : $\begin{bmatrix} 1 & 0 & 0 \\ k & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ k & 0 & 1 \end{bmatrix}, \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & k & 1 \end{bmatrix}$

가우스소거법 다시 보기

가우스소거법과정에서 사용되는 기본행연산에 대응하는 기본행렬은 교환행렬이거나, 하삼각행렬입니다.

다음 예시를 관찰해봅시다.

$[A|\bold b] = \begin{bmatrix} 2 & 1 & 1 & 5 \\ 4 & -6 & 0 & -2 \\ -2 & 7 & 2 & 9 \\ \end{bmatrix}\Rightarrow \begin{bmatrix} 1 & 0.5 & 0.5 & 2.5 \\ 0 & 1 & 0.25 & 0.75\\ 0 & 0 & 1 & 2 \end{bmatrix}$ .

위의 과정에 쓰인 기본행연산들을 코드와 수식을 혼용해서 표현하면 다음과 같습니다.

import numpy as np
row1 = np.array([2, 1, 1, 5])
row2 = np.array([4, -6, 0, -2])
row3 = np.array([-2, 7, 2, 9])
# Triangular Factors
row2 = row2 - 2*row1
row3 = row3+ row1
row3 = row3+ row2
# Normalize
row1 = 0.5*row1
row2 = -0.125*row2
row3 = 0.125*row3

단위행렬 $\begin{bmatrix}1 & 0 & 0\\0& 1&0\\0&0&1\end{bmatrix}$ 에 다음코드를 적용하면,

# Triangular Factors
row2 = row2 - 2*row1 # E_1
row3 = row3+ row1 #E_2
row3 = row3+ row2 #E_3

다음과 같습니다.

$E=E_3E_2E_1=\begin{bmatrix}1 & 0 & 0\\-2& 1&0\\-1&0&1\end{bmatrix}$

여기서, $E$ 가 하삼각행렬 인것은 우연이 아닙니다. 행렬 $E$ 가 하삼각행렬인 이유는, 다음과 같습니다.

$E_1$ 을 단위행렬에 row2 = row2 - 2*row1 # E_1을 적용한 결과,

$E_2$ 을 단위행렬에 row3 = row3+ row1 #E_2을 적용한 결과 그리고

$E_3$ 을 단위행렬에 row3 = row3+ row2 #E_3을 적용한 결과라 합시다. 그러면, 그들은 모두 하삼각 행렬입니다. 그리고 일반적으로 하삼각행렬의 곱은 하삼각행렬이므로, 그들의 곱 $E$ 역시 하삼각 행렬입니다.

그리고 단위행렬 $\begin{bmatrix}1 & 0 & 0\\0& 1&0\\0&0&1\end{bmatrix}$ 에 다음코드를 적용하면,

# Normalize
row1 = 0.5*row1
row2 = -0.125*row2
row3 = 0.125*row3

다음과 같습니다.
$\tilde{D}=\begin{bmatrix}0.5 & 0 & 0\\0& -0.125&0\\0&0&0.125\end{bmatrix}$

따라서, REF $\begin{bmatrix} 1 & 0.5 & 0.5 & 2.5 \\ 0 & 1 & 0.25 & 0.75\\ 0 & 0 & 1 & 2 \end{bmatrix}$ 은 다음과 같이 표현할 수 있습니다.

$\begin{bmatrix} 1 & 0.5 & 0.5 & 2.5 \\ 0 & 1 & 0.25 & 0.75\\ 0 & 0 & 1 & 2 \end{bmatrix}= \begin{bmatrix}0.5 & 0 & 0\\0& -0.125&0\\0&0&0.125\end{bmatrix}\begin{bmatrix}1 & 0 & 0\\-2& 1&0\\-1&0&1\end{bmatrix} \begin{bmatrix} 2 & 1 & 1 & 5 \\ 4 & -6 & 0 & -2 \\ -2 & 7 & 2 & 9 \\ \end{bmatrix}$

요약

LU분해는 가우스소거법을 계수행렬 $A$ 에만 적용한 것입니다. 위의 수식을 $\bold b$ 에 해당하는 부분을 제거하고 다시 적으면 다음과 같습니다. $U=\tilde{D}EA$
$\begin{bmatrix} 1 & 0.5 & 0.5 \\ 0 & 1 & 0.25 \\ 0 & 0 & 1 \end{bmatrix}= \begin{bmatrix}0.5 & 0 & 0\\0& -0.125&0\\0&0&0.125\end{bmatrix}\begin{bmatrix}1 & 0 & 0\\-2& 1&0\\-1&0&1\end{bmatrix} \begin{bmatrix} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \\ \end{bmatrix}$
이 때, (REF)의 정의에 의하면 좌변이 상삼각행렬인 것은 당연합니다. 이를 일반적으로 다음과 같이 적을 수 있습니다.
$U=\tilde{D} E A$ , 여기서 $E$ 는 기본행연산들의 곱으로 주어지는 행렬입니다.

약간의 디테일 : 만약, 가우스소거법에 쓰인 기본행연산에 행교환이 없으면, 가우스소거법의 특성상 대응하는 기본행연산은 언제나 하삼각행렬입니다. 언제, 행교환을 하지 않고 가우스소거법을 수행할 수 있는지, 그 필요충분조건은 행렬 $A$ 의 첫 $k$ 번째 행과 $k$ 번째 열안에 들어있는 원소들로 구성한 $k\times k$ 행렬이 가역행렬인 것입니다.

그리고, (가우스소거법을 통해 얻은) $E$ 가 하삼각행렬인 경우, $E^{-1} D^{-1}$ 역시 하삼각행렬이므로, 행렬 $A$ 는 다음과 같이 표현됩니다.

LDU Decomposition
상삼각 행렬의 대각성분이 모두 1이어야 하는 경우에는, 다음과 같습니다.
$LDU=(E^{-1})\tilde{D}^{-1}(\tilde{D}EA)=A$
LU Decomposition
상삼각 행렬의 대각성분이 1이 아니어도 되는 경우에는, 다음과 같습니다
$LU=(E^{-1})(EA)=A$

$PA=LDU$ 파이썬 코드

LU분해와 알고리즘과 수학적 원리가 거의 동일합니다. 단, 순열행렬 $P$ 가 알고리즘에 작용하는 부분에 대한 설명이 상당히 복잡해서, 정확한 설명은 생략하겠습니다. 그게 필요하면, 다음내용을 보면 됩니다.

Numerical linear algebra, Trefethen & Bau (1997), Lecture 21. Pivoting p. 159

import numpy as np
A = np.matrix('1 2 1; 1 2 2; 2 1 1')
b = np.matrix('-9;-10;-7')
x_answer = np.linalg.solve(A,b)
elementary_low_triangular_matrix_list = [] # It is not used in the algorithm 
eps = 1.0e-10

N = A.shape[1] # Reset N so that A is N by N matrix 
U = A.copy() # It is not yet upper triangular matrix. 
# after the following procedure U will be The upper triangular.
L = np.zeros((N,N)) # It is not yet the desired lower triangular matrix. 
# after the following procedure L will be the lower triangular matrix.
P = np.eye(N) # It is not yet the desired permutation matrix 
# after the following procedure P will be the permutation matrix 
for i in range(0,N-1):
    index_max = np. argmax (np. abs (U[i: , i]) ) +i
        
    """
      index_max can be any number in {i,i+1,...,N-1} 
     such that U[index_max, i] is non zero. Such index_max  always exists unless the matrix is singular. I recommend that set
     ind_max so that U[index_max, i] is the maximum of 
     U[i, i],U[i+1,i],... ,U[N-1,i]
    """
    if abs(U[index_max, i]) < eps:
       # It is not true unless the matrix is singular.  
       raise Exception('index max를 올바르게 설정하세요.')
    # Pi : the ith permutation matrix 
    Pi = np.eye(N)# N by N identity matrix
    if index_max != i:
        Pi[[i, index_max]] = Pi[[index_max, i]] # swap the i-row (i.e. 0-row) and pivot-row (i.e. 2-row)
    U = np.matmul(Pi,U) # swap the i-row (i.e. 0-row) and pivot-row (i.e. 2-row) 
    L = np.matmul(Pi,L)
    P = np.matmul(Pi,P)
    pivot_index = i # now, the pivot_index is changed to i. 
    the_elementary_operation = np.eye(N)# N by N identity matrix
    for k in range(pivot_index+1, N):
        l_el = U[k,pivot_index]/U[pivot_index,pivot_index] 
        the_elementary_operation[k, pivot_index] -= l_el
        L[k, pivot_index] += l_el 
    elementary_low_triangular_matrix_list.append(the_elementary_operation)
    U = np.matmul(the_elementary_operation,U) 
# end of loop
# complete L 
L+=np.eye(N)
# Decompose U by DU
diag = np.diag(U)
D = np.diag(diag)
U /= diag[:,None] #Normalize rows of U
# Check wheter P, L, U, D is deomposition of A or not
if np.allclose ( L @ D @ U , P @ A , rtol = eps ) :
    print (" L , U , D is decomposition of PA ")
else :
    print (" The decomposition results are wrong .")
### Result printing
print('P =\n')
print(P) # P 
print('L =\n')
print(L) 
 
print('U =\n')
print(U) # U
print('D =\n')
print(D) # D    


 ###########################
# Solve Linear System   #
###########################   
 # We should redefine A and b so that DUx =invL b holds.
 
A = np.matmul(P,A)
b = np.matmul(P,b)
 
 # Now, 
 # A x = LDU x = b so, DUx =invL b
 # We can solve x by 
 # x = invU invD invL b
invL = np.eye(N)
for i in range(0, N-1):
    Li = np.eye(N)
    Li[i+1:,[i]] = -L[i+1:,[i]]
    invL = np.matmul(Li,invL)

invU = np.eye(N)
for i in range(1, N):
    Ui = np.eye(N)
    Ui[:i,[i]] = -U[:i,[i]]
    invU = np.matmul(invU,Ui)
diag=np.diag(D)
invD=np.eye(N)
invD/=diag[:None]
x = invU @ invD @ invL @ b
if np.allclose ( x , x_answer , rtol = eps ) :
    print (" correct x")
else :
    print (" wrong !")

김록기

문제풀이를 즐김

이전 포스트