지도가 커지면: LU 분해, 반복법, 조건수

다음 날 아침, 민재 팀장이 도시 지도를 화이트보드에 붙입니다. 해안 도로를 따라 드론 충전소 40곳이 한 줄로 이어져 있습니다. "충전소끼리 전력선이 이웃으로 연결돼 있어서 전력을 나눠 써요. 전력 팀 규칙대로 배분량을 정해야 하는데, 수요가 매시간 바뀌니까 하루에 스물네 번 풀어야 해요. 내년엔 충전소가 400곳이 된대요."

서윤은 어제 만든 np.linalg.solve를 반복문에 넣으면 된다고 생각합니다. 그러다 화면에 찍힌 40 × 40 행렬을 보고 멈칫합니다. 1,600칸 중 0이 아닌 칸이 118개뿐입니다. "이 0들까지 매번 다 소거하는 거예요?"

오 박사가 웃습니다. "좋은 질문이야. 낭비가 두 군데 있어. 행렬이 매번 같은데 소거를 처음부터 다시 하는 것, 그리고 0을 계산하는 것. 어제 곱수 l을 계산하고 버렸지? 그걸 적어 두면 첫 번째 낭비가 사라져. 두 번째는 아예 다른 방법으로 풀지. 그리고 하나 더. 답을 얻으면 그 답을 몇 자리까지 믿을 수 있는지도 말할 줄 알아야 해."

전력 팀의 규칙은 이렇습니다. 충전소 $i$가 받는 전력 $x_i$(kW)는 자기 수요 $d_i$의 1/3에, 양옆 이웃이 받는 전력의 1/3씩을 더한 값입니다. 식으로 쓰고 양변에 3을 곱하면 충전소마다 식이 하나씩 나옵니다.

양 끝 충전소는 이웃이 하나뿐입니다. 충전소 4곳이면 이런 연립방정식입니다. 대각선과 그 양옆에만 수가 있는 이 모양을 삼중대각 행렬(tridiagonal matrix) 이라고 합니다.

8-1. 한 번 분해, 여러 번 풀기: LU 분해 (10분)

어제의 허브 행렬로 시작합니다. 소거하며 쓴 곱수는 $l_{21} = 0.5$, $l_{31} = 0.5$, $l_{32} = 0.5$였습니다. 이 곱수를 대각선이 1인 아래삼각 행렬 $L$에 적어 두고, 소거가 끝난 위삼각 행렬을 $U$라 하면 놀라운 일이 생깁니다.

$L$의 2행과 $U$를 곱해 보면 $0.5 \times (4, 2, 1) + 1 \times (0, 4, 1.5) = (2, 5, 2)$로 $A$의 2행이 나옵니다. "2행에서 1행의 0.5배를 뺐더니 $U$의 2행이 되었다"를 거꾸로 적은 것입니다. 이것이 LU 분해(LU decomposition) 입니다. 소거는 곧 행렬을 두 삼각 행렬의 곱으로 쪼개는 일이었습니다.

>>> A = np.array([[4, 2, 1],
...               [2, 5, 2],
...               [2, 3, 5]], dtype=float)
>>> def lu_nopivot(A):
...     n = len(A)
...     U = A.astype(float)
...     L = np.eye(n)
...     for k in range(n - 1):
...         for i in range(k + 1, n):
...             L[i, k] = U[i, k] / U[k, k]
...             U[i, k:] -= L[i, k] * U[k, k:]
...     return L, U
>>> L, U = lu_nopivot(A)
>>> L
array([[1. , 0. , 0. ],
       [0.5, 1. , 0. ],
       [0.5, 0.5, 1. ]])
>>> U
array([[4.  , 2.  , 1.  ],
       [0.  , 4.  , 1.5 ],
       [0.  , 0.  , 3.75]])
>>> np.allclose(L @ U, A)
True

7교시의 naive_gauss와 다른 점은 곱수를 L[i, k]에 적어 두는 것, 그리고 $\mathbf{b}$를 건드리지 않는 것뿐입니다.

