Lecture 15. 부분공간으로의 투영 — 파이썬 실습 - 선형대수

L15 서술 파트에서 얻은 결론은 두 줄이었다. 오차가 수직이면 그 점이 가장 가깝고, 그 조건을 옮겨 적으면 정규방정식 $A^{\mathsf{T}}A\hat{\vv{x}} = A^{\mathsf{T}}\vv{b}$ 가 된다.

이 노트북에서는 그 두 줄을 코드로 확인한다. 특히 "가장 가깝다"는 주장은 말로만 하지 말고 직선 위의 점을 촘촘히 훑어 정말 최소인지 재 보도록 하자.

서술 파트의 내용	여기서 확인하는 방법
$\hat{x} = (\vv{a}^{\mathsf{T}}\vv{b})/(\vv{a}^{\mathsf{T}}\vv{a})$	직선 위를 훑어 최소점과 대조
정규방정식	`solve(A.T @ A, A.T @ b)` 와 오차의 수직성
$P^2 = P$ , $P^{\mathsf{T}} = P$	행렬을 직접 곱해 본다
두 극단	$\vv{b}$ 를 열공간 / 좌영공간에서 골라 넣는다
$\vv{b} = P\vv{b} + (I-P)\vv{b}$	L14의 분해를 이번에는 공식으로

0. 준비¶

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

from linalg_viz import (COLORS, arrow, layout3d, plane, show_matrix,
                        spin_figure)

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

numpy 2.5.2

1. 직선 위로의 투영¶

서술 파트의 예 $\vv{a} = (1, 2)$ , $\vv{b} = (4, 3)$ 으로 시작한다.

def 직선투영(a, b):
    """a 방향 직선 위로 b 를 투영한다. 계수, 투영점, 오차를 돌려준다."""
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    xhat = (a @ b) / (a @ a)
    p = xhat * a
    return xhat, p, b - p

a = np.array([1.0, 2.0])
b = np.array([4.0, 3.0])
xhat, p, e = 직선투영(a, b)

print(f"xhat = (a . b) / (a . a) = {a @ b:g} / {a @ a:g} = {xhat:g}")
print("p =", p)
print("e =", e)
print()
print("a . e =", a @ e, "  -> 수직이다")

xhat = (a . b) / (a . a) = 10 / 5 = 2
p = [2. 4.]
e = [ 2. -1.]

a . e = 0.0   -> 수직이다

투영행렬은 바깥곱을 수로 나눈 것이다. 위와 아래의 순서가 다르다는 점에 주의하자.

P1 = np.outer(a, a) / (a @ a)          # a a^T / (a^T a)

print(show_matrix(P1, "P = a a^T / (a^T a)"))
print("a a^T 의 크기 :", np.outer(a, a).shape, "  (행렬)")
print("a^T a 의 값   :", a @ a, "  (수)")
print()
print("P b =", P1 @ b, "  ->  p 와 같은가 :", np.allclose(P1 @ b, p))
print("rank(P) :", np.linalg.matrix_rank(P1), "  (랭크 1 이어야 한다)")

P = a a^T / (a^T a)
[  0.2   0.4 ]
[  0.4   0.8 ]
a a^T 의 크기 : (2, 2)   (행렬)
a^T a 의 값   : 5.0   (수)

P b = [2. 4.]   ->  p 와 같은가 : True
rank(P) : 1   (랭크 1 이어야 한다)

2. 정말 가장 가까운가¶

서술 파트에서 피타고라스로 증명했지만, 직선 위를 촘촘히 훑어 직접 재 보자.

t = np.linspace(-1.0, 5.0, 601)
거리제곱 = np.array([(b - s * a) @ (b - s * a) for s in t])

최소자리 = int(np.argmin(거리제곱))
print(f"훑어서 찾은 최소  : t = {t[최소자리]:.3f},  |b - t a|^2 = {거리제곱[최소자리]:.4f}")
print(f"공식이 준 값      : t = {xhat:.3f},  |e|^2 = {e @ e:.4f}")
print("같은가 :", np.isclose(t[최소자리], xhat, atol=1e-2))

훑어서 찾은 최소  : t = 2.000,  |b - t a|^2 = 5.0000
공식이 준 값      : t = 2.000,  |e|^2 = 5.0000
같은가 : True

거리의 제곱은 $t$ 에 대한 이차식이므로 포물선이다. 그 꼭짓점이 $\hat{x}$ 이다.

곡선 = go.Scatter(x=t, y=거리제곱, mode="lines", name="|b - t a|^2",
                 line=dict(color=COLORS["input"], width=3))
꼭짓 = go.Scatter(x=[xhat], y=[e @ e], mode="markers+text", name="xhat",
                 marker=dict(color=COLORS["error"], size=12),
                 text=[f"  xhat = {xhat:g}"], textposition="middle right")

