Lecture 19. 행렬식 공식과 여인수 — 파이썬 실습 - 선형대수

L19 서술 파트에서 두 가지 공식을 유도했다. $n!$ 개의 항으로 된 빅 포뮬러와, 그것을 한 행 기준으로 묶은 여인수 전개이다. 그리고 둘 다 실제 계산에는 쓸 수 없다고 했다.

이 노트북에서는 두 공식을 직접 짜서 np.linalg.det 와 값이 같은지 확인하고, $n$ 을 키워 가며 왜 쓸 수 없는지 시계로 잰다. 값이 맞는데도 쓸 수 없다는 것이 요점이다.

서술 파트의 내용	여기서 확인하는 방법
항 하나 $=$ 순열 하나	`itertools.permutations` 로 전부 돌린다
부호 $=$ 교환 횟수의 홀짝	순열을 제자리로 되돌리며 센다
여인수 전개는 재귀	재귀 함수 다섯 줄
계산량 $n!$ 대 $n^3$	시계로 재고 $n = 20$ 까지 외삽
성질이 공식에서 보인다	살아남는 항의 개수를 센다

0. 준비¶

import math
import time
from itertools import permutations

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, show_matrix

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

numpy 2.5.2

1. 순열과 그 부호¶

먼저 순열 하나의 부호를 구한다. 순열을 제자리로 되돌리는 데 교환이 몇 번 필요한지 세면 된다.

def 순열부호(순열):
    """순열을 (0, 1, ..., n-1) 로 되돌리는 교환 횟수의 홀짝을 돌려준다."""
    남은 = list(순열)
    교환 = 0
    for i in range(len(남은)):
        while 남은[i] != i:                 # 제자리에 올 때까지 맞바꾼다
            j = 남은[i]
            남은[i], 남은[j] = 남은[j], 남은[i]
            교환 += 1
    return (-1) ** 교환

for 순열 in [(0, 1, 2, 3), (1, 0, 2, 3), (1, 2, 3, 0), (3, 2, 1, 0)]:
    보기 = " ".join(str(j + 1) for j in 순열)
    print(f"열 순서 {보기}   부호 {순열부호(순열):+d}")

열 순서 1 2 3 4   부호 +1
열 순서 2 1 3 4   부호 -1
열 순서 2 3 4 1   부호 -1
열 순서 4 3 2 1   부호 +1

서술 파트의 그림에 나온 세 가지와 맞는다. $n = 4$ 이면 순열이 $4! = 24$ 개이고, 부호가 $+$ 인 것과 $-$ 인 것이 정확히 반씩이다.

부호들 = [순열부호(p) for p in permutations(range(4))]
print("순열의 개수 :", len(부호들), " (4! =", math.factorial(4), ")")
print("부호 +1 :", 부호들.count(1), "개,   부호 -1 :", 부호들.count(-1), "개")

순열의 개수 : 24  (4! = 24 )
부호 +1 : 12 개,   부호 -1 : 12 개

2. 빅 포뮬러¶

각 행에서 열을 하나씩, 겹치지 않게 고르고 그 성분들을 곱한다. 그것을 $n!$ 번 한다.

