Lecture 20. 크래머 공식, 역행렬, 부피 — 파이썬 실습 - 선형대수

L20 서술 파트의 심장은 행렬식이 부피라는 것이었다. 가는 길에 여인수로 역행렬과 크래머 공식도 얻었다.

이 노트북에서는 딸림행렬로 역행렬을 직접 만들어 np.linalg.inv 와 대조하고, 크래머 공식이 소거보다 얼마나 느리고 얼마나 덜 정확한지 재고, 단위 정육면체가 납작해지는 과정을 슬라이더로 눌러 본다.

서술 파트의 내용	여기서 확인하는 방법
$AC^{\mathsf{T}} = (\det A)I$	여인수 행렬을 만들어 곱한다
비대각이 0인 이유	행을 갈아끼운 행렬의 행렬식을 직접 계산
$A^{-1} = C^{\mathsf{T}}/\det A$	`np.linalg.inv` 와 대조
크래머 공식	값은 맞지만 느리고 덜 정확하다
$\lvert\det\rvert$ 는 부피	슬라이더로 눌러 0 으로 보낸다
$\det = 0$ 의 네 얼굴	부피·종속·영공간·역행렬을 한 표에
야코비안	극좌표의 $r$ 을 수치미분으로

0. 준비¶

import time

import numpy as np
import plotly.graph_objects as go
from scipy.linalg import null_space

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

np.set_printoptions(precision=3, suppress=True)
rng = np.random.default_rng(20)
print("numpy", np.__version__)

numpy 2.5.2

1. $AC^{\mathsf{T}} = (\det A)\,I$ ¶

여인수를 제자리에 늘어놓은 행렬을 만들고, $A$ 와 그 전치를 곱해 본다.

