배송 허브 세 곳의 물량 맞추기: 가우스 소거와 피벗

입사 둘째 주 아침, 민재 팀장이 노트북을 들고 관제실로 들어옵니다. "다음 주부터 허브 세 곳에서 비행을 나눠 띄워요. 드론이 한 번 뜨면 여러 마을에 상자를 나눠 내려놓거든요. 세 지역 수요가 딱 맞게 허브마다 몇 번씩 띄울지 정해야 하는데, 엑셀에서 하나를 맞추면 다른 하나가 틀어져요. 숫자가 이상해요!"

서윤이 표를 들여다봅니다. "미지수 셋에 식 셋이네요. 중학교 때 푼 연립방정식이에요. 식끼리 빼서 미지수를 하나씩 지우면 되죠." 종이에 10분쯤 끄적이더니 답을 내밉니다. 항구 6번, 읍내 8번, 산마을 4번. 넣어 보니 세 지역 수요가 정확히 맞습니다.

오 박사가 종이를 보고 고개를 끄덕입니다. "방금 한 게 가우스 소거야. 컴퓨터도 똑같이 하지. 그런데 컴퓨터는 0에 아주 가까운 수로 나누는 순간 오류도 없이 틀린 답을 내. 손으로 풀 때는 절대 안 하는 실수야. 오늘은 그걸 피하는 법까지 가자."

이번 교시에는 손으로 푼 방법을 그대로 코드로 옮기고, 그 코드가 망가지는 두 가지 경우(아주 작은 피벗, 해가 없는 행렬)를 직접 만들어 봅니다. 선형대수 강의의 8교시 연립방정식 Ax = b가 "답이 있나, 몇 개인가"를 그림으로 봤다면, 이번 교시는 "컴퓨터로 어떻게, 얼마나 정확하게 구하나"를 봅니다.

7-1. 허브 물량을 식으로: Ax = b (6분)

드론 한 대는 작은 약 상자를 8~10개 싣고 떠서, 정해진 길을 돌며 여러 지역에 나눠 내려놓습니다. 비행 한 번이 지역마다 내려놓는 상자 수는 길마다 정해져 있습니다.

허브별 하루 비행 횟수를 $x_1, x_2, x_3$이라 하면, 지역마다 받는 상자 수가 수요와 같아야 하므로 식이 세 개 나옵니다.

계수만 모아 행렬로 적으면 $A\mathbf{x} = \mathbf{b}$입니다. $A$는 계수 행렬, $\mathbf{x}$는 미지수 벡터, $\mathbf{b}$는 오른쪽 변(수요)입니다.

$A$의 1열 (4, 2, 2)는 "항구 비행 한 번이 세 지역에 내려놓는 상자"입니다. 그러니 $A\mathbf{x}$는 "항구 비행 묶음 $x_1$번 + 읍내 비행 묶음 $x_2$번 + 산마을 비행 묶음 $x_3$번"이 세 지역에 내려놓는 상자 수입니다. 행렬 곱을 열의 일차결합으로 보는 이 그림은 선형대수 강의 5교시에서 다뤘습니다.

소거를 할 때는 $A$ 옆에 $\mathbf{b}$를 붙인 확대행렬(augmented matrix) $[A \mid \mathbf{b}]$ 하나만 다룹니다. 식 하나가 행 하나입니다.

>>> A = np.array([[4, 2, 1],
...               [2, 5, 2],
...               [2, 3, 5]], dtype=float)
>>> b = np.array([44, 60, 56], dtype=float)
>>> A.shape, b.shape
((3, 3), (3,))
>>> Ab = np.column_stack([A, b])
>>> Ab
array([[ 4.,  2.,  1., 44.],
       [ 2.,  5.,  2., 60.],
       [ 2.,  3.,  5., 56.]])

dtype=float를 꼭 붙였습니다. 이유는 '자주 하는 실수'에서 봅니다. np.column_stack은 배열들을 열로 나란히 붙입니다.

7-2. 전진 소거와 후진 대입: 손으로 한 번 (10분)