def 빅포뮬러(A):
    """정의 그대로 n! 개의 항을 전부 더한다. 실전에서는 절대 쓰지 말 것."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    합 = 0.0
    for 순열 in permutations(range(n)):
        곱 = 1.0
        for i, j in enumerate(순열):        # i 행에서 j 열을 골랐다
            곱 *= A[i, j]
        합 += 순열부호(순열) * 곱
    return 합

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

print(show_matrix(A, "A"))
print("빅 포뮬러     :", 빅포뮬러(A))
print("np.linalg.det :", np.linalg.det(A))
print("같은가 :", np.isclose(빅포뮬러(A), np.linalg.det(A)))

A
[  1   2   1 ]
[  3   8   1 ]
[  0   4   1 ]
빅 포뮬러     : 10.0
np.linalg.det : 9.999999999999998
같은가 : True

항 여섯 개를 하나씩 펼쳐 보자. 어느 항이 살아남고 어느 항이 0인지 보인다.

print(f"{'열 순서':>10}{'부호':>6}{'성분들':>22}{'항':>10}")
합계 = 0.0
for 순열 in permutations(range(3)):
    성분 = [A[i, j] for i, j in enumerate(순열)]
    항 = 순열부호(순열) * float(np.prod(성분)) + 0.0      # -0.0 을 0.0 으로
    합계 += 항
    보기 = " x ".join(f"{v:g}" for v in 성분)
    print(f"{' '.join(str(j + 1) for j in 순열):>10}"
          f"{순열부호(순열):>+6d}{보기:>22}{항:>10.0f}")
print(f"{'합계':>48}{합계:>10.0f}")

      열 순서    부호                   성분들         항
     1 2 3    +1             1 x 8 x 1         8
     1 3 2    -1             1 x 1 x 4        -4
     2 1 3    -1             2 x 3 x 1        -6
     2 3 1    +1             2 x 1 x 0         0
     3 1 2    +1             1 x 3 x 4        12
     3 2 1    -1             1 x 8 x 0         0
                                              합계        10

3. 여인수 전개¶

첫 행에서 하나를 고르고, 그 행과 열을 지우고, 남은 것의 행렬식을 다시 같은 방법으로 구한다. 재귀 함수 다섯 줄이면 된다.

def 여인수전개(A):
    """첫 행을 기준으로 재귀 전개한다. 이것도 실전용이 아니다."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    if n == 1:
        return float(A[0, 0])
    합 = 0.0
    for j in range(n):
        소행렬 = np.delete(np.delete(A, 0, axis=0), j, axis=1)   # 1행 j열을 지운다
        합 += (-1) ** j * A[0, j] * 여인수전개(소행렬)
    return 합

print("여인수 전개   :", 여인수전개(A))
print("빅 포뮬러     :", 빅포뮬러(A))
print("np.linalg.det :", np.linalg.det(A))

여인수 전개   : 10.0
빅 포뮬러     : 10.0
np.linalg.det : 9.999999999999998

여인수를 하나씩 꺼내 보자. 서술 파트에서 $C_{11} = 4$ , $C_{12} = -3$ , $C_{13} = 12$ 였다.

def 여인수(A, i, j):
    """(-1)^(i+j) 를 붙인 소행렬식."""
    소행렬 = np.delete(np.delete(np.asarray(A, dtype=float), i, axis=0), j, axis=1)
    return (-1) ** (i + j) * np.linalg.det(소행렬)

for j in range(3):
    print(f"a1{j + 1} = {A[0, j]:g},   C1{j + 1} = {여인수(A, 0, j):+.0f},"
          f"   곱 = {A[0, j] * 여인수(A, 0, j):+.0f}")
print()
print("합 :", sum(A[0, j] * 여인수(A, 0, j) for j in range(3)))

a11 = 1,   C11 = +4,   곱 = +4
a12 = 2,   C12 = -3,   곱 = -6
a13 = 1,   C13 = +12,   곱 = +12

합 : 9.999999999999998

어느 행으로 전개해도 같은 값이 나온다. 둘째 행과 셋째 행으로도 해 보자.

for i in range(3):
    값 = sum(A[i, j] * 여인수(A, i, j) for j in range(3))
    print(f"{i + 1} 행으로 전개 : {값:.4f}")
print()
for j in range(3):
    값 = sum(A[i, j] * 여인수(A, i, j) for i in range(3))
    print(f"{j + 1} 열로 전개 : {값:.4f}")

1 행으로 전개 : 10.0000
2 행으로 전개 : 10.0000
3 행으로 전개 : 10.0000

1 열로 전개 : 10.0000
2 열로 전개 : 10.0000
3 열로 전개 : 10.0000

4. 성질이 공식에서 보인다¶

빅 포뮬러를 보면 L18의 성질들이 왜 성립하는지가 눈에 들어온다. 살아남는 항이 몇 개인지 세어 보자.

def 살아있는항(A):
    """0 이 아닌 항의 개수를 센다."""
    A = np.asarray(A, dtype=float)
    n = A.shape[0]
    return sum(1 for 순열 in permutations(range(n))
               if not np.isclose(np.prod([A[i, j] for i, j in enumerate(순열)]), 0.0))

n = 5
가득 = rng.integers(1, 5, (n, n)).astype(float)
상삼각 = np.triu(가득)
영행 = 가득.copy(); 영행[2] = 0.0