이제 $A\mathbf{x} = \mathbf{b}$는 $L(U\mathbf{x}) = \mathbf{b}$입니다. $U\mathbf{x}$를 $\mathbf{y}$라 부르면 두 번에 나눠 풀 수 있습니다.

1. 전진 대입 $L\mathbf{y} = \mathbf{b}$: 위에서 아래로. $y_i = b_i - \sum_{j<i} l_{ij} y_j$
2. 후진 대입 $U\mathbf{x} = \mathbf{y}$: 아래에서 위로. 7교시와 같다.

둘 다 삼각 행렬이라 각각 $n^2$번 정도면 끝납니다. 비싼 $2n^3/3$은 분해할 때 한 번만 냅니다.

>>> def lu_solve(L, U, b):
...     n = len(b)
...     y = np.zeros(n)
...     for i in range(n):
...         y[i] = b[i] - L[i, :i] @ y[:i]
...     x = np.zeros(n)
...     for i in range(n - 1, -1, -1):
...         x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
...     return x
>>> lu_solve(L, U, np.array([44.0, 60, 56]))
array([6., 8., 4.])
>>> demands = np.array([[44, 60, 56],
...                     [50, 60, 56],
...                     [40, 64, 60]], dtype=float)
>>> for d in demands:
...     print(d, "->", lu_solve(L, U, d))
[44. 60. 56.] -> [6. 8. 4.]
[50. 60. 56.] -> [7.9 7.4 3.6]
[40. 64. 60.] -> [4.2 9.2 4.8]

충전소 40곳($n = 40$)을 하루 24번 푼다면 얼마나 아낄까요?

약 9배입니다. $n$이 클수록 차이는 24배에 가까워집니다.

실제로는 7교시의 부분 피벗팅이 끼어듭니다. 행을 바꾼 기록을 치환 행렬(permutation matrix) $P$에 적어 $PA = LU$로 씁니다. $P$는 단위행렬의 행 순서를 바꾼 것이고, 풀 때 $\mathbf{b}$의 순서를 똑같이 바꿔 주기만 하면 됩니다. NumPy에는 LU 분해 함수가 없습니다. SciPy의 이름은 '실무에서는'에서 봅니다.

8-2. 조금씩 고쳐 가기: 야코비와 가우스-자이델 (10분)

전력 팀의 규칙을 다시 보세요. $x_i = (d_i + x_{i-1} + x_{i+1})/3$. 이웃 값을 알면 내 값이 정해지는 꼴입니다. 그렇다면 아무 값에서나 출발해 이 식으로 계속 고쳐 쓰면 어떨까요? 4교시와 5교시에서 근을 찾을 때 했던 "한 걸음씩 다가가기"를 연립방정식에 쓰는 것입니다.

야코비 방법(Jacobi method). 지난 회차의 값만 써서 모두를 한꺼번에 고칩니다. 일반적인 $A\mathbf{x} = \mathbf{b}$에서는 $i$번째 식을 $x_i$에 대해 풀어 씁니다.

가우스-자이델 방법(Gauss–Seidel method). 같은 식인데, 이번 회차에 이미 고친 값이 있으면 바로 씁니다. $x_1$을 고쳤으면 $x_2$를 고칠 때 새 $x_1$을 씁니다.

충전소 4곳, 수요 $\mathbf{d} = (20, 22, 16, 16)$ kW로 해 봅니다. 0에서 출발하면 야코비 1회차는 $x_1 = (20 + 0)/3 = 6.667$, $x_2 = (22 + 0 + 0)/3 = 7.333$입니다. 가우스-자이델은 $x_2$를 고칠 때 방금 구한 6.667을 써서 $x_2 = (22 + 6.667 + 0)/3 = 9.556$이 됩니다.

