Lecture 3. 행렬 곱셈과 역행렬 — 파이썬 실습 - 선형대수

L3 서술 파트에서 행렬 곱셈의 네 가지 관점과 역행렬을 다루었다. 이 노트북에서는 네 관점을 각각 함수로 구현해 결과가 정말 같은지 확인하고, 역행렬이 없다는 것이 무엇을 뜻하는지 그림으로 보도록 하자.

서술 파트의 내용	여기서 확인하는 방법
곱셈을 보는 네 가지 관점	네 개의 함수로 구현해 `A @ B` 와 대조한다
열 × 행은 랭크 1 행렬이다	외적을 만들어 랭크와 열의 방향을 확인한다
$(AB)^{-1} = B^{-1}A^{-1}$	직접 계산해 순서를 바꾸면 어떻게 되는지 본다
역행렬이 없다 $=$ 공간이 납작해진다	정육면체를 변환해 회전시켜 본다
가우스-조르당	$[A \mid I] \to [I \mid A^{-1}]$ 을 직접 구현한다

0. 준비¶

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, layout3d, spin_figure, show_matrix

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

numpy 2.5.2

A = np.array([[1, 2],
              [3, 4],
              [5, 6]], dtype=float)        # 3 x 2

B = np.array([[7,  8,  9],
              [10, 11, 12]], dtype=float)  # 2 x 3

print("A:", A.shape, "  B:", B.shape, "  AB:", (A @ B).shape)
print(show_matrix(A @ B, "A B ="))

A: (3, 2)   B: (2, 3)   AB: (3, 3)
A B =
[   27    30    33 ]
[   61    68    75 ]
[   95   106   117 ]

$3 \times 2$ 와 $2 \times 3$ 을 곱해 $3 \times 3$ 이 나왔다. 가운데 숫자 2가 서로 맞아야 하고, 결과의 크기는 바깥 숫자를 따라간다. 맞지 않으면 오류가 난다.

try:
    B @ B
except ValueError as err:
    print("곱할 수 없다 :", err)

곱할 수 없다 : matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

1. 네 가지 관점을 각각 구현하기¶

서술 파트에서 본 네 가지 방법을 그대로 함수로 옮긴다. numpy 의 @ 는 쓰지 않고, 각 관점이 말하는 대로만 계산한다.

def by_entry(A, B):
    """(1) 성분 : (i,j) 성분은 A 의 i행과 B 의 j열의 내적"""
    m, n = A.shape
    p = B.shape[1]
    C = np.zeros((m, p))
    for i in range(m):
        for j in range(p):
            C[i, j] = sum(A[i, k] * B[k, j] for k in range(n))
    return C

def by_columns(A, B):
    """(2) 열 : AB 의 각 열은 A 의 열들의 선형결합. 계수는 B 의 그 열."""
    cols = []
    for j in range(B.shape[1]):
        c = np.zeros(A.shape[0])
        for k in range(A.shape[1]):
            c = c + B[k, j] * A[:, k]        # A 의 k열을 B[k, j] 배 해서 더한다
        cols.append(c)
    return np.column_stack(cols)

def by_rows(A, B):
    """(3) 행 : AB 의 각 행은 B 의 행들의 선형결합. 계수는 A 의 그 행."""
    rows = []
    for i in range(A.shape[0]):
        r = np.zeros(B.shape[1])
        for k in range(B.shape[0]):
            r = r + A[i, k] * B[k, :]        # B 의 k행을 A[i, k] 배 해서 더한다
        rows.append(r)
    return np.vstack(rows)

def by_outer(A, B):
    """(4) 열 x 행 : AB 는 (A의 k열)(B의 k행) 들의 합"""
    total = np.zeros((A.shape[0], B.shape[1]))
    for k in range(A.shape[1]):
        total = total + np.outer(A[:, k], B[k, :])
    return total

답 = A @ B
for 이름, 함수 in [("by_entry", by_entry), ("by_columns", by_columns),
                  ("by_rows", by_rows), ("by_outer", by_outer)]:
    print(f"{이름:<11} 일치 : {np.allclose(함수(A, B), 답)}")

by_entry    일치 : True
by_columns  일치 : True
by_rows     일치 : True
by_outer    일치 : True

네 함수가 모두 같은 답을 준다. 계산 과정을 눈으로 확인해 보자.