서윤이 종이에 한 일을 순서대로 적어 봅시다. 쓰는 도구는 하나뿐입니다. 한 행에서 다른 행의 몇 배를 빼기. 식에서 식을 빼도 답은 그대로이므로, 이렇게 바꾼 연립방정식은 처음 것과 답이 같습니다.

전진 소거(forward elimination) 1단계. 1행의 첫 수 $a_{11} = 4$를 피벗(pivot) 으로 삼아, 그 아래 두 수를 0으로 만듭니다. 2행에서 1행의 $l_{21} = 2/4 = 0.5$배를 빼고, 3행에서도 1행의 $l_{31} = 2/4 = 0.5$배를 뺍니다. 이 $l$을 곱수(multiplier) 라고 합니다.

예를 들어 새 2행은 $(2, 5, 2, 60) - 0.5 \times (4, 2, 1, 44) = (0, 4, 1.5, 38)$입니다.

2단계. 이제 2행의 $4$가 피벗입니다. $l_{32} = 2/4 = 0.5$이므로 3행에서 2행의 0.5배를 뺍니다.

왼쪽 아래가 모두 0인 이 모양을 위삼각 행렬(upper triangular matrix) $U$라고 합니다. 여기까지가 전진 소거입니다.

후진 대입(back substitution). 맨 아래 식에는 미지수가 하나뿐입니다. 아래에서 위로 올라가며 하나씩 풉니다.

항구 6번, 읍내 8번, 산마을 4번입니다. 첫 식에 넣으면 $4 \cdot 6 + 2 \cdot 8 + 4 = 44$로 맞습니다.

크기가 $n$인 일반적인 경우로 적으면 이렇습니다. $k$번째 단계에서 피벗은 $a_{kk}$이고, 그 아래 행 $i$마다

후진 대입은 $x_i = \Bigl(b_i - \sum_{j > i} u_{ij} x_j\Bigr) / u_{ii}$입니다. 피벗 $a_{kk}$로 나누는 곳이 두 군데라는 점을 기억해 두세요. 7-5의 사고가 바로 여기서 납니다.

같은 계산을 NumPy로 한 줄씩 해 봅니다. M[1]은 2행 전체이고, M[1] -= l21 * M[0]이 "R2 ← R2 − l21·R1"입니다.

>>> M = Ab.copy()
>>> l21 = M[1, 0] / M[0, 0]
>>> l31 = M[2, 0] / M[0, 0]
>>> l21, l31
(np.float64(0.5), np.float64(0.5))
>>> M[1] -= l21 * M[0]
>>> M[2] -= l31 * M[0]
>>> M
array([[ 4. ,  2. ,  1. , 44. ],
       [ 0. ,  4. ,  1.5, 38. ],
       [ 0. ,  2. ,  4.5, 34. ]])
>>> l32 = M[2, 1] / M[1, 1]
>>> M[2] -= l32 * M[1]
>>> M
array([[ 4.  ,  2.  ,  1.  , 44.  ],
       [ 0.  ,  4.  ,  1.5 , 38.  ],
       [ 0.  ,  0.  ,  3.75, 15.  ]])
>>> x3 = M[2, 3] / M[2, 2]
>>> x2 = (M[1, 3] - M[1, 2] * x3) / M[1, 1]
>>> x1 = (M[0, 3] - M[0, 1] * x2 - M[0, 2] * x3) / M[0, 0]
>>> x1, x2, x3
(np.float64(6.0), np.float64(8.0), np.float64(4.0))

아래 실습대는 확대행렬 [A | b]를 칸마다 보여 줍니다. '처음'을 누른 뒤 '다음'으로 한 단계씩 가 보세요. 호박색 칸이 지금의 피벗, 남색으로 칠해진 행이 지금 바뀌는 행, 회색 행이 빼 주는 쪽 행입니다. 칸 아래 상자에 지금 한 연산(예: R2 ← R2 − 0.5 · R1)과 곱수 계산이 나옵니다. 소거가 끝나면 후진 대입으로 x₃, x₂, x₁이 차례로 채워집니다. 허브 문제에서는 오른쪽 슬라이더로 지역별 수요 b를 바꿀 수 있습니다.

