Lecture 29. 특이값 분해 (SVD) — 파이썬 실습

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, layout2d, show_matrix, slider_figure

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

def 손으로SVD(A, 눈감아=None):
    """정의대로 SVD 를 구한다. np.linalg.svd 가 속으로 하는 일의 재현."""
    A = np.asarray(A, dtype=float)
    고윳값, V = np.linalg.eigh(A.T @ A)              # 대칭이므로 eigh
    순서 = np.argsort(고윳값)[::-1]                   # 큰 것부터
    고윳값, V = 고윳값[순서], V[:, 순서]
    특이값 = np.sqrt(np.maximum(고윳값, 0.0))         # 수치오차로 음수가 될 수 있다
    if 눈감아 is None:
        # A^T A 를 거치면 정밀도가 반으로 준다. 그래서 eps 가 아니라 sqrt(eps).
        눈감아 = (np.sqrt(max(A.shape) * np.finfo(float).eps) * 특이값[0]
                if 특이값.size else 0.0)
    r = int(np.sum(특이값 > 눈감아))
    U = np.zeros((A.shape[0], r))
    for i in range(r):
        U[:, i] = A @ V[:, i] / 특이값[i]            # 부호까지 올바르게 정해진다
    return U, 특이값[:r], V[:, :r]

C = np.array([[1.0, 1.0],
              [1.0, 0.0],
              [0.0, 1.0]])
print(show_matrix(C, "앵커 C  (3x2)"))
print(show_matrix(C.T @ C, "C^T C"))
print("C^T C 의 고윳값 :", np.linalg.eigvalsh(C.T @ C), "  (기대 1, 3)")
print()
U, s, V = 손으로SVD(C)
print("직접 구한 특이값 :", np.round(s, 6), "  = sqrt(3), 1 :",
      np.allclose(s, [np.sqrt(3), 1.0]))
print("np.linalg.svd    :", np.round(np.linalg.svd(C, compute_uv=False), 6))
print()
print("U^T U = I 인가 :", np.allclose(U.T @ U, np.eye(len(s))))
print("V^T V = I 인가 :", np.allclose(V.T @ V, np.eye(len(s))))
print("복원 오차       :", np.abs(U @ np.diag(s) @ V.T - C).max())

print(show_matrix(V, "V  (열이 v1, v2)"))
print("v1 = (1,1)/sqrt2 인가 :", np.allclose(np.abs(V[:, 0]), 1/np.sqrt(2)))
print()
print(show_matrix(U, "U  (열이 u1, u2)"))
print("u1 = (2,1,1)/sqrt6 인가 :",
      np.allclose(np.abs(U[:, 0]), np.abs(np.array([2, 1, 1]) / np.sqrt(6))))
print("u2 = (0,1,-1)/sqrt2 인가 :",
      np.allclose(np.abs(U[:, 1]), np.abs(np.array([0, 1, -1]) / np.sqrt(2))))

Q, _ = np.linalg.qr(rng.normal(size=(4, 4)))
W, _ = np.linalg.qr(rng.normal(size=(4, 4)))
print(f"{'참 sigma_min':>14}{'A^T A 를 거쳐서':>18}{'절대오차':>12}{'상대오차':>12}"
      f"{'np.linalg.svd':>18}{'절대오차':>12}")
for 작은값 in (1e-4, 1e-6, 1e-8, 1e-10, 1e-12):
    A0 = Q @ np.diag([1.0, 0.5, 0.1, 작은값]) @ W.T
    거쳐 = np.sqrt(np.maximum(np.linalg.eigvalsh(A0.T @ A0), 0))[0]
    바로 = np.linalg.svd(A0, compute_uv=False)[-1]
    print(f"{작은값:>14.0e}{거쳐:>18.3e}{abs(거쳐-작은값):>12.1e}"
          f"{abs(거쳐-작은값)/작은값:>12.1e}{바로:>18.3e}{abs(바로-작은값):>12.1e}")
print()
print(f"sqrt(eps) = {np.sqrt(np.finfo(float).eps):.3e}   <- 절대오차가 바닥을 치는 자리")
print()
A0 = Q @ np.diag([1.0, 0.5, 0.1, 1e-8]) @ W.T
print(f"cond(A)     = {np.linalg.cond(A0):.4e}")
print(f"cond(A^T A) = {np.linalg.cond(A0.T @ A0):.4e}   <- 제곱되었다")

