Lecture 16. 투영행렬과 최소제곱 — 파이썬 실습 - 선형대수

L16 서술 파트의 요점은 계산이 아니라 두 그림이 같다는 것이었다. 산점도의 잔차 막대 세 개와 3차원 그림의 오차 벡터 세 성분이 같은 숫자라는 것.

이 노트북에서는 그 두 그림을 나란히 놓고 같은 숫자가 양쪽에 뜨는 것을 확인한다. 그리고 계수를 손으로 흔들어 잔차 제곱합이 정말 커지는지 슬라이더로 재 본다.

서술 파트의 내용	여기서 확인하는 방법
해가 없다	좌영공간 벡터와 $\vv{b}$ 의 내적이 0이 아니다
정규방정식	`solve(A.T @ A, A.T @ b)`
두 개의 그림	2차원과 3차원을 한 그림에 나란히
최소인가	기울기를 흔들며 잔차 제곱합을 잰다
미적분 경로	편미분 두 개가 정규방정식과 같은지 대조
열이 모델의 항	포물선과 다변수로 확장

0. 준비¶

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

from linalg_viz import COLORS, arrow, layout2d, plane, show_matrix, slider_figure

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

numpy 2.5.2

1. 해가 없다는 것부터 확인한다¶

데이터 세 개에 직선 $y = C + Dt$ 를 맞춘다. $A$ 의 첫 열은 상수항, 둘째 열은 $t$ 이다.

t = np.array([1.0, 2.0, 3.0])
y = np.array([1.0, 3.0, 2.0])

A = np.column_stack([np.ones_like(t), t])       # 열 = 모델의 항
b = y.copy()

print(show_matrix(A, "A = [1  t]"))
print("b =", b)
print("크기 :", A.shape, " -> 식 3개, 미지수 2개")

A = [1  t]
[  1   1 ]
[  1   2 ]
[  1   3 ]
b = [1. 3. 2.]
크기 : (3, 2)  -> 식 3개, 미지수 2개

L8의 판정을 그대로 쓴다. 좌영공간의 벡터와 $\vv{b}$ 의 내적이 0이 아니면 해가 없다.

좌영 = null_space(A.T)
print(show_matrix(좌영, "N(A^T) 의 기저"))
print("정수로 보면 :", np.round(좌영[:, 0] / 좌영[0, 0], 6))
print()
print("y^T b =", float(좌영[:, 0] @ b), " -> 0 이 아니므로 해가 없다")
print()
확장 = np.column_stack([A, b])
print("rank(A)     :", np.linalg.matrix_rank(A))
print("rank([A|b]) :", np.linalg.matrix_rank(확장), " -> 늘었으므로 b 는 열공간 밖")

N(A^T) 의 기저
[   0.408 ]
[  -0.816 ]
[   0.408 ]
정수로 보면 : [ 1. -2.  1.]

y^T b = -1.2247448713915885  -> 0 이 아니므로 해가 없다

rank(A)     : 2
rank([A|b]) : 3  -> 늘었으므로 b 는 열공간 밖

2. 정규방정식으로 푼다¶

def 최소제곱(A, b):
    """정규방정식 A^T A xhat = A^T b 를 풀어 계수, 예측값, 잔차를 돌려준다."""
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float)
    xhat = np.linalg.solve(A.T @ A, A.T @ b)
    p = A @ xhat
    return xhat, p, b - p

print(show_matrix(A.T @ A, "A^T A"))
print("A^T b =", A.T @ b)
print()
print("성분의 뜻 :  m =", len(t), ",  sum t =", t.sum(),
      ",  sum t^2 =", (t ** 2).sum())
print("            sum y =", y.sum(), ",  sum t y =", (t * y).sum())

A^T A
[   3    6 ]
[   6   14 ]
A^T b = [ 6. 13.]

성분의 뜻 :  m = 3 ,  sum t = 6.0 ,  sum t^2 = 14.0
            sum y = 6.0 ,  sum t y = 13.0

계수, p, e = 최소제곱(A, b)
C, D = 계수