7-3. 파이썬으로 짜기: naive_gauss와 np.linalg.solve (8분)

손 계산을 반복문 두 겹으로 일반화합니다. 바깥 k는 피벗 열, 안쪽 i는 고칠 행입니다. 피벗을 고르는 장치가 아직 없어서 "순진한(naive)" 가우스 소거라고 부릅니다.

>>> def naive_gauss(A, b):
...     A = A.astype(float)
...     b = b.astype(float)
...     n = len(b)
...     for k in range(n - 1):
...         for i in range(k + 1, n):
...             l = A[i, k] / A[k, k]
...             A[i, k:] -= l * A[k, k:]
...             b[i] -= l * b[k]
...     x = np.zeros(n)
...     for i in range(n - 1, -1, -1):
...         x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
...     return x
>>> naive_gauss(A, b)
array([6., 8., 4.])
>>> np.linalg.solve(A, b)
array([6., 8., 4.])

몇 군데를 짚어 봅니다.

• A.astype(float)는 새 배열을 만들어 돌려줍니다. 그래서 함수 안에서 A를 고쳐도 바깥의 A는 그대로입니다. 이 줄이 없으면 함수를 부를 때마다 원래 행렬이 망가집니다.
• A[i, k:]는 $i$행의 $k$열부터 끝까지입니다. 그 왼쪽은 이미 0이라 계산할 필요가 없습니다.
• 후진 대입의 $\sum_{j>i} u_{ij}x_j$는 A[i, i + 1:] @ x[i + 1:] 한 번으로 계산합니다. 맨 아래 행에서는 빈 배열끼리의 곱이라 0이 됩니다.

답을 받았으면 반드시 넣어 봅니다. $\mathbf{r} = \mathbf{b} - A\mathbf{x}$를 잔차(residual) 라고 합니다. 정확한 답이면 0 벡터이고, 부동소수점 계산에서는 0에 아주 가까워야 합니다. 비교는 == 대신 np.allclose로 합니다(2교시).

>>> b - A @ naive_gauss(A, b)
array([0., 0., 0.])
>>> rng = np.random.default_rng(7)
>>> B = rng.normal(size=(200, 200))
>>> c = rng.normal(size=200)
>>> x = naive_gauss(B, c)
>>> np.allclose(x, np.linalg.solve(B, c))
True
>>> np.allclose(B @ x, c)
True
>>> np.abs(c - B @ x).max() < 1e-10
np.True_

무작위 200 × 200 연립방정식도 np.linalg.solve와 같은 답을 냅니다. 잔차의 가장 큰 값도 $10^{-10}$보다 작습니다. 그렇다면 naive_gauss로 충분할까요? 7-5에서 이 함수가 망가지는 3 × 3 예제를 만들어 봅니다.

7-4. 얼마나 오래 걸리나: 연산량 2n³/3 (5분)

$k$번째 단계에서는 아래 행 $n - k - 1$개마다 길이 약 $n - k$인 행을 곱하고 뺍니다. 곱셈과 뺄셈을 따로 세면 한 행에 약 $2(n-k)$번입니다. 모두 더하면 이렇습니다.

후진 대입은 $n^2$번 정도라서, $n$이 크면 전진 소거에 묻힙니다. 오른쪽 변까지 넣어 정확히 세는 함수로 확인해 봅니다.

>>> def flops(n):
...     k = np.arange(n - 1)
...     return int(((n - 1 - k) * (2 * (n - k) + 3)).sum())
>>> for n in [3, 10, 100, 1000, 10000]:
...     print(n, flops(n), round(flops(n) / (2 * n**3 / 3), 3))
3 25 1.389
10 795 1.193
100 681450 1.022
1000 668164500 1.002
10000 666816645000 1.0