print("[열 관점] AB 의 1열 = 7(A의 1열) + 10(A의 2열)")
print("  ", 7 * A[:, 0] + 10 * A[:, 1], " <-> ", 답[:, 0])
print()
print("[행 관점] AB 의 1행 = 1(B의 1행) + 2(B의 2행)")
print("  ", 1 * B[0, :] + 2 * B[1, :], " <-> ", 답[0, :])

[열 관점] AB 의 1열 = 7(A의 1열) + 10(A의 2열)
   [27. 61. 95.]  <->  [27. 61. 95.]

[행 관점] AB 의 1행 = 1(B의 1행) + 2(B의 2행)
   [27. 30. 33.]  <->  [27. 30. 33.]

2. 열 × 행은 랭크 1 행렬이다¶

열 하나( $3 \times 1$ )와 행 하나( $1 \times 3$ )를 곱하면 숫자가 아니라 $3 \times 3$ 행렬이 나온다. np.outer 가 이 계산을 해 준다.

조각1 = np.outer(A[:, 0], B[0, :])
조각2 = np.outer(A[:, 1], B[1, :])

print(show_matrix(조각1, "(A의 1열)(B의 1행) ="))
print(show_matrix(조각2, "(A의 2열)(B의 2행) ="))
print(show_matrix(조각1 + 조각2, "두 조각의 합 ="))
print("A B 와 같은가 :", np.allclose(조각1 + 조각2, 답))

(A의 1열)(B의 1행) =
[   7    8    9 ]
[  21   24   27 ]
[  35   40   45 ]
(A의 2열)(B의 2행) =
[  20   22   24 ]
[  40   44   48 ]
[  60   66   72 ]
두 조각의 합 =
[   27    30    33 ]
[   61    68    75 ]
[   95   106   117 ]
A B 와 같은가 : True

조각 하나하나가 어떤 행렬인지 보자. 랭크를 확인하고, 열끼리 비례하는지 확인한다.

print("rank(조각1) =", np.linalg.matrix_rank(조각1))
print("rank(A B)   =", np.linalg.matrix_rank(답))
print()
print("조각1 의 열들 :")
for j in range(3):
    print(f"  열{j + 1} = {조각1[:, j]}   =  {B[0, j]:g} x (A의 1열) {A[:, 0]}")

rank(조각1) = 1
rank(A B)   = 2

조각1 의 열들 :
  열1 = [ 7. 21. 35.]   =  7 x (A의 1열) [1. 3. 5.]
  열2 = [ 8. 24. 40.]   =  8 x (A의 1열) [1. 3. 5.]
  열3 = [ 9. 27. 45.]   =  9 x (A의 1열) [1. 3. 5.]

세 열이 모두 $A$ 의 첫 번째 열의 배수이다. 그래서 랭크가 1이고, 세 열이 한 직선 위에 놓인다. 그리고 랭크 1 조각 두 개를 더하니 랭크 2가 되었다.

이 관점은 나중에 행렬을 분해할 때 쓰인다. 지금은 곱을 이렇게도 쪼갤 수 있다는 것만 확인하고 넘어간다.

3. $AB$ 와 $BA$ 는 다르다¶

행렬 곱셈에는 교환법칙이 성립하지 않는다. 크기가 맞아서 둘 다 계산되는 경우에도 결과가 다르다.

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

print(show_matrix(P @ Q, "P Q ="))
print(show_matrix(Q @ P, "Q P ="))
print("같은가 :", np.allclose(P @ Q, Q @ P))

P Q =
[  2   1 ]
[  1   1 ]
Q P =
[  1   1 ]
[  1   2 ]
같은가 : False

크기부터 달라지는 경우도 있다. 우리 예제의 $A$ 와 $B$ 가 그렇다.

print("A B :", (A @ B).shape)
print("B A :", (B @ A).shape)

A B : (3, 3)
B A : (2, 2)

블록 곱셈¶

행렬을 블록으로 나누어 곱해도 결과가 같은지 확인해 보자. $B$ 를 좌우 두 블록으로 나눈다.

B1, B2 = B[:, :1], B[:, 1:]          # 1열 / 나머지 2열

블록결과 = np.hstack([A @ B1, A @ B2])
print(show_matrix(블록결과, "[A B1 | A B2] ="))
print("A B 와 같은가 :", np.allclose(블록결과, 답))

[A B1 | A B2] =
[   27    30    33 ]
[   61    68    75 ]
[   95   106   117 ]
A B 와 같은가 : True

4. 역행렬 (Inverse matrix)¶

서술 파트에서 쓴 $2 \times 2$ 예제이다.

M = np.array([[1, 3],
              [2, 7]], dtype=float)