print(f"맞춘 직선 : y = {C:g} + {D:g} t")
print("p =", p, "  (직선 위의 값)")
print("e =", e, "  (잔차)")
print(f"|e|^2 = {e @ e:g}")
print()
print("A^T e =", A.T @ e, "  -> 수직이다")
print("A xhat 가 b 와 같은가 :", np.allclose(A @ 계수, b), "  <- 해가 아니다")

맞춘 직선 : y = 1 + 0.5 t
p = [1.5 2.  2.5]   (직선 위의 값)
e = [-0.5  1.  -0.5]   (잔차)
|e|^2 = 1.5

A^T e = [0. 0.]   -> 수직이다
A xhat 가 b 와 같은가 : False   <- 해가 아니다

3. 두 개의 그림을 나란히¶

왼쪽은 산점도와 잔차 막대, 오른쪽은 $\R^3$ 의 열공간이다. 빨간 숫자 세 개가 양쪽에 똑같이 나타난다.

그림 = make_subplots(rows=1, cols=2,
                    specs=[[{"type": "xy"}, {"type": "scene"}]],
                    subplot_titles=("데이터 평면 : 잔차 막대",
                                    "열공간 : 오차 벡터"),
                    column_widths=[0.45, 0.55])

# --- 왼쪽 : 데이터 평면 ---
격자t = np.array([0.4, 3.6])
그림.add_trace(go.Scatter(x=격자t, y=C + D * 격자t, mode="lines",
                        line=dict(color=COLORS["input"], width=3),
                        name=f"y = {C:g} + {D:g} t"), row=1, col=1)
막대x, 막대y = [], []
for i in range(3):
    막대x += [t[i], t[i], None]
    막대y += [p[i], y[i], None]
그림.add_trace(go.Scatter(x=막대x, y=막대y, mode="lines",
                        line=dict(color=COLORS["error"], width=5),
                        name="잔차"), row=1, col=1)
그림.add_trace(go.Scatter(x=t, y=y, mode="markers+text",
                        marker=dict(color=COLORS["output"], size=13),
                        text=[f" e={v:+g}" for v in e], textposition="middle right",
                        name="데이터"), row=1, col=1)

# --- 오른쪽 : 열공간 ---
for 트레이스 in ([plane(spans=[A[:, 0], A[:, 1]], extent=3.0,
                      color=COLORS["colspace"], opacity=0.3, name="C(A)")]
               + arrow([0, 0, 0], b, COLORS["output"], "b")
               + arrow([0, 0, 0], p, COLORS["input"], "p")
               + arrow(p, b, COLORS["error"], "e", dashed=True)):
    그림.add_trace(트레이스, row=1, col=2)

축 = dict(range=[-3.2, 3.2], zeroline=True, showbackground=False,
         gridcolor="rgba(0,0,0,0.08)")
그림.update_layout(height=520, margin=dict(l=40, r=10, t=60, b=40),
                 scene=dict(xaxis=dict(title=dict(text="x1"), **축),
                            yaxis=dict(title=dict(text="x2"), **축),
                            zaxis=dict(title=dict(text="x3"), **축),
                            aspectmode="cube",
                            camera=dict(eye=dict(x=-1.7, y=0.6, z=0.9))))
그림.update_xaxes(title_text="t", row=1, col=1)
그림.update_yaxes(title_text="y", row=1, col=1)
그림

왼쪽 막대 옆에 적힌 세 숫자와 오른쪽 빨간 벡터의 성분을 대조해 보자.

print("왼쪽 : 잔차 막대의 길이 (부호 포함)")
for i in range(3):
    print(f"  t = {t[i]:g} 에서  y - (C + D t) = {y[i]:g} - {p[i]:g} = {e[i]:+g}")
print()
print("오른쪽 : 오차 벡터 e =", e)
print("같은 세 숫자인가 :", np.allclose(e, [y[i] - p[i] for i in range(3)]))