print(f"{'':>14}{'전체 항':>10}{'0 이 아닌 항':>14}{'det':>12}")
for 이름, M in (("아무 행렬", 가득), ("상삼각", 상삼각), ("한 행이 0", 영행)):
    print(f"{이름:>14}{math.factorial(n):>10}{살아있는항(M):>14}"
          f"{np.linalg.det(M):>12.1f}")

                    전체 항      0 이 아닌 항         det
         아무 행렬       120           120       -56.0
           상삼각       120             1       384.0
        한 행이 0       120             0         0.0

상삼각행렬에서는 살아남는 항이 하나뿐이다. 0을 피해 고르는 방법이 대각을 따라가는 길밖에 없기 때문이고, 그래서 행렬식이 대각의 곱이 된다. 한 행이 0이면 모든 항에 0이 곱해져 살아남는 항이 없다.

print("상삼각의 대각 :", np.diag(상삼각))
print("대각의 곱     :", np.prod(np.diag(상삼각)))
print("det           :", np.linalg.det(상삼각))

상삼각의 대각 : [3. 4. 4. 2. 4.]
대각의 곱     : 384.0
det           : 384.0

같은 행이 둘이면 왜 0인가. 그 두 행을 맞바꾸는 순열끼리 짝을 지으면 계수는 같고 부호만 반대라 서로 상쇄된다. 짝을 지어 더해 보자.

같은행 = 가득.copy()
같은행[1] = 같은행[0]                       # 1행과 2행을 같게

짝합 = {}
for 순열 in permutations(range(5)):
    곱 = float(np.prod([같은행[i, j] for i, j in enumerate(순열)]))
    바꾼 = list(순열)
    바꾼[0], 바꾼[1] = 바꾼[1], 바꾼[0]      # 1행과 2행의 선택을 맞바꾼다
    열쇠 = tuple(sorted([순열, tuple(바꾼)]))
    짝합[열쇠] = 짝합.get(열쇠, 0.0) + 순열부호(순열) * 곱

print("짝의 개수 :", len(짝합))
print("짝마다의 합이 전부 0 인가 :",
      all(np.isclose(v, 0.0) for v in 짝합.values()))
print("전체 det :", np.linalg.det(같은행))

짝의 개수 : 60
짝마다의 합이 전부 0 인가 : True
전체 det : 0.0

5. 왜 쓸 수 없는가¶

값은 맞다. 문제는 시간이다. $n$ 을 하나씩 키우며 재 보자.

def 피벗곱(A):
    """소거해서 피벗을 곱한다 (L18 에서 만든 방법)."""
    U = np.asarray(A, dtype=float).copy()
    n = U.shape[0]
    교환 = 0
    for k in range(n):
        축 = k + int(np.argmax(np.abs(U[k:, k])))
        if abs(U[축, k]) < 1e-12:
            return 0.0
        if 축 != k:
            U[[k, 축]] = U[[축, k]]
            교환 += 1
        U[k + 1:, k:] -= np.outer(U[k + 1:, k] / U[k, k], U[k, k:])
    return (-1.0) ** 교환 * float(np.prod(np.diag(U)))


잰것 = []
print(f"{'n':>3}{'n!':>10}{'빅 포뮬러':>13}{'여인수':>12}{'피벗 곱':>12}{'numpy':>12}")
for n in range(2, 9):
    M = rng.normal(size=(n, n))
    걸린 = []
    for 함수 in (빅포뮬러, 여인수전개, 피벗곱, np.linalg.det):
        시작 = time.perf_counter()
        값 = 함수(M)
        걸린.append(time.perf_counter() - 시작)
    잰것.append((n, *걸린))
    print(f"{n:>3}{math.factorial(n):>10}"
          f"{걸린[0]:>13.5f}{걸린[1]:>12.5f}{걸린[2]:>12.5f}{걸린[3]:>12.5f}")
print()
print("(초 단위. 값은 네 방법 모두 같다.)")

  n        n!        빅 포뮬러         여인수        피벗 곱       numpy
  2         2      0.00002     0.00005     0.00015     0.00001
  3         6      0.00003     0.00012     0.00014     0.00001
  4        24      0.00006     0.00042     0.00032     0.00001
  5       120      0.00024     0.00191     0.00028     0.00001
  6       720      0.00152     0.01193     0.00048     0.00002

  7      5040      0.01195     0.07934     0.00017     0.00002

  8     40320      0.10444     0.62292     0.00024     0.00003

