5. LU 분해

In [2]:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

In [3]:

# 어떤 행렬을 정의한다.

x = np.matrix([[1,5,3],[2,3,6],[3,1,7]])

5. LU분해¶

1) 기본 개념¶

(1) 기본 개념은, 각각 소블럭의 1행과 1열을 외적하여 행렬로 전개한 값을, 원래 행렬에서 지속적으로 빼서¶

(2) a11 성분으로 나누어 1행과 1열을 원래 행렬의 값과 동일하게 만들어 준 후에¶

(3) 1행과 1열을 0으로 만들어주면서 지속적으로 축소해준다.¶

2) 기본 알고리즘¶

In [4]:

# 1열에 1행을 외적할 경우, 생성되는 행렬의 1행과 1열엔 a00성분이 중복해서 들어가게 된다.
# normalizer는 a00로 나눠줌으로서 중복성을 제거해주고, 원래 행렬에서 잔차 행렬을 뺄 때 1행 1열은 0이 되도록 만들어주는 
# 핵심 역할을 담당한다.
normalizer = 1/(x[0,0])

# 1행에 normalizer를 곱해준다.
temp_l = normalizer * x[0].T
# 1열을 임시로 저장해주고
temp_u = x.T[0]
l = normalizer * x[0].T
u = x.T[0]

print(l,u)

[[1.]
 [5.]
 [3.]] [[1 2 3]]

In [5]:

# 1행과 1열을 외적하여 n*m 잔차행렬을 새로 생성해준다.
temp = np.matmul(temp_l,temp_u).T
print(temp)

[[ 1.  5.  3.]
 [ 2. 10.  6.]
 [ 3. 15.  9.]]

In [6]:

# 원래 행렬에서 잔차행렬을 빼 1행과 1열을 0으로 만들어준다.
x = x - temp
print(x)

[[  0.   0.   0.]
 [  0.  -7.   0.]
 [  0. -14.  -2.]]

In [7]:

normalizer = 1/x[1,1]
temp_l = normalizer * x[1].T
temp_u = x.T[1]

#위 iteration의 반복이나, i=2부터는 하삼각행렬과 상삼각행렬을 numpy의 stack으로 지속적으로 누적해준다.
l = np.hstack([l,temp_l])
u = np.vstack([u,temp_u])
temp = np.matmul(temp_l,temp_u).T
x = x - temp

In [8]:

normalizer = 1/x[2,2]
temp_l = normalizer * x[2].T
temp_u = x.T[2]
l = np.hstack([l,temp_l])
u = np.vstack([u,temp_u])
temp = np.matmul(temp_l,temp_u).T
x = x - temp
# 최종적인 잔차행렬이 0행렬이 되면 성공이다.
print(x)

[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

In [9]:

# 함수로 정의한 LU분해
def LUdecomp(x) :
    if type(x) != type(np.matrix([0])):
        raise ValueError 
    
    #행렬의 행과 열 차원중 작은쪽을 s로 지정해준다.
    s = np.min(x.shape)
    
    normalizer = 1/(x[0,0])
    l = normalizer * x[0].T
    u = x.T[0]
    temp = np.matmul(l,u).T
    x = x - temp
    
    for i in range(1,s):
        normalizer = 1/x[i,i]
        temp_l = normalizer * x[i].T
        temp_u = x.T[i]
        l = np.hstack([l,temp_l])
        u = np.vstack([u,temp_u])
        temp = np.matmul(temp_l,temp_u).T
        x = x - temp
    return (l,u)

In [10]:

x = np.matrix([[1,5,3],[2,3,6],[3,1,7]])

L,U = LUdecomp(x)

print("하삼각행렬")
print(L)
print("상삼각행렬")
print(U)

하삼각행렬
[[ 1. -0. -0.]
 [ 5.  1. -0.]
 [ 3. -0.  1.]]
상삼각행렬
[[  1.   2.   3.]
 [  0.  -7. -14.]
 [  0.   0.  -2.]]

3) a00 = 0일때의 대처법(피보팅)¶

(1) a00 = 0일때의 문제는, 1/a00를 곱하는 과정에서 1/0이 되기 때문에 연산이 불가능해지는 점이다.¶

(2) 이를 해결하기 위해, 0이 아닌 다른 행으로 교체한 후, 앞전 연산의 결과도 동일한 순서로 행교환을 해준다.¶

(3) 알고리즘적으로는, 치환된 순서를 규정하는 '치환행렬'을 만들어 이를 계산에 활용한다.¶

• Ax = y에서 LUx = y로 분해했을 경우, 치환행렬 P를 이용하여 이를 다시 나타내면 PLUx = y,
• 이때, P는 정규직교행렬이므로 inverse(P) = transpose(P)가 성립되고. 우변으로 이항하면 LUx = transpose(T) * y 즉 y쪽의 순서를 변경한다.
• Sxy를 i번째 iteration의 x행과 y행의 교환을 나타내는 치환행렬이라고 할 때
• P = Sxy1 Sxy2 ..... 이다.