$n$이 커질수록 비율이 1에 붙습니다. 연산량이 $O(n^3)$이라는 뜻은 $n$이 2배면 시간이 8배, 10배면 1000배라는 뜻입니다.

np.linalg.solve는 C와 포트란으로 짠 LAPACK을 불러서 파이썬 반복문인 naive_gauss보다 수백 배 빠릅니다. 그래도 $n^3$은 피할 수 없습니다. 직접 재 보려면 아래 코드를 돌려 보세요. 결과는 컴퓨터마다 다릅니다.

import time
rng = np.random.default_rng(0)
for n in [1000, 2000, 4000]:
    M = rng.normal(size=(n, n))
    v = rng.normal(size=n)
    t0 = time.perf_counter()
    np.linalg.solve(M, v)
    print(n, round(time.perf_counter() - t0, 3), "초")

$n$이 커지면 다른 길이 필요합니다. 같은 행렬로 여러 번 풀 때 소거를 한 번만 하는 법, 대부분이 0인 행렬을 빠르게 푸는 법은 8교시의 주제입니다.

7-5. 아주 작은 피벗이 부르는 재앙과 부분 피벗팅 (9분)

이제 naive_gauss를 망가뜨려 봅시다. 피벗 자리에 아주 작은 수가 있는 2 × 2 연립방정식입니다.

둘째 식에서 첫째 식을 빼 보면 답은 거의 정확히 $x_1 = 1,\ x_2 = 1$입니다. 그런데 naive_gauss는 이렇게 답합니다.

>>> E = np.array([[1e-17, 1], [1, 1]])
>>> e = np.array([1.0, 2.0])
>>> naive_gauss(E, e)
array([0., 1.])
>>> np.linalg.solve(E, e)
array([1., 1.])

$x_1 = 0$입니다. 오류도 경고도 없이 틀렸습니다. 손으로 따라가 보면 이유가 보입니다. 곱수가 $l_{21} = 1/10^{-17} = 10^{17}$이고, 2행은 이렇게 바뀝니다.

float64는 유효숫자가 16자리쯤입니다(2교시). $10^{17}$ 옆의 1이나 2는 그 자릿수 밖이라 반올림에서 사라집니다. 2교시의 "큰 수 + 작은 수"입니다.

>>> 1 - 1e17 == -1e17
True
>>> 2 - 1e17 == -1e17
True

2행이 $-10^{17}x_2 = -10^{17}$이 되어 $x_2 = 1$은 맞게 나옵니다. 하지만 2행이 원래 갖고 있던 정보("$x_1 + x_2 = 2$")는 지워졌습니다. 마지막에 $x_1 = (1 - 1)/10^{-17} = 0$이 됩니다. 작은 피벗으로 나누면 곱수가 거대해지고, 거대한 수를 곱해 빼는 순간 원래 행의 정보가 반올림에 먹힙니다.

해결책은 간단합니다. 피벗 자리에 그 열에서 절댓값이 가장 큰 수가 오도록 행을 바꿉니다. 이것을 부분 피벗팅(partial pivoting) 이라고 합니다. 그러면 곱수가 언제나 $|l_{ik}| \le 1$이라, 빼는 행이 원래 행을 덮어 버리는 일이 없습니다. 행을 바꾸는 것은 식의 순서를 바꾸는 것뿐이라 답은 그대로입니다.

>>> def gauss_pivot(A, b):
...     A = A.astype(float)
...     b = b.astype(float)
...     n = len(b)
...     for k in range(n - 1):
...         p = k + np.argmax(np.abs(A[k:, k]))
...         if p != k:
...             A[[k, p]] = A[[p, k]]
...             b[[k, p]] = b[[p, k]]
...         for i in range(k + 1, n):
...             l = A[i, k] / A[k, k]
...             A[i, k:] -= l * A[k, k:]
...             b[i] -= l * b[k]
...     x = np.zeros(n)
...     for i in range(n - 1, -1, -1):
...         x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
...     return x
>>> gauss_pivot(E, e)
array([1., 1.])

