Lecture 2. 소거법과 행렬 — 파이썬 실습 - 선형대수

L2 서술 파트에서 소거법을 다루었다. 이 노트북에서는 소거를 직접 구현해 보고, 소거의 각 단계가 정말 행렬 곱과 같은지 확인해 보도록 하자.

확인할 것은 다음 네 가지이다.

서술 파트의 내용	여기서 확인하는 방법
소거는 $A$ 를 $U$ 로 바꾼다	증강행렬을 한 단계씩 직접 바꿔 본다
후진 대입으로 해를 읽는다	소거와 후진 대입을 함수로 짜서 `np.linalg.solve` 와 대조한다
소거 한 단계 $=$ 행렬 곱 $E_{ij}$	$E_{21}A$ 를 계산해 소거 결과와 비교한다
소거는 해를 바꾸지 않는다	평면을 단계별로 움직이며 교점을 본다

0. 준비¶

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, plane, layout3d, slider_figure, show_matrix

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

numpy 2.5.2

서술 파트에서 쓴 시스템이다.

x + 2y + z = 2, \qquad 3x + 8y + z = 12, \qquad 4y + z = 2

(1)

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 ="))
print("b =", b)

A =
[  1   2   1 ]
[  3   8   1 ]
[  0   4   1 ]
b = [ 2. 12.  2.]

1. 증강행렬로 한 단계씩 소거하기¶

계수행렬만 다루면 우변을 따로 관리해야 하므로, $A$ 와 $\vv{b}$ 를 붙인 증강행렬 $[A \mid \vv{b}]$ 로 시작한다. np.hstack 으로 붙이면 된다.

M = np.hstack([A, b.reshape(-1, 1)])       # [A | b]
print(show_matrix(M, "[A | b] ="))

[A | b] =
[   1    2    1    2 ]
[   3    8    1   12 ]
[   0    4    1    2 ]

첫 번째 피벗은 M[0, 0] 이다. 두 번째 행의 첫 성분을 0으로 만들려면 M[1, 0] / M[0, 0] 배만큼 첫 번째 행을 빼면 된다. 이 배수를 곱수(multiplier)라 한다.

M1 = M.copy()
multiplier = M1[1, 0] / M1[0, 0]           # 3 / 1 = 3
M1[1] = M1[1] - multiplier * M1[0]

print("multiplier =", multiplier)
print(show_matrix(M1, "row2 <- row2 - 3(row1)"))

multiplier = 3.0
row2 <- row2 - 3(row1)
[   1    2    1    2 ]
[   0    2   -2    6 ]
[   0    4    1    2 ]

세 번째 행의 첫 성분은 이미 0이므로 곱수가 0이 되어 아무것도 바뀌지 않는다. 바로 두 번째 열로 넘어간다.

M2 = M1.copy()
multiplier = M2[2, 1] / M2[1, 1]           # 4 / 2 = 2
M2[2] = M2[2] - multiplier * M2[1]

print("multiplier =", multiplier)
print(show_matrix(M2, "row3 <- row3 - 2(row2)"))

multiplier = 2.0
row3 <- row3 - 2(row2)
[    1     2     1     2 ]
[    0     2    -2     6 ]
[    0     0     5   -10 ]

U, c = M2[:, :-1], M2[:, -1]
print(show_matrix(U, "U ="))
print("c      =", c)
print("pivots =", np.diag(U))

U =
[   1    2    1 ]
[   0    2   -2 ]
[   0    0    5 ]
c      = [  2.   6. -10.]
pivots = [1. 2. 5.]

서술 파트에서 손으로 구한 것과 같다. 피벗은 $1, 2, 5$ 이다.

2. 소거를 함수로 만들기¶

위에서 한 일을 일반적인 크기의 행렬에 대해 반복하면 된다. 열을 왼쪽부터 훑으면서, 각 열의 피벗 아래에 있는 행들을 하나씩 정리한다.

나중에 그림을 그릴 때 쓰려고 각 단계의 증강행렬도 함께 기록해 둔다.