경우 = (("3x2", C),
        ("2x5 직사각", rng.normal(size=(2, 5))),
        ("랭크 1", np.outer(rng.normal(size=4), rng.normal(size=3))),
        ("영행렬 3x4", np.zeros((3, 4))),
        ("결함 [[3,1],[0,3]]", np.array([[3.0, 1.0], [0.0, 3.0]])),
        ("대칭 [[2,1],[1,2]]", np.array([[2.0, 1.0], [1.0, 2.0]])))
print(f"{'':>22}{'크기':>10}{'랭크':>6}{'특이값':>30}{'복원 오차':>12}")
for 이름, M in 경우:
    Uu, ss, Vv = 손으로SVD(M)
    복원 = Uu @ np.diag(ss) @ Vv.T if len(ss) else np.zeros_like(M)
    보기 = ", ".join(f"{v:.3f}" for v in ss) if len(ss) else "(없음)"
    print(f"{이름:>22}{str(M.shape):>10}{len(ss):>6}{보기:>30}"
          f"{np.abs(복원 - M).max():>12.2e}")

def 따로SVD(A):
    """U 를 A A^T 에서 따로 구한다. 하면 안 되는 방식."""
    값V, V = np.linalg.eigh(A.T @ A)
    값U, U = np.linalg.eigh(A @ A.T)
    V = V[:, np.argsort(값V)[::-1]]
    U = U[:, np.argsort(값U)[::-1]]
    s = np.sqrt(np.maximum(np.sort(값V)[::-1], 0.0))
    return U, s, V

A = np.array([[3.0, 0.0], [4.0, 5.0]])
B = np.array([[1.0, 2.0], [3.0, 4.0]])
for 이름, M in (("A = [[3,0],[4,5]]", A), ("B = [[1,2],[3,4]]", B)):
    Uf, sf, Vf = 따로SVD(M)
    Uo, so, Vo = 손으로SVD(M)
    print(f"{이름}")
    print(f"   따로 구한 것의 복원 오차   : {np.abs(Uf @ np.diag(sf) @ Vf.T - M).max():.3e}")
    print(f"   정의대로 만든 것의 복원 오차 : {np.abs(Uo @ np.diag(so) @ Vo.T - M).max():.3e}")

Uf, sf, Vf = 따로SVD(B)
print("따로 구한 U :\n", np.round(Uf, 6))
print("정의대로 만든 U :\n", np.round(손으로SVD(B)[0], 6))
print()
print("따로 구한 U 도 정규직교이기는 하다 :", np.allclose(Uf.T @ Uf, np.eye(2)))
print("특이값도 정확히 맞다 :", np.allclose(sf, np.linalg.svd(B, compute_uv=False)))
print()
print("그런데 A v_i = sigma_i u_i 라는 끈이 끊어져 있다.")
for i in range(2):
    print(f"   B v{i+1}      = {np.round(B @ Vf[:, i], 6)}")
    print(f"   sigma{i+1} u{i+1} = {np.round(sf[i] * Uf[:, i], 6)}   "
          f"맞는가 : {np.allclose(B @ Vf[:, i], sf[i] * Uf[:, i])}")

def 복원오차(만드는법, X):
    U, s, V = 만드는법(X)
    return float(np.abs(U @ np.diag(s) @ V.T - X).max())

시험 = [rng.normal(size=(4, 4)) for _ in range(500)]
깨짐 = sum(복원오차(따로SVD, X) > 1e-8 for X in 시험)
바름 = sum(복원오차(손으로SVD, X) > 1e-8 for X in 시험)
print(f"무작위 4x4 행렬 500개")
print(f"  따로 구한 U   : 깨진 것 {깨짐}개 ({깨짐/5:.0f}%)")
print(f"  정의대로 만든 U : 깨진 것 {바름}개")

U뭉치, s뭉치, Vt뭉치 = np.linalg.svd(C, full_matrices=True)   # full 이어야 좌영공간까지
r = int(np.sum(s뭉치 > 1e-10))
m, n = C.shape
print(f"C 는 {m}x{n}, 랭크 r = {r}")
print(f"  행공간   기저 : V[:, :{r}]        ({r} 차원, R^{n} 안)")
print(f"  영공간   기저 : V[:, {r}:]        ({n-r} 차원)")
print(f"  열공간   기저 : U[:, :{r}]        ({r} 차원, R^{m} 안)")
print(f"  좌영공간 기저 : U[:, {r}:]        ({m-r} 차원)")
print()
u3 = U뭉치[:, r:]
print("좌영공간 기저 u3 :", np.round(u3.ravel(), 6))
print("  = (-1,1,1)/sqrt3 인가 :",
      np.allclose(np.abs(u3.ravel()), 1/np.sqrt(3)))