np.argmax(np.abs(A[k:, k]))는 $k$열의 $k$행 아래에서 절댓값이 가장 큰 수가 몇 번째인지 돌려줍니다. 거기에 k를 더하면 행 번호가 됩니다. A[[k, p]] = A[[p, k]]는 두 행을 맞바꿉니다.

허브 문제에서도 같은 일이 생길 수 있습니다. 어느 날 항구 드론이 정비에 들어가 항구 지역 배달을 거의 하지 않는다고 해 봅시다. $a_{11}$이 정확히 0이면 0으로 나누기에서 바로 nan이 나와 차라리 눈에 띕니다. 더 위험한 것은 앞선 계산의 반올림 찌꺼기 때문에 0 대신 $10^{-16}$ 같은 수가 남는 경우입니다.

>>> T = np.array([[1e-16, 2, 1],
...               [2, 5, 2],
...               [2, 3, 5]])
>>> t = np.array([20.0, 60, 56])
>>> naive_gauss(T, t)
array([0., 6., 8.])
>>> gauss_pivot(T, t)
array([6., 8., 4.])
>>> np.linalg.solve(T, t)
array([6., 8., 4.])

답 (6, 8, 4)가 (0, 6, 8)로 바뀌었습니다. 세 값이 모두 그럴듯한 정수라서 결과만 봐서는 틀린 줄 모릅니다.

아래 실습대에서 '작은 피벗 예제'를 누르세요. 1행 1열의 1e-16이 피벗입니다. '부분 피벗팅'이 꺼져 있으면 곱수가 2e16이 되어 R2, R3에 −4e16 같은 거대한 수가 생깁니다. 체크하면 첫 단계에서 R1 ↔ R2 바꾸기(청록 행)가 먼저 일어납니다. 오른쪽 아래에 참값과 상대 오차가 나옵니다.

7-6. 해가 없을 때: 특이 행렬과 잔차 확인 (5분)

이번에는 산마을 드론이 모두 고장 나서, 산마을 비행을 항구 드론이 항구 비행과 똑같은 길로 대신 뛴다고 해 봅시다. 그러면 3열이 1열과 같아집니다.

소거하면 마지막 행이 $[\,0\ \ 0\ \ 0 \mid 15\,]$, 곧 "$0 \cdot x_3 = 15$"가 됩니다. 어떤 $x_3$도 이 식을 맞출 수 없으니 해가 없습니다. 오른쪽이 0이었다면 $x_3$을 아무 값으로 골라도 되어 해가 무수히 많았을 것입니다. 어느 쪽이든 해가 하나로 정해지지 않습니다. 이런 행렬을 특이 행렬(singular matrix) 이라고 합니다. 행을 바꿔도 0이 아닌 피벗을 찾을 수 없다는 것이 소거로 보는 특이 행렬의 모습입니다.

>>> S = np.array([[4, 2, 4],
...               [2, 5, 2],
...               [2, 3, 2]], dtype=float)
>>> np.linalg.det(S)
np.float64(0.0)
>>> np.linalg.matrix_rank(S)
np.int64(2)
>>> np.linalg.solve(S, b)
Traceback (most recent call last):
  ...
numpy.linalg.LinAlgError: Singular matrix

행렬식이 0이고 계수(rank)가 2입니다. 열 두 개가 같으니, 세 열을 아무리 섞어도 3차원 공간을 다 채우지 못하고 평면 하나만 만듭니다. 수요 (44, 60, 56)이 그 평면 밖에 있어서 해가 없습니다. 선형대수 강의 7교시 행렬식과 9교시 계수에서 본 그림 그대로입니다.

이번에는 NumPy가 LinAlgError로 알려 줬습니다. 소거 중에 피벗이 정확히 0이 되었기 때문입니다. 부동소수점에서는 0이어야 할 피벗이 $10^{-16}$ 같은 찌꺼기로 남는 일이 흔합니다. 그러면 오류 없이 엉뚱한 답이 나옵니다. 3행이 1행과 2행의 합인 행렬로 해 봅시다.