>>> A4 = np.array([[3, -1, 0, 0],
...                [-1, 3, -1, 0],
...                [0, -1, 3, -1],
...                [0, 0, -1, 3]], dtype=float)
>>> d4 = np.array([20.0, 22, 16, 16])
>>> np.linalg.solve(A4, d4)
array([12., 16., 14., 10.])
>>> def jacobi(A, b, x0, iters):
...     D = np.diag(A)
...     R = A - np.diag(D)
...     x = x0.astype(float)
...     for _ in range(iters):
...         x = (b - R @ x) / D
...     return x
>>> def gauss_seidel(A, b, x0, iters):
...     x = x0.astype(float)
...     for _ in range(iters):
...         for i in range(len(b)):
...             s = A[i] @ x - A[i, i] * x[i]
...             x[i] = (b[i] - s) / A[i, i]
...     return x
>>> x0 = np.zeros(4)
>>> jacobi(A4, d4, x0, 1)
array([6.66666667, 7.33333333, 5.33333333, 5.33333333])
>>> gauss_seidel(A4, d4, x0, 1)
array([6.66666667, 9.55555556, 8.51851852, 8.17283951])

np.diag(A)는 대각 원소를 뽑고, np.diag(D)는 거꾸로 그 값을 대각선에 놓은 행렬을 만듭니다. 그래서 R은 $A$에서 대각선만 0으로 지운 행렬입니다. 야코비는 x 전체를 새 배열로 바꾸고, 가우스-자이델은 x[i]를 그 자리에서 고쳐 다음 줄이 바로 쓰게 합니다. 차이는 그것뿐입니다.

몇 번 고치면 될까요? 반복 횟수마다 잔차 $\|\mathbf{b} - A\mathbf{x}\|$의 가장 큰 성분을 봅니다.

>>> for k in [1, 2, 4, 8, 16]:
...     rj = np.abs(d4 - A4 @ jacobi(A4, d4, x0, k)).max()
...     rg = np.abs(d4 - A4 @ gauss_seidel(A4, d4, x0, k)).max()
...     print(k, f"{rj:.1e}", f"{rg:.1e}")
1 1.3e+01 9.6e+00
2 6.7e+00 4.0e+00
4 2.0e+00 5.4e-01
8 1.7e-01 3.9e-03
16 1.2e-03 2.0e-07

둘 다 잔차가 일정한 비율로 줄어듭니다(로그로 보면 직선). 4교시 이분법처럼 1차 수렴입니다. 가우스-자이델은 같은 횟수에서 맞는 자릿수가 약 두 배라, 대개 야코비보다 두 배쯤 빠릅니다. 대신 야코비는 모든 $x_i$를 서로 기다리지 않고 동시에 고칠 수 있어서 여러 코어나 GPU에 나눠 맡기기 좋습니다.

반복법의 진짜 장점은 크기입니다. 한 번 고칠 때 0이 아닌 칸만 계산하면 되니, 충전소마다 이웃이 둘뿐인 이 문제에서는 한 회차가 $n$에 비례합니다. 소거의 $n^3$과 견주어 보세요. 정지 조건은 4교시처럼 정합니다. 잔차나 두 회차의 차이가 허용 오차보다 작아지면 멈춥니다.

8-3. 언제 수렴하나: 대각 우세 (6분)

반복법이 늘 수렴하지는 않습니다. 가장 쓰기 쉬운 판정법은 대각 우세(diagonal dominance) 입니다. 모든 행에서 대각 원소의 절댓값이 나머지 원소 절댓값의 합보다 크면 됩니다.

충전소 행렬은 모든 행에서 $3 > 1 + 1$이라 대각 우세이고, 야코비와 가우스-자이델 모두 반드시 수렴합니다. 뜻으로 읽으면 "각 식에서 자기 미지수의 영향이 이웃들 영향을 다 합친 것보다 크다"입니다. 그래서 이웃 값이 조금 틀려도 내 값은 덜 틀리게 고쳐지고, 고칠 때마다 오차가 줄어듭니다.

순서가 중요하다는 것도 보여 드리겠습니다. 실습대와 같은 2 × 2 연립방정식입니다.

