Lecture 4. A = LU 분해 — 파이썬 실습 - 선형대수

L4 서술 파트에서 $A = LU$ 분해를 다루었다. 이 노트북에서는 소거하면서 곱수를 기록해 $L$ 을 직접 만들어 보고, 곱하는 순서에 따라 간섭이 생기는지 아닌지를 숫자로 확인해 보도록 하자.

서술 파트의 내용	여기서 확인하는 방법
$L$ 에는 곱수가 그대로 들어간다	소거 함수가 곱수를 `L[i, j]` 에 적게 고친다
곱하는 순서에 따라 간섭이 생긴다	$E_{32}E_{21}$ 과 $E_{21}^{-1}E_{32}^{-1}$ 의 $(3,1)$ 성분을 비교한다
간섭 조건은 $\vv{e}_j^{\mathsf{T}}\vv{e}_k = 1$	외적으로 소거 행렬을 만들어 확인한다
분해를 재사용하면 빠르다	우변 100개를 두 방식으로 풀어 시간을 잰다
$PA = LU$	행 교환이 필요한 행렬에 `scipy` 를 써 본다

0. 준비¶

import numpy as np

from linalg_viz import show_matrix

np.set_printoptions(precision=3, suppress=True)
print("numpy", np.__version__)

numpy 2.5.2

A = np.array([[1, 2, 1],
              [3, 8, 1],
              [0, 4, 1]], dtype=float)
b = np.array([2, 12, 2], dtype=float)

print(show_matrix(A, "A ="))

A =
[  1   2   1 ]
[  3   8   1 ]
[  0   4   1 ]

1. 소거하면서 $L$ 을 함께 기록하기¶

L2에서 만든 소거 함수와 거의 같다. 달라지는 것은 한 줄뿐이다. 곱수를 계산해서 쓰고 버리는 대신, L[i, j] 에 적어 둔다.

def lu_decompose(A):
    """소거하면서 곱수를 L 에 기록한다. 행 교환은 하지 않는다."""
    U = np.asarray(A, dtype=float).copy()
    n = U.shape[0]
    L = np.eye(n)

    for j in range(n - 1):
        if abs(U[j, j]) < 1e-12:
            raise ValueError(f"{j + 1}번째 피벗이 0이다. 행 교환이 필요하다.")
        for i in range(j + 1, n):
            m = U[i, j] / U[j, j]          # 곱수
            U[i] = U[i] - m * U[j]         # 소거
            L[i, j] = m                    # 버리지 않고 적어 둔다
    return L, U

L, U = lu_decompose(A)

print(show_matrix(L, "L ="))
print(show_matrix(U, "U ="))
print()
print("곱수  m21 =", L[1, 0], "  m31 =", L[2, 0], "  m32 =", L[2, 1])
print("피벗 :", np.diag(U))

L =
[  1   0   0 ]
[  3   1   0 ]
[  0   2   1 ]
U =
[   1    2    1 ]
[   0    2   -2 ]
[   0    0    5 ]

곱수  m21 = 3.0   m31 = 0.0   m32 = 2.0
피벗 : [1. 2. 5.]

서술 파트에서 손으로 구한 것과 같다. $L$ 의 대각선 아래가 곱수 $3, 0, 2$ 이고 $U$ 의 대각선이 피벗 $1, 2, 5$ 이다.

print(show_matrix(L @ U, "L U ="))
print("A 와 같은가 :", np.allclose(L @ U, A))

L U =
[  1   2   1 ]
[  3   8   1 ]
[  0   4   1 ]
A 와 같은가 : True

2. 곱하는 순서에 따라 간섭이 생긴다¶

서술 파트의 핵심이다. 같은 소거 행렬들인데 곱하는 순서에 따라 없던 성분이 생기기도 하고 안 생기기도 한다. 네 가지 경우를 모두 계산해 비교해 보자.

I3 = np.eye(3)

E21 = I3.copy(); E21[1, 0] = -3
E32 = I3.copy(); E32[2, 1] = -2

E21_inv = I3.copy(); E21_inv[1, 0] = 3      # 부호만 바꾼 것
E32_inv = I3.copy(); E32_inv[2, 1] = 2

print("E21 E21^-1 가 I 인가 :", np.allclose(E21 @ E21_inv, I3))
print("E32 E32^-1 가 I 인가 :", np.allclose(E32 @ E32_inv, I3))