>>> S2 = np.array([[4, 2, 1],
...                [2, 5, 2],
...                [6, 7, 3]], dtype=float)
>>> c2 = np.array([44.0, 60, 100])
>>> x2 = np.linalg.solve(S2, c2)
>>> np.abs(x2).max() > 1e10
np.True_
>>> np.allclose(S2 @ x2, c2)
False

오류는 없지만 답의 크기가 $10^{15}$ 근처이고, 넣어 보면 식이 맞지 않습니다. 그래서 답을 받으면 잔차 확인(np.allclose(A @ x, b)) 을 습관으로 합니다. 잔차가 작아도 답이 크게 틀릴 수 있는 경우가 있는데, 그 이야기는 8교시의 조건수에서 합니다.

책의 기호 읽는 법

실무에서는

np.linalg.solve(A, b)는 LAPACK의 gesv를 부릅니다. 하는 일은 이번 교시의 gauss_pivot과 같습니다(부분 피벗팅을 쓰는 LU 분해, 8교시). 직접 짠 함수보다 빠르고, 특이 행렬도 알려 줍니다.

x = np.linalg.solve(A, b)          # 기본
X = np.linalg.solve(A, B)          # B의 열마다 한 번씩: 수요 여러 개를 한꺼번에
np.allclose(A @ x, b)              # 잔차 확인

# 하지 말 것: 역행렬을 만들어 곱하기
x = np.linalg.inv(A) @ b           # 계산이 약 3배 더 들고, 오차도 더 크다

# SciPy (이 사이트 서버에는 없음)
import scipy.linalg
scipy.linalg.solve(A, b, assume_a="pos")   # 대칭·양의 정부호 행렬이면 더 빠른 방법(촐레스키)

inv(A) @ b도 작은 행렬에서는 같은 답이 나와서 습관이 되기 쉽습니다. 하지만 역행렬은 $\mathbf{b}$가 $n$개일 때 연립방정식 $n$개를 푸는 것과 같은 일이라 일이 많고, 반올림 오차도 한 번 더 쌓입니다. 수학책의 $A^{-1}\mathbf{b}$는 "연립방정식을 풀어라"로 읽으세요.

자주 하는 실수 (3분)

정수 배열로 소거하기. np.array([[4, 2], [2, 5]])처럼 숫자를 정수로만 적으면 dtype이 int64입니다. 여기에 A[1] -= 0.5 * A[0]을 하면 NumPy가 실수 결과를 정수 칸에 넣을 수 없다며 오류를 냅니다. 오류 대신 소수점이 잘려 나가는 코드도 있으니, 행렬을 만들 때 dtype=float를 붙이거나 astype(float)로 바꾸세요.

>>> np.array([[4, 2], [2, 5]]).dtype
dtype('int64')

행 바꾸기를 파이썬 식으로 쓰기. 리스트에서 하던 대로 M[0], M[1] = M[1], M[0]이라고 쓰면 두 행이 같아집니다. M[1]은 복사본이 아니라 원래 배열을 들여다보는 창(view)이라, 첫 대입이 끝나는 순간 오른쪽의 M[0]도 바뀌기 때문입니다. 행 번호 목록으로 한 번에 바꾸세요.

>>> M = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> M[0], M[1] = M[1], M[0]
>>> M
array([[3., 4.],
       [3., 4.]])
>>> M = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> M[[0, 1]] = M[[1, 0]]
>>> M
array([[3., 4.],
       [1., 2.]])

답을 넣어 보지 않기. 피벗팅 없는 소거, 거의 특이한 행렬 모두 오류 없이 틀린 답을 냅니다. np.allclose(A @ x, b) 한 줄이면 대부분 걸러집니다.

연습 문제 (4분)