print("  C^T u3 =", np.round(C.T @ u3.ravel(), 12), "  <- 좌영공간이 맞다")

# 랭크가 모자란 경우로 네 칸을 다 채워 보자
D = np.array([[1.0, 2.0, 3.0],
              [2.0, 4.0, 6.0],
              [1.0, 1.0, 1.0],
              [0.0, 1.0, 2.0]])                        # 2행 = 1행의 2배
Ud, sd, Vtd = np.linalg.svd(D, full_matrices=True)
rd = int(np.sum(sd > 1e-10))
print(show_matrix(D, "D  (4x3)"))
print("특이값 :", np.round(sd, 6), "  랭크 :", rd)
print()
영공간 = Vtd[rd:].T                                    # 열이 기저
좌영공간 = Ud[:, rd:]
print(f"영공간 : {영공간.shape[1]} 차원 (3 - {rd})")
for i in range(영공간.shape[1]):
    v = 영공간[:, i]
    print(f"  기저 {np.round(v, 6)}  ->  D v = {np.round(D @ v, 12)}")
print(f"좌영공간 : {좌영공간.shape[1]} 차원 (4 - {rd})")
for i in range(좌영공간.shape[1]):
    u = 좌영공간[:, i]
    print(f"  기저 {np.round(u, 6)}  ->  D^T u = {np.round(D.T @ u, 12)}")
print()
print(f"차원 확인 : {rd} + {3-rd} = 3 (열 개수),  {rd} + {4-rd} = 4 (행 개수)")

A = np.array([[3.0, 0.0], [4.0, 5.0]])
Ua, sa, Vta = np.linalg.svd(A)
print("특이값 :", sa, "  = 3sqrt5, sqrt5 :",
      np.allclose(sa, [3*np.sqrt(5), np.sqrt(5)]))
print()
print("sigma1 =", sa[0], "   np.linalg.norm(A, 2) =", np.linalg.norm(A, 2))
print("곱 =", np.prod(sa), "   |det A| =", abs(np.linalg.det(A)))
print()
# 최대 증폭률을 무작위로 확인
최대 = 0.0
for _ in range(20000):
    x = rng.normal(size=2)
    x = x / np.linalg.norm(x)
    최대 = max(최대, float(np.linalg.norm(A @ x)))
print(f"단위벡터 20000개 중 |Ax| 의 최댓값 : {최대:.6f}   sigma1 = {sa[0]:.6f}")
print("v1 을 넣었을 때 :", np.linalg.norm(A @ Vta[0]))

t = np.linspace(0, 2*np.pi, 300)
원 = np.vstack([np.cos(t), np.sin(t)])
단계 = (("① 단위원", np.eye(2)),
        ("② V^T 로 회전", Vta),
        ("③ Sigma 로 늘이기", np.diag(sa) @ Vta),
        ("④ U 로 회전 -> 최종", Ua @ np.diag(sa) @ Vta))
프레임 = []
for 이름, M in 단계:
    곡선 = M @ 원
    프레임.append([
        go.Scatter(x=원[0], y=원[1], mode="lines", name="단위원",
                   line=dict(color="#bbbbbb", width=1.5, dash="dot")),
        go.Scatter(x=곡선[0], y=곡선[1], mode="lines", name="현재",
                   line=dict(color=COLORS["output"], width=4)),
        go.Scatter(x=[0], y=[-7.0], mode="text", text=[이름],
                   textfont=dict(size=17), showlegend=False),
    ])
배치 = layout2d(f"sigma = {np.round(sa,3)},  곱 = {np.prod(sa):.0f} = |det A|",
               extent=7.8)
배치["height"] = 620
slider_figure(프레임, ["1", "2", "3", "4"], 배치, prefix="단계 ")

def 조각들(A):
    """A = sum sigma_i u_i v_i^T 로 쪼갠다."""
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    return [(s[i], np.outer(U[:, i], Vt[i])) for i in range(len(s)) if s[i] > 1e-12]

쌓음 = np.zeros_like(C)
for k, (sig, 조각) in enumerate(조각들(C), start=1):
    쌓음 = 쌓음 + sig * 조각
    print(f"조각 {k}개까지 (sigma = {sig:.4f})  복원 오차 {np.abs(쌓음 - C).max():.3e}")
print()
for k, (sig, 조각) in enumerate(조각들(C), start=1):
    print(show_matrix(sig * 조각, f"sigma{k} u{k} v{k}^T   (랭크 "
                                 f"{np.linalg.matrix_rank(조각)})"))