>>> def diag_dominant(A):
...     d = np.abs(np.diag(A))
...     rest = np.abs(A).sum(axis=1) - d
...     return bool(np.all(d > rest))
>>> P = np.array([[3.0, 2], [1, 2]])
>>> q = np.array([7.0, 5])
>>> np.linalg.solve(P, q)
array([1., 2.])
>>> gauss_seidel(P, q, np.array([-3.0, -1]), 10)
array([1.00010161, 1.99994919])
>>> P_swap, q_swap = P[[1, 0]], q[[1, 0]]
>>> P_swap
array([[1., 2.],
       [3., 2.]])
>>> diag_dominant(P), diag_dominant(P_swap)
(True, False)
>>> gauss_seidel(P_swap, q_swap, np.array([-3.0, -1]), 10)
array([ 118099., -177145.])

두 식의 순서만 바꿨습니다. 같은 연립방정식이고 답도 (1, 2) 그대로인데, 가우스-자이델은 십만 단위로 튀어 나갑니다. 2 × 2에서는 가우스-자이델이 한 번 돌 때 오차가 $\rho = \left|\dfrac{a_{12}a_{21}}{a_{11}a_{22}}\right|$배가 됩니다. 처음 순서는 $\rho = \frac{2 \cdot 1}{3 \cdot 2} = \frac{1}{3}$, 바꾼 순서는 $\rho = \frac{2 \cdot 3}{1 \cdot 2} = 3$입니다. 3배씩 10번이면 $3^{10} \approx 59{,}000$배입니다.

대각 우세는 충분조건입니다. 대각 우세가 아니어도 수렴하는 행렬이 있습니다. 7교시의 허브 행렬은 3행이 $5 = 2 + 3$이라 엄밀하게는 대각 우세가 아니지만, 가우스-자이델이 잘 수렴합니다. 반대로 대각 우세가 아니면 수렴을 장담할 수 없으니, 행 순서를 바꿔 큰 수를 대각선에 모으는 것이 먼저입니다. 7교시 부분 피벗팅과 같은 생각입니다.

아래 실습대의 두 직선이 연립방정식의 두 식입니다. 청록 직선이 식 ①(x를 고칠 때 쓰는 식), 진회색 직선이 식 ②(y를 고칠 때 쓰는 식)이고, 회색 점 (1, 2)가 교점, 곧 참값입니다. 검은 손잡이로 시작점을 끌어 옮길 수 있습니다(키보드 화살표도 됩니다). 가우스-자이델은 시작점에서 가로로 움직여 직선 ①에 닿고(y는 그대로 두고 x를 고침), 거기서 세로로 움직여 직선 ②에 닿습니다(새 x로 y를 고침). 그래서 남색 길이 계단 모양으로 교점을 향합니다. 호박색이 지금 단계입니다. 야코비로 바꾸면 이전 점에서 두 직선까지 가로·세로로 잰 값을 한꺼번에 써서 대각선으로 건너뜁니다. 슬라이더는 식 ①의 계수 a₁₁, a₁₂를 바꿉니다. 직선 ①은 늘 참값 (1, 2)를 지나도록 b₁이 따라 바뀝니다.

8-4. 0이 대부분인 행렬: 희소 행렬 (6분)

충전소 40곳의 행렬을 실제로 만들어 봅시다. np.eye(n, k=1)은 대각선 바로 위 칸만 1인 행렬입니다.

>>> n = 40
>>> T = 3 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
>>> T[:5, :5]
array([[ 3., -1.,  0.,  0.,  0.],
       [-1.,  3., -1.,  0.,  0.],
       [ 0., -1.,  3., -1.,  0.],
       [ 0.,  0., -1.,  3., -1.],
       [ 0.,  0.,  0., -1.,  3.]])
>>> np.count_nonzero(T), T.size
(np.int64(118), 1600)
>>> T.nbytes
12800