왼쪽 : 잔차 막대의 길이 (부호 포함)
  t = 1 에서  y - (C + D t) = 1 - 1.5 = -0.5
  t = 2 에서  y - (C + D t) = 3 - 2 = +1
  t = 3 에서  y - (C + D t) = 2 - 2.5 = -0.5

오른쪽 : 오차 벡터 e = [-0.5  1.  -0.5]
같은 세 숫자인가 : True

4. 기울기를 흔들어 본다¶

기울기 $D$ 를 바꾸면서 그때그때 가장 좋은 $C$ 를 골라도, 잔차 제곱합은 $D = 0.5$ 에서 최소여야 한다.

D후보 = np.linspace(0.0, 1.0, 21)
프레임, 라벨 = [], []

for d in D후보:
    c = (y.sum() - d * t.sum()) / len(t)        # 그 기울기에서 가장 좋은 절편
    예측 = c + d * t
    잔차 = y - 예측
    막대x, 막대y = [], []
    for i in range(3):
        막대x += [t[i], t[i], None]
        막대y += [예측[i], y[i], None]
    프레임.append([
        go.Scatter(x=격자t, y=c + d * 격자t, mode="lines",
                   line=dict(color=COLORS["input"], width=3), name="직선"),
        go.Scatter(x=막대x, y=막대y, mode="lines",
                   line=dict(color=COLORS["error"], width=5), name="잔차"),
        go.Scatter(x=t, y=y, mode="markers",
                   marker=dict(color=COLORS["output"], size=13), name="데이터"),
    ])
    라벨.append(f"D={d:.2f}  |e|^2={잔차 @ 잔차:.3f}")

바탕 = layout2d("기울기를 바꾸면 잔차 제곱합이 어떻게 되는가", extent=4)
바탕["xaxis"] = dict(title=dict(text="t"), range=[0.2, 4.0])
바탕["yaxis"] = dict(title=dict(text="y"), range=[0.2, 3.6],
                   scaleanchor="x", scaleratio=1)

slider_figure(프레임, 라벨, 바탕, prefix="", initial=10)

슬라이더를 움직여 $|e|^2$ 가 가장 작아지는 자리를 찾아보자. 숫자로도 확인한다.

def 잔차제곱합(d):
    c = (y.sum() - d * t.sum()) / len(t)
    잔차 = y - (c + d * t)
    return 잔차 @ 잔차

훑기 = np.linspace(0.0, 1.0, 1001)
값 = np.array([잔차제곱합(d) for d in 훑기])
최소 = 훑기[int(np.argmin(값))]

print(f"훑어서 찾은 최소 : D = {최소:.3f},  |e|^2 = {값.min():.5f}")
print(f"정규방정식의 답  : D = {D:.3f},  |e|^2 = {e @ e:.5f}")
print("같은가 :", np.isclose(최소, D, atol=1e-3))

훑어서 찾은 최소 : D = 0.500,  |e|^2 = 1.50000
정규방정식의 답  : D = 0.500,  |e|^2 = 1.50000
같은가 : True

5. 미적분으로 풀어도 같은 식¶

편미분 $\partial f/\partial x_k$ 가 $A$ 의 $k$ 열과 오차의 내적의 -2 배라고 했다. 수치미분으로 확인해 보자.

def 오차제곱(x):
    """잔차 제곱합 f(x) = |b - A x|^2."""
    잔차 = b - A @ x
    return 잔차 @ 잔차

x0 = np.array([0.4, 0.9])          # 최적점이 아닌 아무 자리
h = 1e-6

수치미분 = np.array([
    (오차제곱(x0 + h * np.eye(2)[k]) - 오차제곱(x0 - h * np.eye(2)[k])) / (2 * h)
    for k in range(2)
])
공식 = -2 * A.T @ (b - A @ x0)

print("수치미분          :", 수치미분)
print("-2 A^T (b - A x)  :", 공식)
print("같은가 :", np.allclose(수치미분, 공식))

수치미분          : [1.2 4. ]
-2 A^T (b - A x)  : [1.2 4. ]
같은가 : True

따라서 편미분을 0으로 놓는 것과 정규방정식은 같은 식이다. 최적점에서는 둘 다 0이 된다.