def 앞k개(A, k):
    """앞 k 개의 조각만 남긴 최적 근사."""
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    return (U[:, :k] * s[:k]) @ Vt[:k]

M = rng.normal(size=(8, 6))
s = np.linalg.svd(M, compute_uv=False)
print(f"{'k':>4}{'|A - A_k|_2':>16}{'sigma_(k+1)':>16}{'무작위 3000개 최선':>22}{'졌는가':>10}")
for k in (1, 2, 3, 4, 5):
    Ak = 앞k개(M, k)
    오차 = np.linalg.norm(M - Ak, 2)
    최선 = np.inf
    for _ in range(3000):
        X = rng.normal(size=(8, k)) @ rng.normal(size=(k, 6))
        최선 = min(최선, float(np.linalg.norm(M - X, 2)))
    print(f"{k:>4}{오차:>16.6f}{s[k]:>16.6f}{최선:>22.6f}"
          f"{str(최선 < 오차 - 1e-9):>10}")

# 더 똑똑한 후보로도 겨뤄 보자 : A 의 앞 k 열 / 무작위 k 열
k = 3
Ak = 앞k개(M, k)
print(f"랭크 {k} 후보들의 |A - X|_2")
print(f"  A_k (SVD 앞 {k}개)          : {np.linalg.norm(M - Ak, 2):.6f}   <- 최선")
앞열 = M.copy(); 앞열[:, k:] = 0
print(f"  A 의 앞 {k}열만 남기기       : {np.linalg.norm(M - 앞열, 2):.6f}")
Q, _ = np.linalg.qr(M[:, :k])
print(f"  앞 {k}열의 공간으로 투영     : "
      f"{np.linalg.norm(M - Q @ (Q.T @ M), 2):.6f}")
print()
print("'앞 k 열을 남기는 것' 은 저계수 근사가 아니다.")

def 합성사진(n=192, seed=29):
    """자연스러운 구조를 가진 시험용 이미지."""
    rng2 = np.random.default_rng(seed)
    y, x = np.mgrid[0:n, 0:n] / (n - 1)
    그림 = 0.42 + 0.16 * np.cos(2.0*x + 0.5) * np.cos(1.5*y - 0.2)
    그림 += 0.26 * np.exp(-((x-0.28)**2 + (y-0.30)**2) / 0.014)
    그림 -= 0.22 * np.exp(-((x-0.74)**2 + (y-0.70)**2) / 0.010)
    그림 += 0.22 * ((x > 0.58) & (x < 0.92) & (y > 0.08) & (y < 0.34))
    그림 -= 0.28 * (np.abs(y - x - 0.18) < 0.035)
    그림 += 0.16 * (np.abs(y + x - 1.55) < 0.028)
    그림 += 0.09 * np.sin(24*x) * ((y > 0.78) & (y < 0.96))
    그림 += 0.012 * rng2.standard_normal((n, n))
    return np.clip(그림, 0, 1)

사진 = 합성사진()
n = 사진.shape[0]
Ui, Si, Vti = np.linalg.svd(사진, full_matrices=False)
print(f"{n} x {n} = {n*n:,} 개의 수")
print("특이값 상위 8 :", np.round(Si[:8], 3))
print()
print(f"{'k':>5}{'저장할 수':>12}{'원본 대비':>10}{'상대 오차':>12}{'담은 에너지':>12}")
전체 = (Si**2).sum()
for k in (1, 3, 10, 25, 60):
    Ak = 앞k개(사진, k)
    저장 = k * (2*n + 1)
    print(f"{k:>5}{저장:>12,}{저장/(n*n):>10.1%}"
          f"{np.linalg.norm(사진-Ak)/np.linalg.norm(사진):>12.4f}"
          f"{(Si[:k]**2).sum()/전체:>12.4f}")

보기 = [1, 2, 3, 5, 8, 12, 20, 30, 45, 60, 90, 130]
프레임 = []
for k in 보기:
    Ak = np.clip(앞k개(사진, k), 0, 1)
    프레임.append([go.Heatmap(z=Ak[::-1], colorscale="gray", zmin=0, zmax=1,
                             showscale=False)])
배치 = dict(title=dict(text="랭크를 올려 가며"),
           xaxis=dict(visible=False, scaleanchor="y"),
           yaxis=dict(visible=False),
           height=560, margin=dict(l=40, r=40, t=60, b=40))
slider_figure(프레임, 보기, 배치, prefix="랭크 k = ", initial=4)