def 여인수행렬(A):
    """C[i][j] = (-1)^(i+j) 곱하기 i행 j열을 지운 소행렬식."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    return np.array([[(-1.0) ** (i + j)
                      * np.linalg.det(np.delete(np.delete(A, i, 0), j, 1))
                      for j in range(n)] for i in range(n)])

A = np.array([[1.0, 2.0, 1.0],
              [3.0, 8.0, 1.0],
              [0.0, 4.0, 1.0]])            # L18, L19 에서 쓰던 행렬
행렬식 = np.linalg.det(A)

C = 여인수행렬(A)
print(show_matrix(A, "A"))
print(show_matrix(np.round(C), "C  (여인수 행렬)"))
print(show_matrix(np.round(C.T), "C transpose  (딸림행렬)"))
print()
print(show_matrix(np.round(A @ C.T), "A C^T"))
print("det A =", round(행렬식))
print("(det A) I 와 같은가 :", np.allclose(A @ C.T, 행렬식 * np.eye(3)))

A
[  1   2   1 ]
[  3   8   1 ]
[  0   4   1 ]
C  (여인수 행렬)
[   4   -3   12 ]
[   2    1   -4 ]
[  -6    2    2 ]
C transpose  (딸림행렬)
[   4    2   -6 ]
[  -3    1    2 ]
[  12   -4    2 ]

A C^T
[  10    0    0 ]
[  -0   10    0 ]
[  -0    0   10 ]
det A = 10
(det A) I 와 같은가 : True

비대각 성분이 0인 이유를 직접 확인하자. 1행과 2열의 곱은 2행을 1행으로 갈아끼운 행렬의 행렬식과 같아야 한다. 그 행렬은 같은 행이 둘이라 행렬식이 0이다.

i, k = 0, 1                                   # 1행의 성분에 2행의 여인수를 곱한다
합 = sum(A[i, j] * C[k, j] for j in range(3))

갈아낀 = A.copy()
갈아낀[k] = A[i]                               # 2행을 1행으로 바꾼다

print("sum_j a[1][j] C[2][j] =", 합)
print()
print(show_matrix(갈아낀, "2행을 1행으로 갈아끼운 행렬"))
print("그 행렬식 :", np.linalg.det(갈아낀))
print("같은가 :", np.isclose(합, np.linalg.det(갈아낀)))
print("두 행이 같은가 :", np.allclose(갈아낀[0], 갈아낀[1]))

sum_j a[1][j] C[2][j] = 0.0

2행을 1행으로 갈아끼운 행렬
[  1   2   1 ]
[  1   2   1 ]
[  0   4   1 ]
그 행렬식 : 0.0
같은가 : True
두 행이 같은가 : True

2. 딸림행렬로 만든 역행렬¶

def 딸림역행렬(A):
    """A^-1 = C^T / det A. 공식 그대로 만든다."""
    A = np.asarray(A, dtype=float)
    행렬식 = np.linalg.det(A)
    if abs(행렬식) < 1e-12:
        raise ValueError("det 가 0 이라 역행렬이 없다")
    return 여인수행렬(A).T / 행렬식

print(show_matrix(딸림역행렬(A), "C^T / det A"))
print(show_matrix(np.linalg.inv(A), "np.linalg.inv(A)"))
print("같은가 :", np.allclose(딸림역행렬(A), np.linalg.inv(A)))
print()
print(show_matrix(np.round(A @ 딸림역행렬(A), 10) + 0.0, "A 곱하기 그 역행렬"))

C^T / det A
[   0.4    0.2   -0.6 ]
[  -0.3    0.1    0.2 ]
[   1.2   -0.4    0.2 ]
np.linalg.inv(A)
[   0.4    0.2   -0.6 ]
[  -0.3    0.1    0.2 ]
[   1.2   -0.4    0.2 ]
같은가 : True

A 곱하기 그 역행렬
[  1   0   0 ]
[  0   1   0 ]
[  0   0   1 ]

$2 \times 2$ 에서는 L3에서 손으로 유도한 그 공식이 그대로 나온다.

A2 = np.array([[1.0, 2.0], [3.0, 4.0]])
a, b_, c, d = A2.ravel()
손공식 = np.array([[d, -b_], [-c, a]]) / (a * d - b_ * c)

print(show_matrix(딸림역행렬(A2), "딸림행렬로"))
print(show_matrix(손공식, "[d -b ; -c a] / (ad - bc)"))
print("같은가 :", np.allclose(딸림역행렬(A2), 손공식))

딸림행렬로
[    -2      1 ]
[   1.5   -0.5 ]
[d -b ; -c a] / (ad - bc)
[    -2      1 ]
[   1.5   -0.5 ]
같은가 : True

3. 크래머 공식¶

값은 맞다. 문제는 $n+1$ 개의 행렬식이 필요하다는 것이다.

def 크래머(A, b):
    """x_j = det(j열을 b 로 갈아끼운 행렬) / det A."""
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float)
    분모 = np.linalg.det(A)
    해 = np.empty(len(b))
    for j in range(len(b)):
        B = A.copy()
        B[:, j] = b                            # j 열만 b 로
        해[j] = np.linalg.det(B) / 분모
    return 해

b = np.array([2.0, 12.0, 2.0])                # L2 에서 쓰던 우변

print("크래머      :", 크래머(A, b))
print("solve       :", np.linalg.solve(A, b))
print("L2 의 답    : [ 2.  1. -2.]")
print()
for j in range(3):
    B = A.copy(); B[:, j] = b
    print(f"det B{j + 1} = {np.linalg.det(B):6.1f},"
          f"   x{j + 1} = {np.linalg.det(B) / 행렬식:+.1f}")

크래머      : [ 2.  1. -2.]
solve       : [ 2.  1. -2.]
L2 의 답    : [ 2.  1. -2.]

det B1 =   20.0,   x1 = +2.0
det B2 =   10.0,   x2 = +1.0
det B3 =  -20.0,   x3 = -2.0

이제 값이 맞는데도 왜 쓰지 않는지 재 본다. 크래머는 행렬식을 $n+1$ 번 구하고, 행렬식 하나가 이미 소거 한 번만큼 든다. 그래서 대략 $n$ 배 느리다.

def 재기(함수, 반복=20):
    """한 번 예열한 뒤 여러 번 재서 평균을 낸다."""
    함수()
    시작 = time.perf_counter()
    for _ in range(반복):
        결과 = 함수()
    return (time.perf_counter() - 시작) / 반복, 결과

print(f"{'n':>5}{'크래머 ms':>13}{'solve ms':>12}{'느린 배수':>11}"
      f"{'크래머의 오차':>15}{'solve 의 오차':>15}")
for n in (5, 10, 20, 50):
    M = rng.normal(size=(n, n))
    참해 = np.ones(n)
    우변 = M @ 참해

    t1, x1 = 재기(lambda: 크래머(M, 우변))
    t2, x2 = 재기(lambda: np.linalg.solve(M, 우변))
    print(f"{n:>5}{t1 * 1000:>13.3f}{t2 * 1000:>12.3f}{t1 / t2:>11.0f}"
          f"{np.abs(x1 - 참해).max():>15.1e}{np.abs(x2 - 참해).max():>15.1e}")

    n       크래머 ms    solve ms      느린 배수        크래머의 오차     solve 의 오차
    5        0.031       0.007          4        3.3e-16        2.2e-16
   10        0.061       0.008          7        1.8e-15        4.4e-16
   20        0.170       0.011         16        1.4e-14        4.0e-15
   50        1.786       0.024         74        2.9e-14        2.6e-14

“느린 배수” 가 $n$ 을 그대로 따라간다. 행렬식을 $n+1$ 번 구하니 당연한 결과이다.

서술 파트에서 말한 대로 크래머 공식의 쓸모는 계산이 아니라 민감도를 보는 데 있다. 분모가 $\det A$ 이므로, 행렬식이 작으면 $\vv{b}$ 의 작은 흔들림이 크게 증폭된다.

print("분모가 작으면 b 의 작은 흔들림이 크게 증폭된다.")
print(f"{'det A':>12}{'b 를 1e-6 흔들 때 x 의 변화':>28}")
for 눈금 in (1.0, 1e-2, 1e-4):
    M = np.array([[1.0, 1.0], [1.0, 1.0 + 눈금]])       # det = 눈금
    우변 = np.array([1.0, 1.0])
    흔든 = 우변 + np.array([1e-6, 0.0])
    변화 = np.abs(크래머(M, 흔든) - 크래머(M, 우변)).max()
    print(f"{np.linalg.det(M):>12.0e}{변화:>28.3e}")

분모가 작으면 b 의 작은 흔들림이 크게 증폭된다.
       det A        b 를 1e-6 흔들 때 x 의 변화
       1e+00                   2.000e-06
       1e-02                   1.010e-04
       1e-04                   1.000e-02

4. 부피 — 이 강의의 심장¶

단위 정육면체가 어디로 가는지 보고, 세 번째 열을 앞의 두 열의 결합으로 밀어 부피를 0으로 보내 보자.

def 평행육면체(M, 이름=""):
    """단위 정육면체를 M 으로 옮긴 도형과 세 열을 그린다."""
    M = np.asarray(M, dtype=float)
    꼭짓 = np.array([[a, b_, c] for a in (0, 1) for b_ in (0, 1)
                    for c in (0, 1)], dtype=float) @ M.T
    자료 = [go.Mesh3d(x=꼭짓[:, 0], y=꼭짓[:, 1], z=꼭짓[:, 2], alphahull=0,
                     color=COLORS["colspace"], opacity=0.32, name=이름 or "도형")]
    for c, 색, 라벨 in zip(range(3), (COLORS["input"], COLORS["second"],
                                    COLORS["third"]), ("열 1", "열 2", "열 3")):
        자료 += arrow([0, 0, 0], M[:, c], 색, 라벨)
    return 자료

M = np.array([[1.0, 1.0, 0.0],
              [0.0, 1.0, 1.0],
              [1.0, 0.0, 1.0]])

print(show_matrix(M, "M"))
print("det M =", np.linalg.det(M))
print("부피   =", abs(np.linalg.det(M)))

go.Figure(data=평행육면체(M),
          layout=layout3d(f"단위 정육면체의 상 : 부피 = {abs(np.linalg.det(M)):.0f}",
                          extent=2.4))

M
[  1   1   0 ]
[  0   1   1 ]
[  1   0   1 ]
det M = 2.0
부피   = 2.0

세 번째 열을 $(0, 1, 1-2t)$ 로 두고 $t$ 를 0에서 1로 밀어 보자. $\det = 2 - 2t$ 이므로 $t = 1$ 에서 도형이 평면 안으로 눌린다.

프레임, 라벨 = [], []
for t in np.linspace(0.0, 1.0, 11):
    Mt = M.copy()
    Mt[:, 2] = [0.0, 1.0, 1.0 - 2.0 * t]
    프레임.append(평행육면체(Mt))
    라벨.append(f"t={t:.1f}   det={np.linalg.det(Mt):+.2f}")

slider_figure(프레임, 라벨,
              layout3d("세 번째 열을 밀어 부피를 0 으로", extent=2.4),
              initial=0)

print(f"{'t':>6}{'det':>9}{'부피':>9}{'rank':>7}{'dim N(M)':>10}{'역행렬':>9}")
for t in (0.0, 0.5, 0.9, 1.0):
    Mt = M.copy()
    Mt[:, 2] = [0.0, 1.0, 1.0 - 2.0 * t]
    d = np.linalg.det(Mt)
    있음 = "있다" if abs(d) > 1e-12 else "없다"
    print(f"{t:>6.1f}{d:>9.2f}{abs(d):>9.2f}"
          f"{np.linalg.matrix_rank(Mt):>7}{null_space(Mt).shape[1]:>10}{있음:>9}")

     t      det       부피   rank  dim N(M)      역행렬
   0.0     2.00     2.00      3         0       있다
   0.5     1.00     1.00      3         0       있다
   0.9     0.20     0.20      3         0       있다
   1.0     0.00     0.00      2         1       없다

$t = 1$ 에서 네 가지가 한꺼번에 바뀐다. 부피가 0이 되고, 랭크가 줄고, 영공간이 커지고, 역행렬이 사라진다. 서술 파트의 사슬이 그것이다.

납작 = M.copy()
납작[:, 2] = [0.0, 1.0, -1.0]

print("세 번째 열 :", 납작[:, 2])
print("2 열 - 1 열 :", 납작[:, 1] - 납작[:, 0], "  -> 같다, 곧 종속")
print()
print("영공간의 기저 :", null_space(납작).ravel())
print("M 곱하기 그 벡터 :", (납작 @ null_space(납작)).ravel())

세 번째 열 : [ 0.  1. -1.]
2 열 - 1 열 : [ 0.  1. -1.]   -> 같다, 곧 종속

영공간의 기저 : [ 0.577 -0.577  0.577]
M 곱하기 그 벡터 : [-0. -0. -0.]

5. 삼각형과 부호¶

v1 = np.array([4.0, 1.0])
v2 = np.array([1.0, 3.0])
d2 = np.linalg.det(np.column_stack([v1, v2]))

print("평행사변형의 넓이 :", abs(d2))
print("삼각형의 넓이     :", abs(d2) / 2)
print()
# 신발끈 공식과 대조
점들 = np.array([[0.0, 0.0], v1, v2])
x, y = 점들[:, 0], 점들[:, 1]
신발끈 = 0.5 * abs(x[0] * (y[1] - y[2]) + x[1] * (y[2] - y[0])
                 + x[2] * (y[0] - y[1]))
print("신발끈 공식으로도 :", 신발끈)
print("같은가 :", np.isclose(abs(d2) / 2, 신발끈))

평행사변형의 넓이 : 11.000000000000002
삼각형의 넓이     : 5.500000000000001

신발끈 공식으로도 : 5.5
같은가 : True

직교행렬의 행렬식은 $\pm1$ 뿐이다. 회전이면 +1, 반사가 섞이면 -1 이다.

각도 = np.pi / 5
회전 = np.array([[np.cos(각도), -np.sin(각도)],
               [np.sin(각도),  np.cos(각도)]])
반사 = np.array([[1.0, 0.0], [0.0, -1.0]])

for 이름, Q in (("회전", 회전), ("반사", 반사), ("회전 후 반사", 반사 @ 회전)):
    print(f"{이름:>12} : Q^T Q = I 인가 {np.allclose(Q.T @ Q, np.eye(2))},"
          f"   det = {np.linalg.det(Q):+.0f}")

          회전 : Q^T Q = I 인가 True,   det = +1
          반사 : Q^T Q = I 인가 True,   det = -1
     회전 후 반사 : Q^T Q = I 인가 True,   det = -1

6. 야코비안 — 극좌표의 $r$ ¶

$x = r\cos\theta$ , $y = r\sin\theta$ 의 편미분을 모아 놓은 행렬의 행렬식이 $r$ 이다. 수치미분으로 확인하자.

def 야코비행렬(변환, 점, h=1e-6):
    """변환의 편미분을 수치로 구해 행렬로 모은다."""
    점 = np.asarray(점, dtype=float)
    n = len(점)
    기둥 = []
    for k in range(n):
        걸음 = np.zeros(n); 걸음[k] = h
        기둥.append((변환(점 + 걸음) - 변환(점 - 걸음)) / (2 * h))
    return np.column_stack(기둥)

극좌표 = lambda p: np.array([p[0] * np.cos(p[1]), p[0] * np.sin(p[1])])

print(f"{'r':>6}{'theta':>8}{'수치 det':>12}{'공식 r':>10}")
for r, 각 in ((1.0, 0.0), (2.0, 0.7), (3.0, 2.1), (0.5, 4.0)):
    J = 야코비행렬(극좌표, [r, 각])
    print(f"{r:>6.1f}{각:>8.2f}{np.linalg.det(J):>12.6f}{r:>10.1f}")

     r   theta      수치 det      공식 r
   1.0    0.00    1.000000       1.0
   2.0    0.70    2.000000       2.0
   3.0    2.10    3.000000       3.0
   0.5    4.00    0.500000       0.5

곱해 주던 $r$ 은 아주 작은 사각형이 좌표 변환으로 얼마나 늘어나는지의 부피 배율이었다. 원점에서 멀어질수록 같은 $d\theta$ 가 더 긴 호를 그으니 배율이 $r$ 에 비례한다.

# 넓이로 직접 확인 : (r, theta) 의 작은 사각형이 얼마나 커지는가
h = 1e-3
for r in (1.0, 3.0):
    네점 = [극좌표([r, 0.0]), 극좌표([r + h, 0.0]),
           극좌표([r + h, h]), 극좌표([r, h])]
    두변 = np.column_stack([네점[1] - 네점[0], 네점[3] - 네점[0]])
    실제 = abs(np.linalg.det(두변))
    print(f"r = {r:g} :  옮겨진 넓이 {실제:.3e},   h^2 곱하기 r = {h * h * r:.3e}")

r = 1 :  옮겨진 넓이 1.000e-06,   h^2 곱하기 r = 1.000e-06
r = 3 :  옮겨진 넓이 3.000e-06,   h^2 곱하기 r = 3.000e-06

마치며...¶

서술 파트의 내용	이 노트북의 코드
$AC^{\mathsf{T}} = (\det A)I$	`여인수행렬(A)` 를 만들어 곱한다
비대각이 0인 이유	행을 갈아끼운 행렬의 행렬식과 대조
$A^{-1} = C^{\mathsf{T}}/\det A$	`딸림역행렬(A)`
크래머 공식	`크래머(A, b)` — 값은 맞고 속도와 정확도는 진다
민감도	$\det A$ 를 줄이며 $\vv{b}$ 를 흔든다
부피	슬라이더로 $\det$ 를 2에서 0으로
$\det = 0$ 의 네 얼굴	부피·랭크·영공간·역행렬을 한 표에
야코비안	수치미분으로 극좌표의 $r$

더 해 볼 것¶

딸림역행렬 과 np.linalg.inv 의 속도를 $n$ 을 키우며 재 보자. 여인수를 $n^2$ 개 구해야 하므로 얼마나 나빠지는가?
4절의 슬라이더에서 세 번째 열 대신 첫 번째 열을 밀어 보자. 같은 일이 일어나는가?
5절의 반사행렬을 두 번 곱하면 행렬식이 얼마인가? 왜 그런가?
3차원 극좌표(구면좌표)의 야코비안을 야코비행렬 로 구해 보자. 교과서의 $r^2\sin\phi$ 와 맞는가?

다음 강의에서 $\det(A - \lambda I) = 0$ 을 만난다. 이제 이 식을 부피가 0이 되도록 $\lambda$ 를 고른다로 읽을 수 있다. 세 강의를 쓴 이유가 거기서 드러난다.

Lecture 20. 크래머 공식, 역행렬, 부피 — 파이썬 실습