1,600칸 중 118칸만 0이 아닙니다. 이렇게 대부분이 0인 행렬을 희소 행렬(sparse matrix) 이라고 합니다. 충전소마다 이웃이 둘뿐이니 한 행에 0이 아닌 수가 많아야 3개입니다. 도로망, 전력망, 소셜 그래프, 미분방정식을 격자로 푼 행렬(12교시)이 모두 이런 모양입니다.

40곳일 때는 12.8 KB라 상관없지만, 충전소가 10만 곳이면 이야기가 다릅니다.

>>> 100_000**2 * 8 / 1e9
80.0
>>> 3 * 100_000 * 8 / 1e6
2.4

모든 칸을 float64(8바이트)로 저장하면 80 GB, 대각선 세 줄만 저장하면 2.4 MB입니다. 소거는 더 나쁩니다. 연산이 $n^3$이고, 소거하는 동안 0이던 칸이 0이 아닌 수로 채워지기도 합니다(fill-in).

반복법은 행렬을 만들 필요조차 없습니다. "이웃의 값을 더해 3으로 나누기"만 할 수 있으면 됩니다. 아래 jacobi_chain은 T를 쓰지 않고 이웃 값을 한 칸씩 밀어서 가져옵니다. np.r_[0, x[:-1]]은 왼쪽 이웃(맨 앞은 0), np.r_[x[1:], 0]은 오른쪽 이웃입니다.

>>> def jacobi_chain(d, iters):
...     x = np.zeros(len(d))
...     for _ in range(iters):
...         left = np.r_[0, x[:-1]]
...         right = np.r_[x[1:], 0]
...         x = (d + left + right) / 3
...     return x
>>> rng = np.random.default_rng(8)
>>> d = rng.uniform(5, 30, size=n).round(1)
>>> d[:6]
array([13.2, 29.7, 13. , 24.7, 26.7, 14.8])
>>> x = jacobi_chain(d, 60)
>>> np.round(x[:6], 2)
array([10.92, 19.55, 18.03, 21.55, 21.92, 17.5 ])
>>> np.abs(d - T @ x).max() < 1e-6
np.True_
>>> np.allclose(x, np.linalg.solve(T, d))
True

한 회차에 덧셈 몇 번과 나눗셈 한 번이 충전소 수만큼, 60회차면 충분합니다. 충전소가 10만 곳이어도 메모리는 배열 몇 개, 계산은 수백만 번입니다. 소거라면 $\frac{2}{3}(10^5)^3 \approx 7 \times 10^{14}$번입니다.

8-5. 답을 얼마나 믿을 수 있나: 조건수 (10분)

1교시에 오차의 출처로 데이터 오차를 적고 "8교시(조건수)"라고 미뤄 두었습니다. 수요 $\mathbf{b}$는 계량기로 잰 값이라 늘 조금씩 틀립니다. 그 작은 틀림이 답을 얼마나 흔들까요?

두 직선이 거의 평행한 연립방정식을 봅시다. 두 식의 계수가 0.0001만 다릅니다.

답은 $(1, 2)$입니다. 이제 첫 식의 오른쪽 5를 0.01% 바꿔 5.0005로 해 봅니다.

>>> C = np.array([[1.0, 2.0], [1.0001, 2.0]])
>>> c = C @ np.array([1.0, 2.0])
>>> c
array([5.    , 5.0001])
>>> np.linalg.solve(C, c)
array([1., 2.])
>>> c_bumped = c * np.array([1.0001, 1.0])
>>> c_bumped
array([5.0005, 5.0001])
>>> np.round(np.linalg.solve(C, c_bumped), 4)
array([-4.    ,  4.5002])

입력을 만분의 일 바꿨더니 답이 $(1, 2)$에서 $(-4, 4.5)$로 갔습니다. 계산이 틀린 것이 아닙니다. 두 직선이 거의 평행하면, 한 직선을 아주 조금만 옮겨도 교점이 직선을 따라 멀리 미끄러집니다. 문제 자체가 예민한 것입니다. 이런 문제를 나쁜 조건(ill-conditioned) 이라고 합니다.