E21 E21^-1 가 I 인가 : True
E32 E32^-1 가 I 인가 : True

경우 = [
    ("E32 E21          (소거 순서)",        E32 @ E21),
    ("E21 E32          (순서를 바꾸면)",     E21 @ E32),
    ("E21^-1 E32^-1    (L 을 만드는 순서)",  E21_inv @ E32_inv),
    ("E32^-1 E21^-1    (순서를 바꾸면)",     E32_inv @ E21_inv),
]

for 이름, M in 경우:
    print(f"{이름}   ->  (3,1) 성분 = {M[2, 0]:>4.0f}")
    print(show_matrix(M))
    print()

E32 E21          (소거 순서)   ->  (3,1) 성분 =    6
[   1    0    0 ]
[  -3    1    0 ]
[   6   -2    1 ]

E21 E32          (순서를 바꾸면)   ->  (3,1) 성분 =    0
[   1    0    0 ]
[  -3    1    0 ]
[   0   -2    1 ]

E21^-1 E32^-1    (L 을 만드는 순서)   ->  (3,1) 성분 =    0
[  1   0   0 ]
[  3   1   0 ]
[  0   2   1 ]

E32^-1 E21^-1    (순서를 바꾸면)   ->  (3,1) 성분 =    6
[  1   0   0 ]
[  3   1   0 ]
[  6   2   1 ]

네 경우 중 두 개에서 $(3,1)$ 자리에 $\pm 6$ 이 나타났다. 서술 파트에서 본 대로, 왼쪽 행렬의 열 번호와 오른쪽 행렬의 행 번호가 같을 때만 간섭이 생긴다.

곱	왼쪽 $(i,j)$	오른쪽 $(k,l)$	$j = k$ 인가	$(3,1)$
$E_{32}E_{21}$	$(3,2)$	$(2,1)$	예	6
$E_{21}E_{32}$	$(2,1)$	$(3,2)$	아니오	0
$E_{21}^{-1}E_{32}^{-1}$	$(2,1)$	$(3,2)$	아니오	0
$E_{32}^{-1}E_{21}^{-1}$	$(3,2)$	$(2,1)$	예	6

$L$ 을 만들 때 쓰는 순서가 간섭이 없는 쪽이라는 것을 확인할 수 있다.

print("E21^-1 E32^-1 가 lu_decompose 의 L 과 같은가 :",
      np.allclose(E21_inv @ E32_inv, L))

E21^-1 E32^-1 가 lu_decompose 의 L 과 같은가 : True

3. 외적으로 확인하기¶

서술 파트에서 소거 행렬을 $E_{ij} = I - m\,\vv{e}_i\vv{e}_j^{\mathsf{T}}$ 로 썼다. 간섭 항이 살아남는 조건은 가운데에 생기는 내적 $\vv{e}_j^{\mathsf{T}}\vv{e}_k$ 가 1이 되는 것, 즉 $j = k$ 였다. 직접 만들어 확인해 보자.

def e(i, n=3):
    """i 번째 표준기저벡터. 파이썬 관례대로 0부터 센다."""
    v = np.zeros(n)
    v[i] = 1.0
    return v

print(show_matrix(np.outer(e(2), e(0)), "e3 e1^T ="))
print()
print("e1 . e3 =", e(0) @ e(2), "  -> 0 이므로 간섭 항이 사라진다")
print("e2 . e2 =", e(1) @ e(1), "  -> 1 이므로 간섭 항이 살아난다")

e3 e1^T =
[  0   0   0 ]
[  0   0   0 ]
[  1   0   0 ]

e1 . e3 = 0.0   -> 0 이므로 간섭 항이 사라진다
e2 . e2 = 1.0   -> 1 이므로 간섭 항이 살아난다

소거 행렬도 이 형태로 만들 수 있다.

E21_outer = I3 - 3 * np.outer(e(1), e(0))       # I - m e2 e1^T
E32_outer = I3 - 2 * np.outer(e(2), e(1))       # I - m e3 e2^T

print("E21 과 같은가 :", np.allclose(E21_outer, E21))
print("E32 와 같은가 :", np.allclose(E32_outer, E32))

E21 과 같은가 : True
E32 와 같은가 : True

그러면 $L$ 은 간섭 항 없이 두 외적을 더한 것이어야 한다.