print("최적점에서의 기울기 :", -2 * A.T @ (b - A @ 계수))
print()
print("편미분을 0 으로 둔 두 식")
print(f"  d/dC = 0 :  {len(t):g} C + {t.sum():g} D = {y.sum():g}")
print(f"  d/dD = 0 :  {t.sum():g} C + {(t ** 2).sum():g} D = {(t * y).sum():g}")
print()
print("정규방정식 A^T A xhat = A^T b")
print(show_matrix(np.column_stack([A.T @ A, A.T @ b]), "[A^T A | A^T b]"))

최적점에서의 기울기 : [0. 0.]

편미분을 0 으로 둔 두 식
  d/dC = 0 :  3 C + 6 D = 6
  d/dD = 0 :  6 C + 14 D = 13

정규방정식 A^T A xhat = A^T b
[A^T A | A^T b]
[   3    6    6 ]
[   6   14   13 ]

6. 열이 곧 모델의 항이다¶

열을 바꾸면 모델이 바뀐다. 정규방정식은 그대로이다. 점을 다섯 개로 늘려서 비교해 보자.

t5 = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y5 = np.array([1.0, 3.0, 2.0, 4.0, 3.0])

모델 = {
    "y = C":             np.column_stack([np.ones_like(t5)]),
    "y = C + Dt":        np.column_stack([np.ones_like(t5), t5]),
    "y = C + Dt + Et^2": np.column_stack([np.ones_like(t5), t5, t5 ** 2]),
}

print(f"{'model':<20}{'columns':>9}{'|e|^2':>10}")
결과 = {}
for 이름, M in 모델.items():
    계수M, pM, eM = 최소제곱(M, y5)
    결과[이름] = (M, 계수M, pM, eM)
    print(f"{이름:<20}{M.shape[1]:>9}{eM @ eM:>10.4f}")

model                 columns     |e|^2
y = C                       1    5.2000
y = C + Dt                  2    2.7000
y = C + Dt + Et^2           3    2.0571

열을 늘릴수록 잔차가 줄어든다. 그러나 열이 점의 개수만큼 되면 잔차가 0이 되는데, 그것은 잘 맞춘 것이 아니라 점을 전부 지나는 곡선을 그린 것이다.

가득 = np.column_stack([t5 ** k for k in range(5)])      # 4차 다항식, 열 5개
계수가득, p가득, e가득 = 최소제곱(가득, y5)

print("열 5개, 점 5개")
print("  rank :", np.linalg.matrix_rank(가득))
print(f"  |e|^2 = {e가득 @ e가득:.2e}   -> 사실상 0")
print("  예측값 :", p가득)
print("  데이터 :", y5)
print()
print("L13 의 표에서 r = m = n 인 경우이다. 오차가 0이라고 좋은 모델은 아니다.")

열 5개, 점 5개
  rank : 5
  |e|^2 = 4.17e-23   -> 사실상 0
  예측값 : [1. 3. 2. 4. 3.]
  데이터 : [1. 3. 2. 4. 3.]

L13 의 표에서 r = m = n 인 경우이다. 오차가 0이라고 좋은 모델은 아니다.

격자 = np.linspace(-0.3, 4.3, 200)
곡선들 = [go.Scatter(x=t5, y=y5, mode="markers", name="데이터",
                   marker=dict(color=COLORS["output"], size=13))]
for (이름, (M, 계수M, _, eM)), 색 in zip(결과.items(),
                                      [COLORS["nullspace"], COLORS["input"],
                                       COLORS["second"]]):
    차수 = M.shape[1]
    값 = sum(계수M[k] * 격자 ** k for k in range(차수))
    곡선들.append(go.Scatter(x=격자, y=값, mode="lines", name=f"{이름}  (열 {차수})",
                          line=dict(color=색, width=3)))
곡선들.append(go.Scatter(
    x=격자, y=sum(계수가득[k] * 격자 ** k for k in range(5)), mode="lines",
    name="4차 (열 5개)", line=dict(color=COLORS["error"], width=2, dash="dash")))