M_inv = np.linalg.inv(M)
print(show_matrix(M_inv, "M^-1 ="))
print(show_matrix(M @ M_inv, "M M^-1 ="))
print(show_matrix(M_inv @ M, "M^-1 M ="))

M^-1 =
[   7   -3 ]
[  -2    1 ]
M M^-1 =
[  1   0 ]
[  0   1 ]
M^-1 M =
[  1   0 ]
[  0   1 ]

양쪽 모두 단위행렬이 된다. 정방행렬에서는 한쪽만 확인해도 나머지가 따라온다.

곱의 역행렬은 순서가 뒤집힌다¶

N = np.array([[2, 1],
              [1, 1]], dtype=float)

왼쪽 = np.linalg.inv(M @ N)
바른순서 = np.linalg.inv(N) @ np.linalg.inv(M)      # B^-1 A^-1
틀린순서 = np.linalg.inv(M) @ np.linalg.inv(N)      # A^-1 B^-1

print(show_matrix(왼쪽, "(M N)^-1 ="))
print(show_matrix(바른순서, "N^-1 M^-1 ="))
print(show_matrix(틀린순서, "M^-1 N^-1 ="))
print()
print("(M N)^-1 == N^-1 M^-1 :", np.allclose(왼쪽, 바른순서))
print("(M N)^-1 == M^-1 N^-1 :", np.allclose(왼쪽, 틀린순서))

(M N)^-1 =
[    9    -4 ]
[  -11     5 ]
N^-1 M^-1 =
[    9    -4 ]
[  -11     5 ]
M^-1 N^-1 =
[   10   -13 ]
[   -3     4 ]

(M N)^-1 == N^-1 M^-1 : True
(M N)^-1 == M^-1 N^-1 : False

역행렬로 방정식을 풀지 않는 이유¶

서술 파트에서 $\vv{x} = A^{-1}\vv{b}$ 는 해가 무엇인지 말해 주는 식이지 해를 구하는 방법은 아니라고 하였다. 실제로 얼마나 차이가 나는지 확인해 보자.

힐베르트 행렬은 계산이 까다로운 것으로 잘 알려진 행렬이다. 정답이 전부 1인 문제를 만들어 놓고 두 방법으로 풀어 본다.

from scipy.linalg import hilbert

n = 10
H = hilbert(n)
x_true = np.ones(n)                   # 정답을 정해 놓고
b_h = H @ x_true                      # 그에 맞는 우변을 만든다

x_inv = np.linalg.inv(H) @ b_h        # 역행렬을 구해서 곱하기
x_solve = np.linalg.solve(H, b_h)     # 소거로 바로 풀기

오차_inv = np.max(np.abs(x_inv - x_true))
오차_solve = np.max(np.abs(x_solve - x_true))

print(f"조건수          : {np.linalg.cond(H):.1e}")
print(f"inv @ b 의 오차 : {오차_inv:.3e}")
print(f"solve   의 오차 : {오차_solve:.3e}")
print(f"inv 쪽이 {오차_inv / 오차_solve:.0f} 배 더 부정확하다.")

조건수          : 1.6e+13
inv @ b 의 오차 : 5.547e-03
solve   의 오차 : 1.692e-04
inv 쪽이 33 배 더 부정확하다.

같은 문제를 푸는데 오차가 다르다. 역행렬을 따로 구한 다음 곱하면 반올림 오차가 두 번 쌓이기 때문이다. 조건수가 큰 행렬일수록 차이가 벌어진다. np.linalg.solve 를 쓰는 것이 맞다.

조건수(condition number)는 입력의 작은 오차가 답에서 얼마나 커지는지를 나타내는 값이다. 여기서는 그런 것이 있다는 정도만 알아 두면 된다.

5. 역행렬이 없는 경우¶

$A\vv{x} = \vv{0}$ 을 만족하는 $\vv{x} \neq \vv{0}$ 이 있으면 역행렬은 존재하지 않는다. 서술 파트의 예제로 확인해 보자.

S = np.array([[1, 3],
              [2, 6]], dtype=float)     # 2열 = 3 x (1열)

x0 = np.array([3.0, -1.0])
print("S x0 =", S @ x0, " (영벡터인데 x0 는 영벡터가 아니다)")
print("det(S) =", np.linalg.det(S))

try:
    np.linalg.inv(S)
except np.linalg.LinAlgError as err:
    print("역행렬 없음 :", err)