L_outer = I3 + 3 * np.outer(e(1), e(0)) + 2 * np.outer(e(2), e(1))

print(show_matrix(L_outer, "I + 3 e2 e1^T + 2 e3 e2^T ="))
print("L 과 같은가 :", np.allclose(L_outer, L))

I + 3 e2 e1^T + 2 e3 e2^T =
[  1   0   0 ]
[  3   1   0 ]
[  0   2   1 ]
L 과 같은가 : True

간섭이 생기는 쪽도 식으로 재현해 보자. $E_{32}E_{21}$ 은 두 항에 더해 $(-2)(-3)\,\vv{e}_3\vv{e}_1^{\mathsf{T}}$ 가 붙어야 한다.

예상 = I3 - 3 * np.outer(e(1), e(0)) - 2 * np.outer(e(2), e(1)) \
       + (-2) * (-3) * np.outer(e(2), e(0))

print(show_matrix(예상, "간섭 항까지 넣어 예상한 값 ="))
print("E32 @ E21 과 같은가 :", np.allclose(예상, E32 @ E21))

간섭 항까지 넣어 예상한 값 =
[   1    0    0 ]
[  -3    1    0 ]
[   6   -2    1 ]
E32 @ E21 과 같은가 : True

4. $LDU$ 분해¶

$U$ 의 각 행을 그 행의 피벗으로 나누면 대각이 1인 위삼각행렬이 남고, 피벗은 대각행렬 $D$ 로 모인다.

pivots = np.diag(U)
D = np.diag(pivots)
U_unit = U / pivots[:, None]          # 각 행을 그 행의 피벗으로 나눈다

print(show_matrix(D, "D ="))
print(show_matrix(U_unit, "U' ="))
print()
print("D U' 가 U 와 같은가   :", np.allclose(D @ U_unit, U))
print("L D U' 가 A 와 같은가 :", np.allclose(L @ D @ U_unit, A))

D =
[  1   0   0 ]
[  0   2   0 ]
[  0   0   5 ]
U' =
[   1    2    1 ]
[   0    1   -1 ]
[   0    0    1 ]

D U' 가 U 와 같은가   : True
L D U' 가 A 와 같은가 : True

5. 두 번의 대입으로 풀기¶

$A\vv{x} = \vv{b}$ 를 $L\vv{c} = \vv{b}$ 와 $U\vv{x} = \vv{c}$ 로 나누어 푼다. $L$ 은 하삼각이므로 위에서 아래로, $U$ 는 위삼각이므로 아래에서 위로 하나씩 확정된다.

def forward_substitute(L, b):
    """하삼각 시스템 L c = b 를 첫 식부터 차례로 푼다."""
    n = len(b)
    c = np.zeros(n)
    for i in range(n):
        이미_구한_항 = L[i, :i] @ c[:i]
        c[i] = (b[i] - 이미_구한_항) / L[i, i]
    return c


def back_substitute(U, c):
    """위삼각 시스템 U x = c 를 마지막 식부터 거꾸로 푼다."""
    n = len(c)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        이미_구한_항 = U[i, i + 1:] @ x[i + 1:]
        x[i] = (c[i] - 이미_구한_항) / U[i, i]
    return x

c = forward_substitute(L, b)
x = back_substitute(U, c)

print("b =", b)
print("c =", c, "  <- L c = b 를 풀어 얻은 값")
print("x =", x, "  <- U x = c 를 풀어 얻은 값")
print()
print("np.linalg.solve :", np.linalg.solve(A, b))
print("같은가 :", np.allclose(x, np.linalg.solve(A, b)))

b = [ 2. 12.  2.]
c = [  2.   6. -10.]   <- L c = b 를 풀어 얻은 값
x = [ 2.  1. -2.]   <- U x = c 를 풀어 얻은 값

np.linalg.solve : [ 2.  1. -2.]
같은가 : True

여기서 나온 $\vv{c} = (2, 6, -10)$ 은 L2에서 증강행렬을 소거했을 때 우변이 도달한 값과 같다. 소거가 우변에 하던 일이 곧 $L\vv{c} = \vv{b}$ 를 푸는 일이었다.

6. 분해를 재사용하면 얼마나 빨라지는가¶

서술 파트에서 분해는 $n^3/3$ , 우변 하나당 대입은 $n^2$ 이라고 하였다. 우변이 여러 개일 때 차이가 얼마나 나는지 직접 재 보자.