go.Figure(data=곡선들, layout=go.Layout(
    title=dict(text="열을 늘릴수록 잘 맞지만, 끝까지 늘리면 그냥 점을 잇는다"),
    xaxis=dict(title=dict(text="t"), range=[-0.3, 4.3]),
    yaxis=dict(title=dict(text="y"), range=[0.0, 5.0]),
    height=460, margin=dict(l=60, r=20, t=60, b=50)))

7. numpy 의 `lstsq` 와 대조¶

실무에서는 정규방정식을 직접 세우지 않고 lstsq 를 쓴다. 답이 같은지 확인하자.

lstsq답 = np.linalg.lstsq(A, b, rcond=None)[0]

print("정규방정식 :", 계수)
print("lstsq      :", lstsq답)
print("같은가 :", np.allclose(계수, lstsq답))
print()
print("polyfit    :", np.polyfit(t, y, 1)[::-1], "  (차수가 낮은 항부터)")

정규방정식 : [1.  0.5]
lstsq      : [1.  0.5]
같은가 : True

polyfit    : [1.  0.5]   (차수가 낮은 항부터)

답은 같지만 계산 방법이 다르다. lstsq 는 $A^{\mathsf{T}}A$ 를 만들지 않는데, 서술 파트에서 예고한 대로 그쪽이 수치적으로 안전하기 때문이다. 이 이야기는 L33에서 마무리한다.

나쁜t = np.array([1.0, 1.0 + 1e-7, 1.0 + 2e-7])        # t 값이 거의 같다
나쁜A = np.column_stack([np.ones(3), 나쁜t])
나쁜b = np.array([1.0, 2.0, 3.0])

print("A^T A =", (나쁜A.T @ 나쁜A).ravel())
print("det(A^T A) =", np.linalg.det(나쁜A.T @ 나쁜A), "  -> 0 에 가깝다")
print("조건수 =", f"{np.linalg.cond(나쁜A.T @ 나쁜A):.3e}")
print()
print("t 값이 거의 같으면 열이 거의 종속이고, A^T A 를 뒤집기가 위험해진다.")

A^T A = [3. 3. 3. 3.]
det(A^T A) = 5.995204932496273e-14   -> 0 에 가깝다
조건수 = 5.971e+14

t 값이 거의 같으면 열이 거의 종속이고, A^T A 를 뒤집기가 위험해진다.

마치며...¶

서술 파트의 내용	이 노트북의 코드
해가 없다	`null_space(A.T)` 와 `b` 의 내적
정규방정식	`최소제곱(A, b)`
두 개의 그림	`make_subplots` 로 2차원과 3차원을 나란히
최소인가	기울기 슬라이더와 1001점 훑기
미적분 경로	수치미분과 $-2A^{\mathsf{T}}(\vv{b}-A\vv{x})$ 대조
열이 모델의 항	상수 / 직선 / 포물선 / 4차 비교
뒤집기의 위험	$t$ 가 거의 같을 때의 조건수

더 해 볼 것¶

3절의 그림을 마우스로 돌려 빨간 오차 벡터가 청록 평면에 정말 수직인지 확인해 보자. 어느 각도에서 가장 잘 보이는가?
데이터 $y$ 를 바꿔 가며 두 그림이 함께 움직이는지 보자. $y$ 를 직선 위의 값으로 정확히 주면 오차 벡터는 어떻게 되는가?
6절의 표에 "3차 다항식"을 추가해 보자. 잔차는 얼마나 줄어드는가?
4절의 슬라이더에서 절편 $C$ 를 최적값으로 두지 말고 1로 고정해 보자. 최소가 되는 $D$ 가 달라지는가?

다음 강의에서는 $A^{\mathsf{T}}A$ 를 뒤집는 일 자체를 피하는 방법을 배운다. 열들이 애초에 서로 직교하면 $A^{\mathsf{T}}A$ 가 대각행렬이 되기 때문이다.

Lecture 16. 투영행렬과 최소제곱 — 파이썬 실습