In [125]:

A = np.matrix([[0,5,3,1],[2,3,6,5],[3,1,7,3],[2,1,5,6]])
print(x)

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

In [126]:

LUdecomp(A) 

C:\Users\never\Anaconda3\lib\site-packages\ipykernel_launcher.py:9: RuntimeWarning: divide by zero encountered in long_scalars
  if __name__ == '__main__':
C:\Users\never\Anaconda3\lib\site-packages\ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in multiply
  # Remove the CWD from sys.path while we load stuff.
C:\Users\never\Anaconda3\lib\site-packages\ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in multiply

Out[126]:

(matrix([[nan, nan, nan, nan],
         [inf, nan, nan, nan],
         [inf, nan, nan, nan],
         [inf, nan, nan, nan]]), matrix([[  0.,   2.,   3.,   2.],
         [ nan, -inf, -inf, -inf],
         [ nan,  nan,  nan,  nan],
         [ nan,  nan,  nan,  nan]]))

• a00성분이 0이기 때문에, 무한으로 발산(inf)한다.

(4) 핵심요소는 다음과 같다.¶

• 대각성분이 1인 피봇행렬을 생성한다.

In [127]:

# 그리고, 행 교환을 원하는 행의 '피봇 행렬의 행'을 서로 바꿔준다. 아래에서 그 예를 확인할 수 있다.

print(A)

pivot_table = np.matrix([[0,1,0,0],[1,0,0,0],[0,0,1,0],[0,0,0,1]]) #1행과 2행의 순서가 바뀐 pivottable을 만든다.

print(np.dot(pivot_table,x))# 2행과 1행의 순서가 바뀐것을 확인할 수 있다.

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

• 0이 아닌 행을 찾는다.

In [128]:

# X의 전치 행렬에서, 0이 아닌 가장 가까운 행을 np.where로 찾는다.
ind = np.where(A.T[0] != 0)[1][0]
print(ind) # 해당 잔차블록행렬의 첫 번째 원소가 0이 아님을 보여준다.

1

In [129]:

def making_pivot(x):
    rows = np.shape(x)[0]
    pivot_table = np.diag([1 for i in range(0,rows)])
    return pivot_table

In [130]:

def pivoting(pivot_table,where,ind):
    temp = pivot_table[where].copy()
    pivot_table[where] = pivot_table[ind].copy()
    pivot_table[ind] = temp.copy()
    return pivot_table

In [131]:

# 함수로 정의한 LU분해
def LUdecomp_pivot(x) :
    if type(x) != type(np.matrix([0])):
        raise ValueError 
    
    #행렬의 행과 열 차원중 작은쪽을 s로 지정해준다.
    s = np.min(x.shape)
    
    # X의 전치행렬에서 0이 아닌 가장 가까운 행을 np.where로 찾는다.
    ind = np.where(x.T[0] != 0)[1][0]
    # 행의 순서 교환을 기록한 pivot_table을 생성한다.
    pivot_table = pivoting(making_pivot(x),0,ind)
    x = np.matmul(pivot_table,x)
    pivot_table_summary = pivot_table.copy()
    
    normalizer = 1/(x[0,0])
    l = normalizer * x[0].T
    u = x.T[0]
    temp = np.matmul(l,u).T
    x = np.matrix(np.round(x - temp,decimals = 5))
    
    for i in range(1,s):
        # X의 전치행렬에서 0이 아닌 가장 가까운 행을 np.where로 찾는다.
        ind = np.where(x.T[i] != 0)[1][0]
        # 행의 순서 교환을 기록한 pivot_table을 생성한다.
        pivot_table = pivoting(making_pivot(x),i,ind)
        x = np.matmul(pivot_table,x)
        pivot_table_summary = np.matmul(pivot_table_summary,pivot_table)
        
        normalizer = 1/x[i,i]
        temp_l = normalizer * x[i].T
        temp_u = x.T[i]
        l = np.hstack([l,temp_l])
        u = np.vstack([u,temp_u])
        temp = np.matmul(temp_l,temp_u).T
        x = np.matrix(np.round(x - temp,decimals = 5))
    return (l,u,pivot_table_summary)

In [132]:

l,u,pivot = LUdecomp_pivot(A)

In [133]:

A

Out[133]:

matrix([[0, 5, 3, 1],
        [2, 3, 6, 5],
        [3, 1, 7, 3],
        [2, 1, 5, 6]])

In [134]:

LU = np.matmul(l,u).T

피보팅된 LU분해가 끝난후, 원 행렬과 비교하면 원래와는 달리 순서가 뒤바뀌어 있음을 확인할 수 있다.¶