scipy 의 lu_factor 는 분해를 한 번 해서 결과를 들고 있고, lu_solve 는 그 결과로 대입만 한다.

import time

from scipy.linalg import lu_factor, lu_solve

rng = np.random.default_rng(0)
n, k = 200, 50
A_big = rng.standard_normal((n, n))
B_many = rng.standard_normal((n, k))        # 우변 50개

시작 = time.perf_counter()
X1 = np.column_stack([np.linalg.solve(A_big, B_many[:, j]) for j in range(k)])
시간_매번 = time.perf_counter() - 시작

시작 = time.perf_counter()
분해 = lu_factor(A_big)                                        # 한 번만
X2 = np.column_stack([lu_solve(분해, B_many[:, j]) for j in range(k)])
시간_재사용 = time.perf_counter() - 시작

print(f"크기 n = {n}, 우변 {k}개")
print(f"  이론적으로 예상되는 배수 : 약 {n / 3:.0f}")
print(f"  매번 처음부터 solve : {시간_매번:.3f} 초")
print(f"  한 번 분해 후 대입   : {시간_재사용:.3f} 초")
print(f"  배수                : {시간_매번 / 시간_재사용:.1f}")
print()
print("두 결과가 같은가 :", np.allclose(X1, X2))

크기 n = 200, 우변 50개
  이론적으로 예상되는 배수 : 약 67
  매번 처음부터 solve : 0.020 초
  한 번 분해 후 대입   : 0.004 초
  배수                : 5.2

두 결과가 같은가 : True

측정한 시간은 실행할 때마다 다르고 컴퓨터에 따라서도 다르다. 중요한 것은 절대값이 아니라 우변의 개수가 늘어날 때 두 방식이 다르게 늘어난다는 점이다.

실제로 재 보면 연산량으로 예상한 $n/3$ 배보다 차이가 더 크게 나온다. 연산 횟수 말고도 함수를 호출할 때마다 붙는 비용이 있기 때문이다. np.linalg.solve 는 부를 때마다 행렬을 복사하고 LAPACK 루틴을 준비하는데, lu_solve 는 이미 분해된 결과를 받아 대입만 한다. 연산량 계산은 큰 그림을 알려 주는 것이고, 실제 속도는 늘 그 위에 구현 비용이 얹힌다.

7. 행 교환이 필요한 경우 : $PA = LU$ ¶

서술 파트 6절의 행렬이다. 소거 도중 두 번째 피벗이 0이 되어 그대로는 분해되지 않는다.

A_swap = np.array([[1, 2, 1],
                   [3, 6, 1],
                   [0, 4, 1]], dtype=float)

try:
    lu_decompose(A_swap)
except ValueError as err:
    print("분해 실패 :", err)

분해 실패 : 2번째 피벗이 0이다. 행 교환이 필요하다.

2행과 3행을 미리 바꿔 두면 소거가 끝까지 진행된다.

P = np.array([[1, 0, 0],
              [0, 0, 1],
              [0, 1, 0]], dtype=float)

L_p, U_p = lu_decompose(P @ A_swap)

print(show_matrix(P @ A_swap, "P A ="))
print(show_matrix(L_p, "L ="))
print(show_matrix(U_p, "U ="))
print("L U 가 P A 와 같은가 :", np.allclose(L_p @ U_p, P @ A_swap))

P A =
[  1   2   1 ]
[  0   4   1 ]
[  3   6   1 ]
L =
[  1   0   0 ]
[  0   1   0 ]
[  3   0   1 ]
U =
[   1    2    1 ]
[   0    4    1 ]
[   0    0   -2 ]
L U 가 P A 와 같은가 : True

`scipy` 는 어떻게 하는가¶

scipy.linalg.lu 는 부분 피벗팅을 쓰므로 우리가 고른 것과 다른 $P$ 를 고를 수 있다. 각 열에서 절댓값이 가장 큰 성분을 피벗 자리로 올리기 때문이다.

표기 관례도 다르다. scipy 는 $PA = LU$ 가 아니라 $A = PLU$ 형태로 돌려준다. 치환행렬은 $P^{-1} = P^{\mathsf{T}}$ 이므로 양변에 $P^{\mathsf{T}}$ 를 곱하면 서로 옮겨 갈 수 있다.

from scipy.linalg import lu