go.Figure(data=[곡선, 꼭짓], layout=go.Layout(
    title=dict(text="직선 위를 훑은 거리의 제곱"),
    xaxis=dict(title=dict(text="t")),
    yaxis=dict(title=dict(text="|b - t a|^2")),
    height=380, margin=dict(l=60, r=20, t=50, b=50)))

$\vv{q} = \vv{a}$ 로 잡았을 때의 피타고라스도 확인해 보자.

q = a.copy()
print(f"|b - p|^2 = {(b - p) @ (b - p):g}")
print(f"|p - q|^2 = {(p - q) @ (p - q):g}")
print(f"합        = {(b - p) @ (b - p) + (p - q) @ (p - q):g}")
print(f"|b - q|^2 = {(b - q) @ (b - q):g}")
print("맞는가 :", np.isclose((b - p) @ (b - p) + (p - q) @ (p - q), (b - q) @ (b - q)))

|b - p|^2 = 5
|p - q|^2 = 5
합        = 10
|b - q|^2 = 10
맞는가 : True

3. 정규방정식¶

이제 직선이 아니라 평면으로 떨어뜨린다. 조건은 그대로 "오차가 수직일 것"이다.

def 투영(A, b):
    """C(A) 위로 b 를 투영한다. A 의 열은 독립이어야 한다."""
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float)
    if np.linalg.matrix_rank(A) < A.shape[1]:
        raise ValueError("A 의 열이 종속이라 A^T A 가 가역이 아니다")
    xhat = np.linalg.solve(A.T @ A, A.T @ b)      # 정규방정식
    p = A @ xhat
    return xhat, p, b - p

A = np.array([[1.0, 0.0],
              [0.0, 1.0],
              [1.0, 1.0]])
b3 = np.array([3.0, 2.0, 2.0])

print(show_matrix(A.T @ A, "A^T A"))
print("A^T b =", A.T @ b3)

xhat3, p3, e3 = 투영(A, b3)
print()
print("xhat =", xhat3)
print("p    =", p3)
print("e    =", e3)

A^T A
[  2   1 ]
[  1   2 ]
A^T b = [5. 4.]

xhat = [2. 1.]
p    = [2. 1. 3.]
e    = [ 1.  1. -1.]

오차가 각 열과 수직인지 확인한다. A.T @ e 가 바로 그 내적표이다.

print("A^T e =", A.T @ e3, "  -> 영벡터여야 한다")
print()
for j in range(A.shape[1]):
    print(f"  열 {j + 1} = {A[:, j]} 와의 내적 : {A[:, j] @ e3:+.1e}")

A^T e = [0. 0.]   -> 영벡터여야 한다

  열 1 = [1. 0. 1.] 와의 내적 : +0.0e+00
  열 2 = [0. 1. 1.] 와의 내적 : +0.0e+00

$\vv{b}$ 를 여러 개 바꿔 넣어도 오차는 언제나 수직이다.

print(f"{'b':>22}{'|e|':>9}{'A^T e 의 최대 절댓값':>22}")
for _ in range(6):
    bb = rng.integers(-5, 6, 3).astype(float)
    _, pp, ee = 투영(A, bb)
    print(f"{str(bb):>22}{np.linalg.norm(ee):>9.3f}{np.abs(A.T @ ee).max():>22.2e}")

                     b      |e|        A^T e 의 최대 절댓값
            [5. 2. 2.]    2.887              4.44e-16
         [ 3. -3. -2.]    1.155              4.44e-16
         [-3. -5. -1.]    4.041              4.44e-16
         [ 1.  5. -4.]    5.774              4.44e-16
            [5. 2. 0.]    4.041              4.44e-16
         [-2.  2.  0.]    0.000              0.00e+00

4. 투영행렬의 성질¶

def 투영행렬(A):
    """C(A) 로의 투영행렬 P = A (A^T A)^-1 A^T 를 만든다."""
    A = np.asarray(A, dtype=float)
    if np.linalg.matrix_rank(A) < A.shape[1]:
        raise ValueError("A 의 열이 종속이라 A^T A 가 가역이 아니다")
    return A @ np.linalg.inv(A.T @ A) @ A.T

P = 투영행렬(A)
print(show_matrix(P, "P = A (A^T A)^-1 A^T"))
print("P b 가 p 와 같은가 :", np.allclose(P @ b3, p3))
print()
print("대칭인가  P^T = P :", np.allclose(P.T, P))
print("멱등인가  P^2 = P :", np.allclose(P @ P, P))
print("rank(P) :", np.linalg.matrix_rank(P), "  (= rank(A) =",
      np.linalg.matrix_rank(A), ")")