얼마나 예민한지를 숫자 하나로 나타낸 것이 조건수(condition number) $\kappa(A)$입니다. 입력의 상대 오차가 답에서 최대 몇 배로 커질 수 있는지를 말합니다.

$\|\cdot\|$는 벡터의 길이(노름)입니다. 행렬의 노름은 "벡터를 최대 몇 배로 늘이는가"입니다. 그래서 $\kappa(A)$는 $A$가 가장 많이 늘이는 배율을 가장 적게 늘이는 배율로 나눈 값입니다. 선형대수 강의 14교시 SVD의 말로 하면 가장 큰 특이값 ÷ 가장 작은 특이값입니다. 가장 작은 쪽이 0이면(특이 행렬) 조건수는 무한대입니다.

>>> f"{np.linalg.cond(C):.1e}"
'5.0e+04'
>>> round(np.linalg.cond(P), 2)
np.float64(4.27)
>>> round(np.linalg.cond(A), 2), round(np.linalg.cond(A4), 2), round(np.linalg.cond(T), 2)
(np.float64(3.87), np.float64(3.34), np.float64(4.96))

C는 5만이고, 허브 행렬·충전소 행렬·실습대의 기본 행렬은 모두 5 아래입니다. 실제로 C에서 입력의 상대 변화는 약 $7 \times 10^{-5}$, 답의 상대 변화는 2.5(250%)였으니 약 35,000배로 커졌습니다. 조건수 5만 안쪽입니다.

여기서 7교시에 미뤄 둔 이야기가 나옵니다. 흔들린 답 $(-4, 4.5)$를 원래 식에 넣어 보면 잔차가 아주 작습니다.

>>> x_b = np.linalg.solve(C, c_bumped)
>>> np.round(np.abs(C @ x_b - c).max(), 6)
np.float64(0.0005)

잔차 0.0005인데 답은 250% 틀렸습니다. 조건수가 크면 잔차가 작아도 답이 맞다는 보장이 없습니다. 잔차 확인과 조건수 확인은 짝으로 합니다.

자릿수로 읽는 법. 2교시에서 float64의 기계 엡실론이 $\varepsilon \approx 2.2 \times 10^{-16}$, 곧 유효숫자 약 16자리라고 했습니다. 입력을 저장하는 순간 생기는 반올림 오차만으로도 답의 상대 오차는 $\kappa \cdot \varepsilon$까지 커질 수 있습니다. 그래서 실무의 어림은 이렇습니다.

악명 높은 예가 힐베르트 행렬(Hilbert matrix) $H_{ij} = 1/(i + j + 1)$입니다(파이썬처럼 0부터 셈). 생긴 것은 얌전한데 크기를 키우면 조건수가 폭발합니다. 답이 모두 1이 되게 $\mathbf{b}$를 만들고 풀어 봅니다.

>>> H4 = 1 / (np.arange(4)[:, None] + np.arange(4) + 1)
>>> np.round(H4, 3)
array([[1.   , 0.5  , 0.333, 0.25 ],
       [0.5  , 0.333, 0.25 , 0.2  ],
       [0.333, 0.25 , 0.2  , 0.167],
       [0.25 , 0.2  , 0.167, 0.143]])
>>> for n in [4, 8, 12]:
...     H = 1 / (np.arange(n)[:, None] + np.arange(n) + 1)
...     x = np.linalg.solve(H, H @ np.ones(n))
...     print(n, f"{np.linalg.cond(H):.1e}", f"{np.abs(x - 1).max():.1e}")
4 1.6e+04 6.6e-14
8 1.5e+10 1.2e-07
12 1.8e+16 7.0e-01