In [135]:

np.matmul(pivot,LU)

Out[135]:

matrix([[0., 5., 3., 1.],
        [2., 3., 6., 5.],
        [3., 1., 7., 3.],
        [2., 1., 5., 6.]])

피보팅 행렬과 LU분해된 행렬을 곱하면, 원래 행렬이 비로소 복원됨을 확인할 수 있다.¶

4) LU분해된 행렬을 이용해 해 x를 구하기¶

(1) 개요¶

• 순서는 다음과 같다.

Ax = y일 때 : pivot행렬 (L U)^T * x = y 인 상황이다.

열공간에 속하는 벡터 y에, pivot행렬의 역행렬을 곱해 알맞는 순서로 바꾸어준다.

(U^T L^T)x = pivot^(-1) y = y'

- U^T(L^T)x = y' 이므로, (L^T)x = z로 놓으면, 이 방정식은 U^T * Z = y'가 된다. 우선, 중간 단계인 z를 구하는 것을 우선으로 한다.¶

In [114]:

def sol_U(U,y,pivot,debug):
    # U를 전치해준다.
    U = U.T
    # 정답벡터 y에, pivot을 우변으로 이항하여(즉, 역행렬로 만들어) y에 행렬곱해준다.
    y_pivot = np.matmul(np.linalg.inv(pivot),y)
    # 행렬의 모양을 가져온다.
    shape = U.shape[0]
    # 차원 만큼 0로 차있는, 임시 벡터 0를 만들어준다.
    z = np.zeros(shape)
    # 초깃값 z[0]는 정답벡터 y'를 상삼각행렬의 첫번째 대각성분으로 나눠준 값으로 한다.
    z[0] = y_pivot[0]/U[0,0]

    if debug:
        print(z)

    for i in range(1,shape):
        # 여기가 핵심이다.
        # 노리는 효과는, z는 기본적으로 i-1번째 요소까지만 값이 저장되어있고, 나머지는 0인 상태를 이용하는 것이다.
        # 상삼각행렬 U의 i번째 행과 i-1번째 요소까지만 값이 저장되어있는 z를 행렬곱하면, 구하고자 하는 성분인 i번째 열에 속하는 
        # 변수 z[i]는 0으로 처리되어 계산에 전혀 영향을 주지 않는다.(디버그 모드를 키면 무슨 말인지 이해할 수 있다)
        # z[i]에 속하는 대각성분인 U[i,i]를 계산한 값에 나누기 해주면, z[i]를 구할 수 있다.
        z[i] = (y_pivot[i] - np.matmul(U[i],z))/U[i,i]
        if debug:
            print(z)
    if debug:
        print("상삼각행렬 연산 종료")
    return z

- 중간단계 벡터인 z를 구했으므로, 이제 (L^T)x = z가 되도록 하는 x를 구하면 연산이 완료된다.¶

In [117]:

def sol_L(L,z,debug):
    L = L.T
    shape = L.shape[0] - 1
    x = np.zeros(shape + 1)
    x[shape] = z[shape]/L[shape,shape]
    if debug:
        print(x)
    for i in reversed(range(0,shape)):
        x[i] = z[i] - np.matmul(L[i],x)/L[i,i]
        if debug:
            print(x)
    if debug:
        print("하삼각행렬 연산 종료")
    return x

In [118]:

def solving_pivot(L,U,pivot,y,debug = False):
    z = sol_U(U,y,pivot,debug)
    x = sol_L(L,z,debug)
    return x

In [139]:

y = np.array([2,3,4,1])
x = solving_pivot(l,u,pivot,y,debug=True)
print(x)

[1.5 0.  0.  0. ]
[1.5 0.4 0.  0. ]
[1.5 0.4 9.  0. ]
[ 1.5         0.4         9.         -0.33333333]
상삼각행렬 연산 종료
[ 0.          0.          0.         -0.33333333]
[ 0.          0.         -3.66666667 -0.33333333]
[ 0.          2.66666667 -3.66666667 -0.33333333]
[ 9.33333333  2.66666667 -3.66666667 -0.33333333]
하삼각행렬 연산 종료
[ 9.33333333  2.66666667 -3.66666667 -0.33333333]

(2) 검산¶

In [141]:

# 도출된 x와 원래 벡터 A를 내적하면

np.matmul(A,x)

# 원래 의도했던 y = [2,3,4,1]이 정상적으로 출력되는 것을 볼 수 있다.

Out[141]:

matrix([[2., 3., 4., 1.]])

(3) LU분해에서 계산까지¶

In [142]:

def solver(A,y,debug=False):
    L,U,pivot = LUdecomp_pivot(A)
    x = solving_pivot(L,U,pivot,y)
    return x

In [143]:

solver(A,y)

Out[143]:

array([ 9.33333333,  2.66666667, -3.66666667, -0.33333333])