무작위 = rng.random((n, n))
S무 = np.linalg.svd(무작위, compute_uv=False)
전체무 = (S무**2).sum()

print(f"{'':>10}{'99% 에너지에 필요한 k':>24}{'그때 저장량':>14}")
for 이름, sv, 합 in (("사진", Si, 전체), ("무작위", S무, 전체무)):
    k = int(np.searchsorted(np.cumsum(sv**2)/합, 0.99) + 1)
    print(f"{이름:>10}{k:>24}{k*(2*n+1)/(n*n):>14.1%}")
print()
print("무작위 잡음은 SVD 로 압축하면 원본보다 커진다.")
print("자연의 데이터가 압축되는 것은 그것이 무작위가 아니기 때문이다.")

go.Figure(
    data=[go.Scatter(x=np.arange(1, n+1), y=Si, mode="lines", name="사진",
                     line=dict(color=COLORS["output"], width=3)),
          go.Scatter(x=np.arange(1, n+1), y=S무, mode="lines", name="무작위 잡음",
                     line=dict(color="#999999", width=2, dash="dash"))],
    layout=go.Layout(title=dict(text="특이값 스펙트럼"),
                     xaxis=dict(title=dict(text="i")),
                     yaxis=dict(title=dict(text="sigma_i"), type="log"),
                     height=420, margin=dict(l=70, r=20, t=60, b=50)))

U3, s3, Vt3 = np.linalg.svd(D, full_matrices=False)   # D 는 4x3
print("np.linalg.svd 가 돌려주는 세 번째 것의 크기 :", Vt3.shape)
print("  이것은 V 가 아니라 V^T 이다.")
print()
print("올바른 복원 : U @ diag(s) @ Vt        오차",
      np.abs(U3 @ np.diag(s3) @ Vt3 - D).max())
print("틀린 복원   : U @ diag(s) @ Vt.T      오차",
      np.abs(U3 @ np.diag(s3) @ Vt3.T - D).max())
print()
print("정사각이 아니면 아예 곱해지지도 않아 금방 들킨다.")
E = rng.normal(size=(2, 5))
Ue, se, Vte = np.linalg.svd(E, full_matrices=False)
try:
    Ue @ np.diag(se) @ Vte.T
except ValueError as 오류:
    print("  2x5 에서는 :", str(오류).split(",")[0])
print("정사각일 때가 오히려 위험하다. 조용히 틀린 답이 나온다.")
print()
print("full_matrices 의 차이")
for 값 in (True, False):
    Uf, sf, Vtf = np.linalg.svd(C, full_matrices=값)
    print(f"  full_matrices={값!s:>5} : U {Uf.shape}, Vt {Vtf.shape}"
          f"   좌영공간까지 주는가 : {Uf.shape[1] == C.shape[0]}")

print("조건수 = sigma1 / sigma_n")
for 이름, M in (("[[2,1],[1,2]]", np.array([[2.0, 1.0], [1.0, 2.0]])),
                ("[[3,0],[4,5]]", A),
                ("힐베르트 6x6", np.array([[1.0/(i+j+1) for j in range(6)]
                                          for i in range(6)]))):
    sv = np.linalg.svd(M, compute_uv=False)
    print(f"  {이름:>16} : sigma1 {sv[0]:>10.4g}, sigma_n {sv[-1]:>10.4g}, "
          f"조건수 {sv[0]/sv[-1]:>12.4g}")
print()
print("np.linalg.cond 와 같은가 :",
      np.isclose(np.linalg.cond(A), np.linalg.svd(A, compute_uv=False)[0]
                 / np.linalg.svd(A, compute_uv=False)[-1]))

print(f"{'':>22}{'|고윳값|':>20}{'특이값':>20}{'같은가':>8}{'A A^T = A^T A':>16}")
for 이름, M in (("대칭 [[2,1],[1,2]]", np.array([[2.0, 1.0], [1.0, 2.0]])),
                ("회전 90도", np.array([[0.0, -1.0], [1.0, 0.0]])),
                ("결함 [[3,1],[0,3]]", np.array([[3.0, 1.0], [0.0, 3.0]])),
                ("[[3,0],[4,5]]", A)):
    고 = np.sort(np.abs(np.linalg.eigvals(M)))[::-1]
    특 = np.linalg.svd(M, compute_uv=False)
    print(f"{이름:>22}{str(np.round(고,4)):>20}{str(np.round(특,4)):>20}"
          f"{str(np.allclose(고, 특)):>8}{str(np.allclose(M @ M.T, M.T @ M)):>16}")