S x0 = [0. 0.]  (영벡터인데 x0 는 영벡터가 아니다)
det(S) = 0.0
역행렬 없음 : Singular matrix

정육면체가 어디로 가는지 보기¶

행렬을 변환으로 보면 역행렬이 없다는 것이 무엇인지 눈으로 볼 수 있다. 단위 정육면체의 여덟 꼭짓점에 행렬을 곱해서 어떤 모양이 되는지 그려 보자.

CUBE = np.array([[i, j, k] for i in (0, 1) for j in (0, 1) for k in (0, 1)], dtype=float)
EDGES = [(a, b) for a in range(8) for b in range(a + 1, 8)
         if np.sum(np.abs(CUBE[a] - CUBE[b])) == 1]      # 길이 1인 모서리만


def cube_trace(M, color, name, width=5):
    """정육면체를 M 으로 변환한 결과를 선분 하나로 묶어 만든다."""
    P = CUBE @ M.T                     # 꼭짓점 하나하나에 M 을 곱한 것과 같다
    xs, ys, zs = [], [], []
    for a, b in EDGES:
        xs += [P[a, 0], P[b, 0], None]   # None 을 끼우면 선이 끊긴다
        ys += [P[a, 1], P[b, 1], None]
        zs += [P[a, 2], P[b, 2], None]
    return go.Scatter3d(x=xs, y=ys, z=zs, mode="lines",
                        line=dict(color=color, width=width), name=name)

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

print("det =", np.linalg.det(A_good))

spin_figure([cube_trace(np.eye(3), "#bbbbbb", "원래 정육면체", width=3),
             cube_trace(A_good, COLORS["output"], "변환 결과")],
            title="가역인 경우 : 부피가 있다", extent=3)

det = 2.0

이번에는 세 번째 열을 앞의 두 열의 합으로 바꾼다. 열 하나만 바뀌었다.

A_flat = A_good.copy()
A_flat[:, 2] = A_good[:, 0] + A_good[:, 1]      # col3 = col1 + col2

print(show_matrix(A_flat, "A_flat ="))
print("det =", np.linalg.det(A_flat))
print("rank =", np.linalg.matrix_rank(A_flat))

spin_figure([cube_trace(np.eye(3), "#bbbbbb", "원래 정육면체", width=3),
             cube_trace(A_flat, COLORS["output"], "변환 결과")],
            title="특이인 경우 : 납작해진다", extent=3)

A_flat =
[  1   1   2 ]
[  0   1   1 ]
[  1   0   1 ]
det = 0.0
rank = 2

moved = CUBE @ A_flat.T
겹친쌍 = [(a, b) for a in range(8) for b in range(a + 1, 8)
          if np.allclose(moved[a], moved[b])]

print("변환 후 같은 자리로 간 꼭짓점 쌍 :")
for a, b in 겹친쌍:
    print(f"  {CUBE[a]} 와 {CUBE[b]}  ->  {moved[a]}")

a, b = 겹친쌍[0]
차이 = CUBE[a] - CUBE[b]
print()
print("두 점의 차이 :", 차이)
print("A_flat 을 곱하면 :", A_flat @ 차이, " (영벡터)")

변환 후 같은 자리로 간 꼭짓점 쌍 :
  [0. 0. 1.] 와 [1. 1. 0.]  ->  [2. 1. 1.]

두 점의 차이 : [-1. -1.  1.]
A_flat 을 곱하면 : [0. 0. 0.]  (영벡터)

서로 다른 두 점이 같은 곳으로 갔다. 도착점만 보고 어느 점에서 왔는지 알 수 없으므로 되돌리는 변환을 만들 수 없다. 이것이 역행렬이 없다는 것의 의미이다.

6. 가우스-조르당 직접 구현하기¶

$[A \mid I]$ 를 놓고 소거해서 왼쪽을 $I$ 로 만들면 오른쪽에 $A^{-1}$ 이 남는다. 지난 강의의 소거와 다른 점은 두 가지이다. 피벗을 1로 맞추고, 위쪽도 함께 소거한다.