12 × 12에서는 맞는 자릿수가 하나도 남지 않습니다. np.linalg.solve는 오류도 없이 답을 냈고, 잔차도 작습니다. 7교시의 부분 피벗팅은 계산 방법 때문에 생기는 오차를 막는 장치이고, 조건수는 문제 자체의 예민함입니다. 피벗팅으로는 조건수를 고칠 수 없습니다. 실제 오차가 어림보다 조금 작게 나온 것은 $\kappa \cdot \varepsilon$이 "최악이면 이만큼"이라는 위쪽 한계이기 때문입니다.

실습대에서 직접 확인해 보세요. 오른쪽의 큰 숫자가 지금 행렬의 조건수입니다. 30을 넘으면 호박색, 1000을 넘으면 붉은색으로 바뀝니다. 'b₁을 1% 바꿔 보기'를 켜면 직선 ①이 b₁ +1%만큼 옮겨진 붉은 점선과, 그때의 새 교점(붉은 점)이 나타납니다.

책의 기호 읽는 법

실무에서는

import scipy.linalg
lu, piv = scipy.linalg.lu_factor(A)      # PA = LU 를 한 번만 (부분 피벗팅 포함)
for b in hourly_demands:                 # 매시간 대입만
    x = scipy.linalg.lu_solve((lu, piv), b)

import scipy.sparse, scipy.sparse.linalg
T = scipy.sparse.diags([-1, 3, -1], [-1, 0, 1], shape=(n, n), format="csr")  # 0이 아닌 칸만 저장
x = scipy.sparse.linalg.spsolve(T, d)    # 희소 행렬 전용 소거
x, info = scipy.sparse.linalg.cg(T, d)   # 켤레 기울기법: 대칭·양의 정부호일 때 쓰는 반복법
scipy.linalg.solve_banded((1, 1), ab, d) # 띠 행렬 전용 (삼중대각이면 O(n))

np.linalg.cond(A)                        # 조건수 (SVD로 계산해서 큰 행렬에는 비싸다)

• 같은 행렬로 여러 번 풀면 lu_factor + lu_solve. 오른쪽 변을 한꺼번에 모을 수 있으면 np.linalg.solve(A, B)도 내부에서 한 번만 분해합니다.
• 큰 희소 행렬은 scipy.sparse로 저장하고, 직접법(spsolve)이나 반복법(cg, gmres)으로 풉니다. 실무의 반복법은 야코비·가우스-자이델보다 훨씬 빠른 크릴로프(Krylov) 방법들이지만, "한 번에 조금씩 고친다"는 생각은 같습니다.
• 조건수가 크게 나오면 계산을 의심하기 전에 모델을 의심하세요. 거의 같은 식이 두 개 들어 있지 않은지, 단위가 너무 다른 변수가 섞여 있지 않은지(예: kW와 W) 봅니다.

자주 하는 실수 (3분)

LU 분해를 해 놓고 매번 solve를 부르기. np.linalg.solve(A, b)는 부를 때마다 처음부터 분해합니다. 같은 행렬이면 분해 결과를 저장해 두거나, 오른쪽 변을 열로 모아 한 번에 부르세요.

반복 횟수를 고정해 두고 끝내기. for _ in range(10)처럼 횟수만 정하면 발산해도, 아직 덜 수렴했어도 모릅니다. 잔차나 두 회차의 차이로 멈추고, 최대 횟수를 넘으면 경고를 내세요.

식의 순서를 아무렇게나 두기. 같은 연립방정식도 순서에 따라 가우스-자이델이 수렴하기도, 발산하기도 합니다. 반복법을 쓰기 전에 큰 수가 대각선에 오도록 행을 정렬하세요.

잔차가 작으니 답이 맞다고 믿기. 조건수가 크면 잔차가 작아도 답이 크게 틀릴 수 있습니다. 답에서 믿을 수 있는 자릿수는 16 − log₁₀ κ 정도입니다.

연습 문제 (5분)