def eliminate(A, b):
    """소거를 수행해 (U, c, history) 를 돌려준다. 행 교환은 하지 않는다."""
    M = np.hstack([np.asarray(A, float), np.asarray(b, float).reshape(-1, 1)])
    n = M.shape[0]
    history = [(M.copy(), "start")]

    for col in range(n - 1):
        pivot = M[col, col]
        if abs(pivot) < 1e-12:
            raise ValueError(f"{col + 1}번째 피벗이 0이다. 행 교환이 필요하다.")
        for row in range(col + 1, n):
            m = M[row, col] / pivot            # 곱수
            M[row] = M[row] - m * M[col]
            history.append((M.copy(), f"row{row + 1} - {m:g}(row{col + 1})"))

    return M[:, :-1], M[:, -1], history

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

U, c, history = eliminate(A, b)

print(show_matrix(U, "U ="))
print("c =", c)
print()

x = back_substitute(U, c)
print("back_substitute   :", x)
print("np.linalg.solve   :", np.linalg.solve(A, b))
print("같은가 :", np.allclose(x, np.linalg.solve(A, b)))

U =
[   1    2    1 ]
[   0    2   -2 ]
[   0    0    5 ]
c = [  2.   6. -10.]

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

history 에는 시작 상태와 각 소거 단계가 들어 있다. 곱수가 0이어서 아무것도 바뀌지 않은 단계도 그대로 기록된다.

for M_step, 설명 in history:
    print(f"[{설명}]")
    print(show_matrix(M_step))
    print()

[start]
[   1    2    1    2 ]
[   3    8    1   12 ]
[   0    4    1    2 ]

[row2 - 3(row1)]
[   1    2    1    2 ]
[   0    2   -2    6 ]
[   0    4    1    2 ]

[row3 - 0(row1)]
[   1    2    1    2 ]
[   0    2   -2    6 ]
[   0    4    1    2 ]

[row3 - 2(row2)]
[    1     2     1     2 ]
[    0     2    -2     6 ]
[    0     0     5   -10 ]

3. 소거 한 단계는 행렬 곱이다¶

서술 파트의 핵심이다. "2행에서 1행의 3배를 뺀다"는 동작은 단위행렬의 $(2,1)$ 자리를 -3 으로 바꾼 행렬 $E_{21}$ 과 같다.

I = np.eye(3)

E21 = I.copy()
E21[1, 0] = -3

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

E21 =
[   1    0    0 ]
[  -3    1    0 ]
[   0    0    1 ]
E21 A =
[   1    2    1 ]
[   0    2   -2 ]
[   0    4    1 ]

1절에서 손으로 만든 첫 단계 결과와 같다. 두 번째 단계도 마찬가지로 쓸 수 있다.

E32 = I.copy()
E32[2, 1] = -2

E = E32 @ E21
print(show_matrix(E, "E = E32 E21 ="))
print(show_matrix(E @ A, "E A ="))
print()
print("E A 가 U 와 같은가 :", np.allclose(E @ A, U))
print("E b 가 c 와 같은가 :", np.allclose(E @ b, c))

E = E32 E21 =
[   1    0    0 ]
[  -3    1    0 ]
[   6   -2    1 ]
E A =
[   1    2    1 ]
[   0    2   -2 ]
[   0    0    5 ]

E A 가 U 와 같은가 : True
E b 가 c 와 같은가 : True

소거의 전 과정이 행렬 하나에 담겼다.

$E$ 의 세 번째 행이 $(6, -2, 1)$ 인 것을 직접 확인해 보자. 서술 파트에서 본 것처럼 $EA$ 의 각 행은 $A$ 의 행들의 선형결합이고, 그 계수가 $E$ 의 해당 행이다.

print("E 의 3행   :", E[2])
print("직접 결합  :", E[2, 0] * A[0] + E[2, 1] * A[1] + E[2, 2] * A[2])
print("(E A) 의 3행:", (E @ A)[2])

E 의 3행   : [ 6. -2.  1.]
직접 결합  : [0. 0. 5.]
(E A) 의 3행: [0. 0. 5.]

왼쪽 곱과 오른쪽 곱은 다르다¶

$EA$ 는 행에 대한 연산이고 $AE$ 는 열에 대한 연산이다. 결과가 전혀 다르다.

print(show_matrix(E21 @ A, "E21 A   (2행에서 1행의 3배를 뺀다)"))
print(show_matrix(A @ E21, "A E21   (1열에서 2열의 3배를 뺀다)"))
print("같은가 :", np.allclose(E21 @ A, A @ E21))