print("det(P)  :", np.linalg.det(P), "  -> 가역이 아니다")

P = A (A^T A)^-1 A^T
[   0.667   -0.333    0.333 ]
[  -0.333    0.667    0.333 ]
[   0.333    0.333    0.667 ]
P b 가 p 와 같은가 : True

대칭인가  P^T = P : True
멱등인가  P^2 = P : True
rank(P) : 2   (= rank(A) = 2 )
det(P)  : 0.0   -> 가역이 아니다

서술 파트에서 $P = A(A^{\mathsf{T}}A)^{-1}A^{\mathsf{T}}$ 를 $I$ 로 약분하면 안 된다고 했다. $A$ 가 정사각이 아니면 $A^{-1}$ 이 아예 없기 때문이다. 실제로 시켜 보자.

print("A 의 크기 :", A.shape, " -> 정사각이 아니다")
try:
    np.linalg.inv(A)
except np.linalg.LinAlgError as 오류:
    print("np.linalg.inv(A) ->", type(오류).__name__, ":", 오류)

print()
print("A^T A 의 크기 :", (A.T @ A).shape, " -> 이쪽은 가역이다")
print("det(A^T A) :", np.linalg.det(A.T @ A))

A 의 크기 : (3, 2)  -> 정사각이 아니다
np.linalg.inv(A) -> LinAlgError : Last 2 dimensions of the array must be square

A^T A 의 크기 : (2, 2)  -> 이쪽은 가역이다
det(A^T A) : 2.9999999999999996

$A$ 가 정사각이고 가역이면 정말로 $P = I$ 가 된다. 투영할 것이 없는 경우이다.

정사각 = np.array([[1.0, 2.0], [3.0, 4.0]])           # 가역인 2 x 2
print(show_matrix(정사각, "정사각 A"))
print(show_matrix(투영행렬(정사각), "그때의 P"))
print("I 인가 :", np.allclose(투영행렬(정사각), np.eye(2)))

정사각 A
[  1   2 ]
[  3   4 ]
그때의 P
[          1   -1.78e-15 ]
[          0           1 ]
I 인가 : True

5. 두 극단과 $I - P$ ¶

$\vv{b}$ 가 이미 열공간 안이면 $\vv{p} = \vv{b}$ , 좌영공간에 있으면 $\vv{p} = \vv{0}$ 이다.

b_안 = A @ np.array([2.0, -1.0])         # 열의 결합이므로 반드시 열공간 안
b_밖 = null_space(A.T)[:, 0]             # 좌영공간의 벡터

for 이름, bb in [("b in C(A)", b_안), ("b in N(A^T)", b_밖)]:
    _, pp, ee = 투영(A, bb)
    print(f"{이름}")
    print(f"  b = {np.round(bb, 6) + 0.0}")
    print(f"  p = {np.round(pp, 6) + 0.0}")
    print(f"  e = {np.round(ee, 6) + 0.0}")
    print(f"  p = b 인가 : {np.allclose(pp, bb)},   p = 0 인가 : {np.allclose(pp, 0)}")
    print()

b in C(A)
  b = [ 2. -1.  1.]
  p = [ 2. -1.  1.]
  e = [0. 0. 0.]
  p = b 인가 : True,   p = 0 인가 : False

b in N(A^T)
  b = [-0.577 -0.577  0.577]
  p = [0. 0. 0.]
  e = [-0.577 -0.577  0.577]
  p = b 인가 : False,   p = 0 인가 : True

$I - P$ 도 대칭이고 멱등이다. 이쪽은 좌영공간으로 투영한다.

Q = np.eye(3) - P

print(show_matrix(Q, "I - P"))
print("대칭인가 :", np.allclose(Q.T, Q))
print("멱등인가 :", np.allclose(Q @ Q, Q))
print("rank     :", np.linalg.matrix_rank(Q), "  (= m - r = 3 - 2 =", 3 - 2, ")")
print()
print("(I - P) b = ", Q @ b3, "  ->  e 와 같은가 :", np.allclose(Q @ b3, e3))
print("A^T (I-P) b =", A.T @ (Q @ b3), "  -> 좌영공간에 들어간다")

I - P
[   0.333    0.333   -0.333 ]
[   0.333    0.333   -0.333 ]
[  -0.333   -0.333    0.333 ]
대칭인가 : True
멱등인가 : True
rank     : 1   (= m - r = 3 - 2 = 1 )

(I - P) b =  [ 1.  1. -1.]   ->  e 와 같은가 : True
A^T (I-P) b = [0. 0.]   -> 좌영공간에 들어간다

L14에서 존재와 유일성만 보였던 분해를 이제 공식으로 구할 수 있다.