def gauss_jordan_inverse(A):
    """[A | I] 를 소거해 [I | A^-1] 을 만든다. 필요하면 행을 바꾼다."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    M = np.hstack([A.copy(), np.eye(n)])

    for col in range(n):
        if abs(M[col, col]) < 1e-12:                      # 피벗이 0이면
            아래 = np.flatnonzero(np.abs(M[col + 1:, col]) > 1e-12)
            if len(아래) == 0:
                raise ValueError("역행렬이 존재하지 않는다.")
            바꿀행 = col + 1 + 아래[0]
            M[[col, 바꿀행]] = M[[바꿀행, col]]            # 행 교환

        M[col] = M[col] / M[col, col]                     # 피벗을 1로 맞춘다
        for row in range(n):
            if row != col:
                M[row] = M[row] - M[row, col] * M[col]    # 위아래 모두 소거

    return M[:, n:]

print(show_matrix(gauss_jordan_inverse(M), "gauss_jordan_inverse(M) ="))
print(show_matrix(np.linalg.inv(M), "np.linalg.inv(M) ="))

gauss_jordan_inverse(M) =
[   7   -3 ]
[  -2    1 ]
np.linalg.inv(M) =
[   7   -3 ]
[  -2    1 ]

서술 파트에서 손으로 구한 것과 같다. 3차 행렬로도 확인해 보자.

mine = gauss_jordan_inverse(A_good)
print(show_matrix(mine, "gauss_jordan_inverse(A_good) ="))
print("np.linalg.inv 와 같은가 :", np.allclose(mine, np.linalg.inv(A_good)))
print("A_good @ mine 이 I 인가 :", np.allclose(A_good @ mine, np.eye(3)))

gauss_jordan_inverse(A_good) =
[   0.5   -0.5    0.5 ]
[   0.5    0.5   -0.5 ]
[  -0.5    0.5    0.5 ]
np.linalg.inv 와 같은가 : True
A_good @ mine 이 I 인가 : True

try:
    gauss_jordan_inverse(A_flat)
except ValueError as err:
    print("특이행렬 :", err)

특이행렬 : 역행렬이 존재하지 않는다.

오른쪽에 왜 $A^{-1}$ 이 남는가¶

지난 강의에서 소거 전체를 행렬 $E$ 하나로 쓸 수 있었다. 증강행렬에 $E$ 를 곱하는 것은 블록마다 곱하는 것과 같으므로 $E[A \mid I] = [EA \mid E]$ 이다. 왼쪽이 $I$ 가 되도록 소거했으니 $EA = I$ 이고, 따라서 $E = A^{-1}$ 이다.

직접 확인해 보자. 소거 전체에 해당하는 행렬은 결국 우리가 얻은 그 역행렬이어야 한다.

E = gauss_jordan_inverse(A_good)      # 소거 전 과정에 해당하는 행렬
print("E A 가 I 인가 :", np.allclose(E @ A_good, np.eye(3)))
print("즉 E 가 A^-1 인가 :", np.allclose(E, np.linalg.inv(A_good)))

E A 가 I 인가 : True
즉 E 가 A^-1 인가 : True

마치며...¶

서술 파트의 내용	이 노트북의 코드
성분 관점	`by_entry` 의 이중 `for` 문
열 관점	`by_columns` — `B[k, j] * A[:, k]` 의 합
행 관점	`by_rows` — `A[i, k] * B[k, :]` 의 합
열 × 행 관점	`by_outer` — `np.outer(A[:, k], B[k, :])` 의 합
랭크 1	`np.linalg.matrix_rank(조각1)`
$(AB)^{-1} = B^{-1}A^{-1}$	4절의 `바른순서` / `틀린순서`
역행렬을 구해 풀지 않는다	힐베르트 행렬의 오차 비교
역행렬이 없다 $=$ 납작해진다	5절의 정육면체와 겹친 꼭짓점
가우스-조르당	`gauss_jordan_inverse`

더 해 볼 것¶

by_outer 에서 조각을 하나만 더해 보고, 두 개를 다 더했을 때와 랭크가 어떻게 달라지는지 확인해 보자. 조각을 앞에서부터 $k$ 개만 더하면 랭크가 몇이 되는가?
A_flat 의 열 대신 행 하나를 다른 두 행의 합으로 바꿔 보자. 정육면체는 어떻게 되는가? 납작해지는가?
gauss_jordan_inverse 가 피벗을 1로 맞추는 줄을 지우면 무엇이 잘못되는가? 직접 지워 보고 결과를 확인해 보자.
$2 \times 2$ 행렬 $\begin{bmatrix} a & b \\ c & d \end{bmatrix}$ 에 대해 gauss_jordan_inverse 를 손으로 따라가면 잘 알려진 공식이 나온다. $ad - bc$ 가 어디에서 나오는지 찾아보자.

다음 강의에서는 소거 행렬과 역행렬을 합쳐서 이 책의 첫 번째 분해인 $A = LU$ 를 만든다.