1. 항구 지역 수요가 44에서 50으로 늘었습니다. 허브별 비행 횟수를 구하고, 정수가 아니면 가장 가까운 정수로 반올림했을 때 지역별 상자 수가 수요와 얼마나 어긋나는지 보세요.
2. 피벗이 정확히 0인 [[0, 1], [1, 1]]로 naive_gauss를 부르면 무슨 일이 일어날까요? np.linalg.solve는 어떤가요?
3. 내 컴퓨터에서 1000 × 1000 연립방정식을 푸는 데 0.05초가 걸렸습니다. 2000 × 2000은 대략 얼마나 걸릴까요?
4. 특이 행렬 S(7-6)에서 항구·읍내 수요는 44, 60 그대로일 때, 산마을 수요가 얼마이면 해가 무수히 많아질까요? 해 두 개를 찾아 확인하세요.

풀이

>>> x = np.linalg.solve(A, np.array([50.0, 60, 56]))
>>> x
array([7.9, 7.4, 3.6])
>>> A @ np.array([8, 7, 4])
array([50., 59., 57.])
>>> Z = np.array([[0.0, 1], [1, 1]])
>>> with np.errstate(all="ignore"):
...     print(naive_gauss(Z, np.array([1.0, 2])))
[nan nan]
>>> np.linalg.solve(Z, np.array([1.0, 2]))
array([1., 1.])
>>> 2000**3 / 1000**3
8.0
>>> M = np.column_stack([S, [44.0, 60, 41]])
>>> M[1] -= 0.5 * M[0]
>>> M[2] -= 0.5 * M[0]
>>> M[2] -= 0.5 * M[1]
>>> M
array([[ 4.,  2.,  4., 44.],
       [ 0.,  4.,  0., 38.],
       [ 0.,  0.,  0.,  0.]])
>>> S @ np.array([6.25, 9.5, 0])
array([44., 60., 41.])
>>> S @ np.array([2.25, 9.5, 4])
array([44., 60., 41.])

1번은 (7.9, 7.4, 3.6)으로 정수가 아닙니다. (8, 7, 4)로 반올림하면 읍내가 1상자 모자라고 산마을이 1상자 남습니다. 비행 횟수처럼 정수여야 하는 답을 정확히 맞추는 일은 정수 계획법이라는 다른 분야의 문제입니다. 현장에서는 이렇게 반올림한 뒤 남는 상자를 다음 비행에 싣는 식으로 맞춥니다.

2번은 첫 곱수가 1/0 = inf가 되어 모든 값이 nan이 됩니다. np.errstate는 그 경고를 잠시 끄는 장치입니다. np.linalg.solve는 부분 피벗팅으로 두 행을 바꿔 바로 (1, 1)을 냅니다. 행렬은 멀쩡한데 피벗 자리만 0이었던 것입니다.

3번은 $n^3$에 비례하므로 8배, 약 0.4초입니다. 실제로는 캐시와 병렬 처리 때문에 정확히 8배가 되지는 않습니다.

4번은 마지막 행의 오른쪽이 0이 되어야 합니다. 손으로 소거하면 그 값은 $b_3 - 0.25,b_1 - 0.5,b_2 = b_3 - 41$이므로 $b_3 = 41$입니다. 이때 $x_2 = 38/4 = 9.5$로 정해지고, $x_1 + x_3 = 6.25$를 만족하는 짝은 모두 답입니다. 1열과 3열이 같으니 항구 비행과 산마을 비행을 얼마로 나누든 결과가 같은 것입니다.

저녁 무렵, 민재 팀장이 다시 들어옵니다. "허브 표는 잘 돌아가요. 그런데 다음 달부터 도시 전체 충전소 마흔 곳의 전력 배분도 이렇게 맞춰야 한대요. 수요가 매시간 바뀌니까 하루에 스물네 번씩요." 서윤이 계산해 봅니다. 40 × 40 소거를 하루 24번, 그리고 충전소는 계속 늘어날 겁니다. 다음 교시에는 같은 행렬로 여러 번 푸는 법(LU 분해), 대부분이 0인 큰 행렬을 조금씩 고쳐 가며 푸는 법(반복법), 그리고 답을 얼마나 믿을 수 있는지 알려 주는 조건수를 배웁니다.