E21 A   (2행에서 1행의 3배를 뺀다)
[   1    2    1 ]
[   0    2   -2 ]
[   0    4    1 ]
A E21   (1열에서 2열의 3배를 뺀다)
[   -5     2     1 ]
[  -21     8     1 ]
[  -12     4     1 ]
같은가 : False

치환행렬¶

행 교환도 같은 방식으로 쓸 수 있다. 단위행렬의 행 순서를 바꾸면 된다.

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

print(show_matrix(P @ A, "P A   (1행과 2행 교환)"))
print(show_matrix(A @ P, "A P   (1열과 2열 교환)"))

P A   (1행과 2행 교환)
[  3   8   1 ]
[  1   2   1 ]
[  0   4   1 ]
A P   (1열과 2열 교환)
[  2   1   1 ]
[  8   3   1 ]
[  4   0   1 ]

4. 소거해도 교점은 움직이지 않는다¶

각 단계의 증강행렬에서 행 하나가 평면 하나이다. 단계를 넘길 때마다 평면은 바뀌지만 세 평면이 만나는 점은 그대로여야 한다. 슬라이더로 단계를 넘기며 확인해 보자.

해가 모든 단계의 세 평면 위에 있으므로, 평면 조각을 해가 있는 자리를 중심으로 그린다.

solution = np.linalg.solve(A, b)
plane_colors = ["#17becf", "#9467bd", "#8c564b"]

frames_data, labels = [], []
for M_step, 설명 in history:
    traces = [plane(normal=M_step[i, :-1], through=solution, extent=4.5,
                    color=plane_colors[i], opacity=0.35, name=f"식 {i + 1}")
              for i in range(3)]
    traces.append(go.Scatter3d(x=[solution[0]], y=[solution[1]], z=[solution[2]],
                               mode="markers", marker=dict(size=8, color="red"),
                               name="해"))
    frames_data.append(traces)
    labels.append(설명)

slider_figure(frames_data, labels,
              layout3d("소거 단계별 세 평면", extent=5),
              prefix="", initial=0)

print("마지막 단계의 3행 :", history[-1][0][2])
print("즉  0x + 0y + 5z = -10   ->   z =", history[-1][0][2][-1] / history[-1][0][2][2])

마지막 단계의 3행 : [  0.   0.   5. -10.]
즉  0x + 0y + 5z = -10   ->   z = -2.0

5. 피벗이 0이 되는 경우¶

5-1. 행을 바꾸면 계속할 수 있는 경우¶

두 번째 행을 조금 바꿔서 두 번째 피벗이 0이 되도록 만들어 보자.

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

try:
    eliminate(A_swap, b)
except ValueError as err:
    print("소거 중단 :", err)

소거 중단 : 2번째 피벗이 0이다. 행 교환이 필요하다.

두 번째 행에서 첫 번째 행의 3배를 빼면 $(0, 0, -2)$ 가 되어 두 번째 피벗 자리가 0이 된다. 그런데 세 번째 행의 두 번째 성분은 4이므로 0이 아니다. 두 행을 바꾸면 소거를 계속할 수 있다.

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

U_swap, c_swap, _ = eliminate(P23 @ A_swap, P23 @ b)
print(show_matrix(U_swap, "U ="))
print("pivots =", np.diag(U_swap))
print("det(A_swap) =", np.linalg.det(A_swap))

U =
[   1    2    1 ]
[   0    4    1 ]
[   0    0   -2 ]
pivots = [ 1.  4. -2.]
det(A_swap) = 8.000000000000002

피벗이 세 개 나왔고 행렬식도 0이 아니다. 행을 바꾼 것뿐이므로 해도 그대로이다.

print("행 교환 후 푼 해 :", back_substitute(U_swap, c_swap))
print("np.linalg.solve  :", np.linalg.solve(A_swap, b))

행 교환 후 푼 해 : [ 2.5   1.25 -3.  ]
np.linalg.solve  : [ 2.5   1.25 -3.  ]

5-2. 소거가 끝까지 가지 못하는 경우¶

이번에는 세 번째 행을 첫 번째 행의 2배로 두어 보자.