P_sp, L_sp, U_sp = lu(A_swap)

print(show_matrix(P_sp, "scipy 의 P ="))
print(show_matrix(L_sp, "scipy 의 L ="))
print(show_matrix(U_sp, "scipy 의 U ="))
print()
print("A = P L U 인가        :", np.allclose(P_sp @ L_sp @ U_sp, A_swap))
print("P^T A = L U 인가      :", np.allclose(P_sp.T @ A_swap, L_sp @ U_sp))
print("P^T P 가 I 인가       :", np.allclose(P_sp.T @ P_sp, np.eye(3)))

scipy 의 P =
[  0   0   1 ]
[  1   0   0 ]
[  0   1   0 ]
scipy 의 L =
[      1       0       0 ]
[      0       1       0 ]
[  0.333       0       1 ]
scipy 의 U =
[      3       6       1 ]
[      0       4       1 ]
[      0       0   0.667 ]

A = P L U 인가        : True
P^T A = L U 인가      : True
P^T P 가 I 인가       : True

scipy 가 고른 순서를 확인해 보자. $P^{\mathsf{T}}A$ 의 각 행이 원래 $A$ 의 몇 번째 행인지 보면 된다.

print(show_matrix(P_sp.T @ A_swap, "P^T A ="))
print()
for i, 행 in enumerate(P_sp.T @ A_swap):
    원래 = int(np.flatnonzero([np.allclose(행, r) for r in A_swap])[0])
    print(f"  P^T A 의 {i + 1}행  =  A 의 {원래 + 1}행")

P^T A =
[  3   6   1 ]
[  0   4   1 ]
[  1   2   1 ]

  P^T A 의 1행  =  A 의 2행
  P^T A 의 2행  =  A 의 3행
  P^T A 의 3행  =  A 의 1행

첫 열에서 절댓값이 가장 큰 값이 3 이므로 그 행을 맨 위로 올렸고, 그 결과 행 순서가 $(2, 3, 1)$ 이 되었다. 우리가 고른 $P$ 는 $(1, 3, 2)$ 였다. 둘은 다르지만 모두 옳은 분해이다. 행 교환을 어떻게 고르든 결과가 위삼각이 되기만 하면 된다.

scipy 의 $U$ 를 보면 피벗이 $3, 4, 0.667$ 로 우리 것과 다르다. 어느 행을 피벗으로 쓰느냐가 달라졌기 때문이다. 그래도 두 분해 모두 같은 $A$ 를 복원한다.

마치며...¶

서술 파트의 내용	이 노트북의 코드
곱수를 기록해 $L$ 을 만든다	`lu_decompose` 의 `L[i, j] = m` 한 줄
순서에 따라 간섭이 생긴다	2절의 네 가지 곱 비교
간섭 조건은 $\vv{e}_j^{\mathsf{T}}\vv{e}_k = 1$	3절의 `e(0) @ e(2)` 와 `e(1) @ e(1)`
$L = I + \sum m_{ij}\vv{e}_i\vv{e}_j^{\mathsf{T}}$	`L_outer`
$A = LDU'$	`U / pivots[:, None]`
두 번의 대입	`forward_substitute`, `back_substitute`
분해의 재사용	`lu_factor` 한 번 + `lu_solve` 여러 번
$PA = LU$	7절

더 해 볼 것¶

lu_decompose 가 곱수를 $U$ 의 0이 된 자리에 바로 적도록 고쳐 보자. 어차피 그 자리는 0이니 $L$ 과 $U$ 를 행렬 하나에 담을 수 있다. 실제 라이브러리가 이렇게 한다.
무작위 행렬 여러 개에 대해 lu_decompose 를 돌리고 np.diag(U) 를 모두 곱해 보자. 그 값과 np.linalg.det(A) 는 어떤 관계인가?
대칭행렬 $A = A^{\mathsf{T}}$ 를 하나 만들어 $LDU'$ 로 분해해 보자. $U'$ 과 $L^{\mathsf{T}}$ 를 비교하면 무엇이 보이는가?
6절에서 np.linalg.solve(A_big, B_many) 처럼 우변을 통째로 넘기면 시간이 어떻게 되는가? numpy 가 내부에서 무엇을 하고 있는지 짐작해 보자.

다음 강의에서는 전치와 순열을 정리한 뒤 벡터공간(vector space)을 다룬다.