(초 단위. 값은 네 방법 모두 같다.)

$n = 8$ 에서 이미 빅 포뮬러가 소거보다 수천 배 느리다. 그러면 $n = 20$ 은 얼마나 걸리는가. $n = 8$ 의 항당 비용으로 외삽해 보자.

n8, 빅8 = 잰것[-1][0], 잰것[-1][1]
항당 = 빅8 / math.factorial(n8)
예상 = 항당 * math.factorial(20)

print(f"n = {n8} 에서 항 하나에 든 시간 : {항당 * 1e9:.1f} 나노초")
print(f"n = 20 의 항 개수            : {math.factorial(20):,}")
print(f"예상 소요 시간               : {예상:.3e} 초")
print(f"                             = {예상 / (365 * 24 * 3600):,.0f} 년")
print()
print(f"같은 n = 20 을 소거로 하면    : 곱셈 약 {20 ** 3 // 3:,} 번, 눈 깜짝할 사이")

n = 8 에서 항 하나에 든 시간 : 2590.4 나노초
n = 20 의 항 개수            : 2,432,902,008,176,640,000
예상 소요 시간               : 6.302e+12 초
                             = 199,838 년

같은 n = 20 을 소거로 하면    : 곱셈 약 2,666 번, 눈 깜짝할 사이

n값 = np.arange(2, 21)
팩 = np.array([math.factorial(int(v)) for v in n값], dtype=float)
세제곱 = n값.astype(float) ** 3
잰n = np.array([r[0] for r in 잰것])
잰빅 = np.array([r[1] for r in 잰것])

go.Figure(
    data=[
        go.Scatter(x=n값, y=팩, mode="lines", name="n!  (빅 포뮬러)",
                   line=dict(color=COLORS["error"], width=3)),
        go.Scatter(x=n값, y=세제곱, mode="lines", name="n^3  (소거)",
                   line=dict(color=COLORS["input"], width=3)),
        go.Scatter(x=잰n, y=잰빅 / 항당, mode="markers", name="실제로 잰 항의 개수",
                   marker=dict(color=COLORS["output"], size=9)),
    ],
    layout=go.Layout(
        title=dict(text="항의 개수 : n! 과 n^3"),
        xaxis=dict(title=dict(text="n")),
        yaxis=dict(title=dict(text="연산 횟수"), type="log"),
        height=460, margin=dict(l=70, r=20, t=60, b=50)))

주황 점은 실제로 잰 시간을 항당 비용으로 나눈 것이다. 빨간 곡선 위에 정확히 얹힌다. $n!$ 이 늘어나는 속도는 그래프를 로그로 그려야 겨우 담긴다.

마치며...¶

서술 파트의 내용	이 노트북의 코드
순열의 부호	`순열부호(순열)` — 제자리로 되돌리며 센다
빅 포뮬러	`빅포뮬러(A)` — $n!$ 개의 항을 전부
여인수 전개	`여인수전개(A)` — 재귀 다섯 줄
어느 행으로 전개해도 같다	세 행과 세 열로 각각 전개
상삼각은 항이 하나뿐	`살아있는항(M)` 이 1
같은 행이 둘이면 상쇄	짝을 지어 더하면 각 짝이 0
$n!$ 대 $n^3$	시계로 재고 $n = 20$ 으로 외삽

더 해 볼 것¶

여인수전개 에서 첫 행 대신 0이 가장 많은 행을 골라 전개하도록 고쳐 보자. 얼마나 빨라지는가? 그래도 $n!$ 인가?
살아있는항 을 하삼각행렬과 대각행렬에 대해 돌려 보자. 몇 개가 살아남는가?
크기가 같은 두 행렬 $A$ , $B$ 에 대해 $\det(AB)$ 를 빅 포뮬러로 구하는 것과 $\det A \cdot \det B$ 를 구하는 것 중 어느 쪽이 빠른가?
순열부호 를 쓰지 않고 순열행렬을 실제로 만들어 np.linalg.det 로 부호를 구해 보자. 같은 답이 나오는가?

다음 강의에서는 이 여인수들이 모여 역행렬의 공식을 이루는 것을 본다. 그리고 행렬식의 진짜 정체인 부피를 제대로 만난다.