A_sing = np.array([[1, 2, 1],
                   [3, 8, 1],
                   [2, 4, 2]], dtype=float)      # 3행 = 2 x 1행

U_sing, c_sing, _ = eliminate(A_sing, b)
print(show_matrix(U_sing, "U ="))
print("대각 성분 :", np.diag(U_sing))
print("0이 아닌 피벗 개수 :", int(np.sum(np.abs(np.diag(U_sing)) > 1e-12)))
print("det =", np.linalg.det(A_sing))

U =
[   1    2    1 ]
[   0    2   -2 ]
[   0    0    0 ]
대각 성분 : [1. 2. 0.]
0이 아닌 피벗 개수 : 2
det = 0.0

이번에는 소거가 중간에 멈추지 않았다. eliminate 는 마지막 열은 검사하지 않기 때문이다. 그러나 결과를 보면 마지막 대각 성분이 0이라 피벗은 두 개뿐이다. 행렬식도 0이다. 소거를 끝냈다고 끝난 것이 아니라, 대각 성분을 확인해야 한다.

우변까지 함께 보면 마지막 행이 무엇을 말하는지 알 수 있다.

print("마지막 행 (계수 | 우변) :", U_sing[2], "|", c_sing[2])
print()
print(f"즉  0 = {c_sing[2]:g}  라는 식이 되었다.")
print("좌변은 0인데 우변은 0이 아니므로 이 시스템은 해가 없다.")

마지막 행 (계수 | 우변) : [0. 0. 0.] | -2.0

즉  0 = -2  라는 식이 되었다.
좌변은 0인데 우변은 0이 아니므로 이 시스템은 해가 없다.

지난 강의의 언어로 말하면, $\vv{b}$ 가 $A$ 의 열들이 만드는 평면 밖에 있는 경우이다. 우변을 그 평면 안으로 옮기면 마지막 행이 $0 = 0$ 이 되고, 이때는 해가 무수히 많아진다.

b_ok = A_sing @ np.array([1.0, 1.0, 1.0])     # 열들의 결합이므로 반드시 평면 안에 있다
_, c_ok, _ = eliminate(A_sing, b_ok)

print("b_ok =", b_ok)
print("소거 후 마지막 우변 :", c_ok[2], " -> 0 = 0 이므로 모순이 없다")

b_ok = [ 4. 12.  8.]
소거 후 마지막 우변 : 0.0  -> 0 = 0 이므로 모순이 없다

마치며...¶

이번 실습에서 확인한 것을 정리하면 다음과 같다.

서술 파트의 내용	이 노트북의 코드
증강행렬 $[A \mid \vv{b}]$	`np.hstack([A, b.reshape(-1, 1)])`
소거	`eliminate(A, b)` 의 이중 `for` 문
곱수(multiplier)	`M[row, col] / pivot`
후진 대입	`back_substitute(U, c)` 의 역순 `for` 문
소거 한 단계 $= E_{ij}$	`E21 = np.eye(3); E21[1, 0] = -3`
소거 전체 $= E$	`E = E32 @ E21`, `E @ A == U`
소거는 해를 바꾸지 않는다	4절의 슬라이더
피벗이 모자라면 특이행렬	`np.diag(U)` 확인

더 해 볼 것¶

eliminate 가 마지막 대각 성분까지 검사하도록 고쳐 보자. 피벗의 개수를 함께 돌려주면 특이행렬을 바로 알아볼 수 있다.
A 의 값을 바꿔 가며 피벗을 관찰해 보자. 피벗을 모두 곱한 값과 np.linalg.det(A) 는 어떤 관계인가? (행 교환을 하지 않은 경우로 한정한다.)
4절의 슬라이더 그림에서 through=solution 을 through=None 으로 바꾸면 무엇이 달라지는가? 왜 해를 중심으로 평면을 그렸는지 생각해 보자.
$E_{21}$ 의 $(2,1)$ 성분을 -3 이 아니라 +3 으로 두면 $E_{21}A$ 는 무엇이 되는가? 그 행렬은 소거를 되돌리는 것과 어떤 관계인가?

마지막 문제가 다음 강의로 이어진다. 다음 강의에서는 행렬 곱셈을 네 가지 관점에서 정리하고, 역행렬(Inverse matrix)에 대해 알아보도록 하자.