print("b       =", b3)
print("P b     =", P @ b3, "  (열공간 성분)")
print("(I-P) b =", Q @ b3, "  (좌영공간 성분)")
print()
print("합이 b 인가   :", np.allclose(P @ b3 + Q @ b3, b3))
print("서로 직교한가 :", np.isclose((P @ b3) @ (Q @ b3), 0, atol=1e-12))
print("피타고라스    :",
      np.isclose(b3 @ b3, (P @ b3) @ (P @ b3) + (Q @ b3) @ (Q @ b3)))

b       = [3. 2. 2.]
P b     = [2. 1. 3.]   (열공간 성분)
(I-P) b = [ 1.  1. -1.]   (좌영공간 성분)

합이 b 인가   : True
서로 직교한가 : True
피타고라스    : True

6. 직각을 여러 각도에서 보기¶

빨간 오차가 청록 평면에 정말 수직인지, 평면을 옆에서 보는 각도로 돌려 확인해 보자.

그림 = [plane(spans=[A[:, 0], A[:, 1]], extent=2.9, color=COLORS["colspace"],
             opacity=0.35, name="C(A)")]
그림 += arrow([0, 0, 0], A[:, 0], COLORS["colspace"], "열 1")
그림 += arrow([0, 0, 0], A[:, 1], COLORS["second"], "열 2")
그림 += arrow([0, 0, 0], b3, COLORS["output"], "b")
그림 += arrow([0, 0, 0], p3, COLORS["input"], "p = Pb")
그림 += arrow(p3, b3, COLORS["error"], "e = b - p", dashed=True)

spin_figure(그림, title="투영 : 열공간 위에서 b 에 가장 가까운 점",
            extent=3.4, n_frames=48, radius=2.6)

$\vv{b}$ 를 바꿔 가며 오차가 언제나 수직인지 보자. 슬라이더 대신 여러 개를 한 번에 그린다.

후보 = [np.array([3.0, 2, 2]), np.array([-2.0, 3, 1]), np.array([1.0, -3, 4])]

여러개 = [plane(spans=[A[:, 0], A[:, 1]], extent=3.4, color=COLORS["colspace"],
               opacity=0.28, name="C(A)")]
for k, bb in enumerate(후보, 1):
    _, pp, ee = 투영(A, bb)
    여러개 += arrow([0, 0, 0], bb, COLORS["output"], f"b{k}", legend=(k == 1))
    여러개 += arrow([0, 0, 0], pp, COLORS["input"], f"p{k}", legend=(k == 1))
    여러개 += arrow(pp, bb, COLORS["error"], f"e{k}", dashed=True, legend=(k == 1))
    print(f"b{k} = {bb}   A^T e = {np.round(A.T @ ee, 12) + 0.0}")

go.Figure(data=여러개, layout=layout3d("b 를 바꿔도 오차는 언제나 수직이다", extent=4.2))

b1 = [3. 2. 2.]   A^T e = [0. 0.]
b2 = [-2.  3.  1.]   A^T e = [0. 0.]
b3 = [ 1. -3.  4.]   A^T e = [0. 0.]

마치며...¶

서술 파트의 내용	이 노트북의 코드
직선 위로의 투영	`직선투영(a, b)`
$P = \vv{a}\vv{a}^{\mathsf{T}} / (\vv{a}^{\mathsf{T}}\vv{a})$	`np.outer(a, a) / (a @ a)`
가장 가깝다	직선 위를 훑어 최소점과 대조
정규방정식	`solve(A.T @ A, A.T @ b)`
오차가 수직	`A.T @ e` 가 영벡터
$P$ 의 성질	`P.T == P`, `P @ P == P`, `rank P = r`
약분 금지	`inv(A)` 는 `LinAlgError`
$I - P$	좌영공간으로의 투영

더 해 볼 것¶

투영 은 열이 종속이면 오류를 낸다. 열이 종속인 $A$ 를 넣어 확인하고, 그래도 투영 자체는 잘 정의되는지 생각해 보자. 무엇이 유일하지 않은가?
투영행렬(A) 와 투영행렬(A @ M) 을 비교해 보자. $M$ 이 가역인 $2 \times 2$ 행렬이면 두 $P$ 가 같은가? 왜 그런가?
$P$ 를 여러 번 곱해 보자. $P^{10}$ 은 무엇인가?
$\vv{b}$ 를 무작위로 100개 뽑아 $\|P\vv{b}\|^2 + \|(I-P)\vv{b}\|^2 = \|\vv{b}\|^2$ 이 언제나 성립하는지 확인해 보자.

다음 강의에서 이 도구가 통계학을 떠받치고 있다는 것을 보게 된다. 표에 찍힌 점들에 직선을 맞추는 일이 정확히 이번 강의의 투영이다.

Lecture 15. 부분공간으로의 투영 — 파이썬 실습