1. 충전소 4곳의 행렬 A4를 lu_nopivot으로 분해해 보세요. L과 U에도 0이 많이 남나요? 이것이 큰 삼중대각 행렬에 어떤 의미가 있을까요?
2. A4를 0에서 출발한 가우스-자이델로 풀 때, 잔차의 가장 큰 성분이 10⁻⁶ 아래가 되려면 몇 번 고쳐야 하나요?
3. 6 × 6 힐베르트 행렬의 조건수를 구하고, 답에서 몇 자리가 남을지 어림한 뒤 실제 오차와 견주어 보세요.
4. Q = [[1, 3], [2, 1]], 답 (1, 2)인 연립방정식을 가우스-자이델로 10번 풀어 보세요. 행 순서를 바꾸면 어떻게 될까요?

풀이

>>> L4, U4 = lu_nopivot(A4)
>>> np.round(L4, 3)
array([[ 1.   ,  0.   ,  0.   ,  0.   ],
       [-0.333,  1.   ,  0.   ,  0.   ],
       [ 0.   , -0.375,  1.   ,  0.   ],
       [ 0.   ,  0.   , -0.381,  1.   ]])
>>> np.round(U4, 3)
array([[ 3.   , -1.   ,  0.   ,  0.   ],
       [ 0.   ,  2.667, -1.   ,  0.   ],
       [ 0.   ,  0.   ,  2.625, -1.   ],
       [ 0.   ,  0.   ,  0.   ,  2.619]])
>>> x = np.zeros(4)
>>> k = 0
>>> while np.abs(d4 - A4 @ x).max() > 1e-6:
...     x = gauss_seidel(A4, d4, x, 1)
...     k += 1
>>> k
15
>>> H6 = 1 / (np.arange(6)[:, None] + np.arange(6) + 1)
>>> f"{np.linalg.cond(H6):.1e}"
'1.5e+07'
>>> x6 = np.linalg.solve(H6, H6 @ np.ones(6))
>>> f"{np.abs(x6 - 1).max():.1e}"
'2.4e-10'
>>> Q = np.array([[1.0, 3], [2, 1]])
>>> q = Q @ np.array([1.0, 2])
>>> gauss_seidel(Q, q, np.zeros(2), 10)
array([ 6.0466177e+07, -1.2093235e+08])
>>> gauss_seidel(Q[[1, 0]], q[[1, 0]], np.zeros(2), 10)
array([1.0000001 , 1.99999997])

1번은 L에 대각선 바로 아래, U에 대각선과 바로 위만 남습니다. 삼중대각 행렬은 소거해도 0이 채워지지 않아서, 0이 아닌 칸만 다루면 LU 분해와 대입이 모두 $n$에 비례합니다. 이것을 토마스 알고리즘(Thomas algorithm)이라고 부르고, SciPy의 solve_banded가 하는 일입니다.

2번은 15번입니다. 8-2의 표에서 16번에 잔차가 2.0e-07이었던 것과 맞습니다. 이런 식으로 "잔차가 기준 아래가 될 때까지"를 정지 조건으로 씁니다.

3번은 조건수가 약 1.5 × 10⁷이니 약 7자리를 잃고 9자리쯤 남는다고 어림합니다. 실제 오차는 2.4 × 10⁻¹⁰으로 9~10자리가 맞습니다.

4번은 Q가 대각 우세가 아니라(1행 |1| < |3|) 1억 단위로 발산합니다. 행을 바꾸면 [[2, 1], [1, 3]]이 되어 두 행 모두 대각 우세이고, 10번 만에 (1, 2)에 거의 붙습니다. 2 × 2 비율로는 ρ = 6에서 1/6이 됩니다.

그날 오후, 정비 팀이 관제실에 노트 한 권을 놓고 갑니다. 3호기 배터리의 방전 기록입니다. 1분마다 잰 남은 전력이 점으로 찍혀 있는데, 점들이 한 줄로 늘어서지 않고 들쭉날쭉합니다. "이 기록으로 남은 비행 시간을 미리 알 수 있을까요?" 점은 수십 개인데 직선을 정하는 미지수는 두 개뿐입니다. 식이 미지수보다 많은 연립방정식, 다음 교시의 최소제곱 회귀입니다.