보강 1. 주성분분석과 SVD — 파이썬 실습

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, show_matrix, slider_figure

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(35)
print("numpy", np.__version__)
print("sklearn 은 이 환경에 없다. 전부 numpy 로 직접 만든다.")

X = np.array([[1.0, 3.0],
              [3.0, 1.0],
              [4.0, 3.0],
              [5.0, 3.0],
              [7.0, 5.0]])
m, n = X.shape
mu = X.mean(axis=0)
Xc = X - mu
print(show_matrix(X, "X  (표본 5개, 변수 2개)"))
print("평균 mu =", mu, "   (4, 3) 인가 :", np.allclose(mu, [4, 3]))
print(show_matrix(Xc, "Xc = X - mu"))
print("열의 합이 0 인가 :", np.allclose(Xc.sum(axis=0), 0))

C = Xc.T @ Xc / (m - 1)
print(show_matrix(Xc.T @ Xc, "Xc^T Xc"))
print(show_matrix(C, "C = Xc^T Xc / (m-1)"))
print("서술의 [[5,2],[2,2]] 인가 :", np.allclose(C, [[5, 2], [2, 2]]))
print("np.cov 와 같은가 :", np.allclose(C, np.cov(X, rowvar=False)))
print()
값, Q = np.linalg.eigh(C)
순 = np.argsort(값)[::-1]
값, Q = 값[순], Q[:, 순]
print("고윳값 :", 값, "   (6, 1) 인가 :", np.allclose(값, [6, 1]))
q1 = np.array([2.0, 1.0]) / np.sqrt(5)
q2 = np.array([-1.0, 2.0]) / np.sqrt(5)
print("q1 = (2,1)/sqrt5 인가 :", np.allclose(np.abs(Q[:, 0]), np.abs(q1)))
print("q2 = (-1,2)/sqrt5 인가 :", np.allclose(np.abs(Q[:, 1]), np.abs(q2)))
print("q1 . q2 =", float(q1 @ q2), "  <- 직교는 스펙트럼 정리가 준다")
print()
print("제1주성분의 각도 :", f"{np.degrees(np.arctan2(1, 2)):.3f} 도")
print("대각합 =", np.trace(C), " = 5+2 = 6+1 :", np.isclose(np.trace(C), 7))

각도 = np.linspace(0, np.pi, 721)
분산 = np.array([float(np.array([np.cos(t), np.sin(t)]) @ C
                      @ np.array([np.cos(t), np.sin(t)])) for t in 각도])
최대, 최소 = int(np.argmax(분산)), int(np.argmin(분산))
print(f"가장 큰 분산 : {분산[최대]:.6f} at {np.degrees(각도[최대]):.3f} 도")
print(f"가장 작은 분산 : {분산[최소]:.6f} at {np.degrees(각도[최소]):.3f} 도")
print(f"두 각도의 차이 : {np.degrees(각도[최소] - 각도[최대]):.3f} 도  <- 90 도")
print(f"두 값의 합 : {분산[최대] + 분산[최소]:.6f}  = trace(C) = {np.trace(C)}")
print()
print("고윳값과 같은가 :", np.isclose(분산[최대], 6, atol=1e-4),
      np.isclose(분산[최소], 1, atol=1e-4))

프레임, 이름표 = [], []
반지름 = np.sqrt(분산)
for i in range(0, 721, 40):
    t = 각도[i]
    v = np.array([np.cos(t), np.sin(t)])
    투영 = (Xc @ v)[:, None] * v[None, :]
    프레임.append([
        go.Scatter(x=반지름*np.cos(각도), y=반지름*np.sin(각도), mode="lines",
                   line=dict(color="#cccccc", width=1.5), name="모든 방향"),
        go.Scatter(x=Xc[:, 0], y=Xc[:, 1], mode="markers",
                   marker=dict(size=11, color=COLORS["input"]), name="Xc"),
        go.Scatter(x=[-4*v[0], 4*v[0]], y=[-4*v[1], 4*v[1]], mode="lines",
                   line=dict(color=COLORS["output"], width=3), name="방향 v"),
        go.Scatter(x=투영[:, 0], y=투영[:, 1], mode="markers",
                   marker=dict(size=9, color="#2ca02c", symbol="square"),
                   name=f"분산 {분산[i]:.3f}"),
    ])
    이름표.append(f"{np.degrees(t):.0f}")
배치 = dict(title=dict(text="자를 돌려 가며 그림자의 분산을 잰다"),
           xaxis=dict(range=[-4.5, 4.5], title=dict(text="x")),
           yaxis=dict(range=[-4.5, 4.5], scaleanchor="x", title=dict(text="y")),
           height=560, margin=dict(l=60, r=30, t=60, b=50))
slider_figure(프레임, 이름표, 배치, prefix="각도 ", initial=0)

U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
print("Xc 의 특이값 :", s)
print("sqrt(24), 2 인가 :", np.allclose(s, [np.sqrt(24), 2.0]))
print()
print("sigma^2 / (m-1) =", s**2 / (m - 1), "  고윳값 :", 값)
print("같은가 :", np.allclose(s**2 / (m - 1), 값))
print()
print("V 의 행이 고유벡터인가 :",
      np.allclose(np.abs(Vt[0]), np.abs(q1)), np.allclose(np.abs(Vt[1]), np.abs(q2)))
print()
점수_eigh = Xc @ Q
점수_svd = U * s
print(show_matrix(점수_eigh, "점수 (eigh 경유)  Xc Q"))
print(show_matrix(점수_svd, "점수 (svd 경유)  U Sigma"))
print("부호까지 같은가 :", np.allclose(점수_eigh, 점수_svd))
print("부호를 빼면 같은가 :", np.allclose(np.abs(점수_eigh), np.abs(점수_svd)))
print()
print("-> 부호는 규약이다. 고유벡터의 방향은 어차피 정해지지 않는다(L21).")

비율 = s**2 / np.sum(s**2)
print("설명분산 비율 :", 비율, "  = 6/7, 1/7 :",
      np.allclose(비율, [6/7, 1/7]))
print("누적 :", np.cumsum(비율))
print()
print("고윳값으로 계산해도 같은가 :", np.allclose(값 / np.sum(값), 비율))
print()
# 에크하르트-영 : 랭크 1 근사가 최선인가
A1 = np.outer(U[:, 0], Vt[0]) * s[0]
print("랭크 1 근사의 오차 :", f"{np.linalg.norm(Xc - A1):.6f}",
      "  sigma_2 =", f"{s[1]:.6f}")
print("같은가 :", np.isclose(np.linalg.norm(Xc - A1, 2), s[1]))
최선 = np.inf
for _ in range(3000):
    Z = np.outer(rng.normal(size=m), rng.normal(size=n))
    최선 = min(최선, float(np.linalg.norm(Xc - Z, 2)))
print(f"무작위 랭크 1 행렬 3000개 중 최선 : {최선:.6f}  (진 것인가 : {최선 < s[1]})")

def 자료(kappa, mm=400, nn=8, seed=0):
    """특이값을 정확히 아는 중심화된 자료를 만든다.

    먼저 중심화한 뒤 QR 을 하면 Q 의 열도 이미 중심화되어 있다.
    """
    r = np.random.default_rng(seed)
    B = r.normal(size=(mm, nn))
    B = B - B.mean(axis=0)
    Q1, _ = np.linalg.qr(B)
    Vv, _ = np.linalg.qr(r.normal(size=(nn, nn)))
    시그마 = np.logspace(0, -np.log10(kappa), nn)
    return Q1 @ np.diag(시그마) @ Vv.T, 시그마

print("조건수가 정확히 제곱되는가")
print(f"{'k(Xc)':>9}{'k(C) 실측':>14}{'k(Xc)^2':>14}{'비':>8}   비고")
for k in (1e2, 1e4, 1e6, 1e7, 1e8):
    Xk, _ = 자료(k)
    kX, kC = np.linalg.cond(Xk), np.linalg.cond(Xk.T @ Xk)
    비고 = "1/eps 를 넘었다" if kX**2 > 1/np.finfo(float).eps else ""
    print(f"{k:>9.0e}{kC:>14.3e}{kX**2:>14.3e}{kC/kX**2:>8.4f}   {비고}")
print()
print("-> 1e7 까지는 비가 1.0000 이다. 정확히 제곱된다.")
print("   마지막 줄만 0.76 으로 어긋나는데, 이것이 오차가 아니라 이 절의 결론이다.")
print("   kappa(Xc)^2 이 1/eps 를 넘어 버려 kappa(C) 를 '재는 일' 자체가 이미")
print("   부정확해졌다. 제곱이 정확히 되지 않는 것이 아니라, 제곱된 값을 컴퓨터가")
print("   더 이상 붙들지 못하는 것이다.")

print("가장 작은 특이값을 얼마나 정확히 되찾는가  (400 x 8)")
print(f"{'k(Xc)':>9}{'k(C)':>12}{'공분산 경유':>14}{'SVD 직접':>13}   비고")
for k in (1e2, 1e4, 1e6, 1e8, 1e10):
    Xk, 참 = 자료(k)
    Ck = Xk.T @ Xk
    람 = np.linalg.eigvalsh(Ck)
    경유 = np.sqrt(max(람[0], 0.0))
    직접 = np.linalg.svd(Xk, compute_uv=False)[-1]
    비고 = "고윳값이 음수" if 람[0] < 0 else ""
    print(f"{k:>9.0e}{np.linalg.cond(Ck):>12.1e}"
          f"{abs(경유-참[-1])/참[-1]:>14.1e}{abs(직접-참[-1])/참[-1]:>13.1e}   {비고}")
print()
print(f"1/eps = {1/np.finfo(float).eps:.2e}")
print("-> k(Xc) 가 1e8 을 넘으면 k(C) 가 1/eps 를 넘어 공분산이 무너진다.")
print("   양의 준정부호인데 고윳값이 음수로 나온다. 수학이 아니라 반올림의 일이다.")

print(f"eps(float64) = {np.finfo(np.float64).eps:.2e}  ->  1/sqrt(eps) = "
      f"{1/np.sqrt(np.finfo(np.float64).eps):.1e}")
print(f"eps(float32) = {np.finfo(np.float32).eps:.2e}  ->  1/sqrt(eps) = "
      f"{1/np.sqrt(np.finfo(np.float32).eps):.0f}   (약 10^3.5)")
print()
r2 = np.random.default_rng(1)
U2, _ = np.linalg.qr(r2.normal(size=(200, 200)))
V2, _ = np.linalg.qr(r2.normal(size=(2, 2)))
print(f"{'sigma_min':>11}{'kappa':>9}{'공분산 경유':>16}{'상대오차':>11}{'SVD 직접':>16}")
for sm in (1e-3, 3e-4, 1e-4):
    Xf = (U2[:, :2] @ np.diag([1.0, sm]) @ V2.T).astype(np.float32)
    람 = np.linalg.eigvalsh((Xf.T @ Xf).astype(np.float64))
    경유 = float(np.sqrt(max(람[0], 0.0)))
    직접 = float(np.linalg.svd(Xf, compute_uv=False)[-1])
    print(f"{sm:>11.0e}{1/sm:>9.0e}{경유:>16.4e}{abs(경유-sm)/sm:>11.1%}{직접:>16.4e}")
print()
print("-> kappa = 1e4 에서 공분산 경유가 0 을 준다. 제2주성분의 크기가 0 으로 보고된다.")
print("   기계학습에서 float32 는 흔하고, 조건수 1e4 도 드물지 않다.")

def 제1주성분(D, 중심화=True):
    E = D - D.mean(axis=0) if 중심화 else D
    _, _, Vt2 = np.linalg.svd(E, full_matrices=False)
    return Vt2[0]

바르게 = 제1주성분(X, True)
빼먹고 = 제1주성분(X, False)
평균방향 = mu / np.linalg.norm(mu)
사이 = lambda a, b: float(np.degrees(np.arccos(np.clip(abs(a @ b), 0, 1))))
print("제대로 한 제1주성분 :", 바르게, f"  (각도 {np.degrees(np.arctan2(*바르게[::-1])):.2f} 도)")
print("중심화를 뺀 제1주성분 :", 빼먹고)
print("평균 방향           :", 평균방향)
print()
print(f"빼먹은 것과 평균 방향의 사잇각   : {사이(빼먹고, 평균방향):.3f} 도  <- 거의 붙었다")
print(f"빼먹은 것과 제대로 한 것의 사잇각 : {사이(빼먹고, 바르게):.3f} 도")
print()
print("-> 중심화를 빼면 X^T X 에 m * mu mu^T 가 얹힌다.")
print("   그 항은 랭크 1 이고 방향이 mu 이므로, 평균이 클수록 답을 그쪽으로 끌어당긴다.")

print("평균을 멀리 옮길수록 얼마나 끌려가는가")
print(f"{'평균 이동':>10}{'중심화한 쪽의 각도':>20}{'뺀 쪽의 각도':>16}"
      f"{'평균 방향과의 각':>18}")
for 배 in (0.0, 1.0, 3.0, 10.0, 50.0):
    이동 = X + 배 * np.array([4.0, 3.0])
    빼 = 제1주성분(이동, False)
    참 = 제1주성분(이동, True)
    if 참[0] < 0: 참 = -참
    if 빼[0] < 0: 빼 = -빼
    mu2 = 이동.mean(axis=0); mu2 = mu2 / np.linalg.norm(mu2)
    print(f"{배:>10.0f}{np.degrees(np.arctan2(참[1], 참[0])):>20.4f}"
          f"{np.degrees(np.arctan2(빼[1], 빼[0])):>16.4f}"
          f"{사이(빼, mu2):>18.3f}")
print()
print("-> 가운데 칸이 26.5651 로 꿈쩍도 하지 않는다. 중심화한 PCA 는 옮김에 불변이다.")
print("   오른쪽 두 칸은 평균을 옮길수록 평균 방향에 달라붙는다.")
print("   평균을 50배 옮기면 사잇각이 0.001 도다. 사실상 평균 방향 그 자체다.")

# y 를 x 로 회귀 / x 를 y 로 회귀 / PCA
A회귀 = np.column_stack([np.ones(m), Xc[:, 0]])
기울기_yx = np.linalg.lstsq(A회귀, Xc[:, 1], rcond=None)[0][1]
B회귀 = np.column_stack([np.ones(m), Xc[:, 1]])
기울기_xy = np.linalg.lstsq(B회귀, Xc[:, 0], rcond=None)[0][1]
기울기_pca = 바르게[1] / 바르게[0]
print(f"y 를 x 로 회귀   : 기울기 {기울기_yx:.6f}   (세로 거리의 제곱합을 최소로)")
print(f"x 를 y 로 회귀   : 기울기 {1/기울기_xy:.6f}   (가로 거리의 제곱합을 최소로)")
print(f"PCA 제1주성분    : 기울기 {기울기_pca:.6f}   (수직 거리의 제곱합을 최소로)")
print()
def 거리합(기울기, 방식):
    v = np.array([1.0, 기울기]); v = v / np.linalg.norm(v)
    if 방식 == "수직":
        return float(np.sum((Xc - (Xc @ v)[:, None] * v)**2))
    if 방식 == "세로":
        return float(np.sum((Xc[:, 1] - 기울기 * Xc[:, 0])**2))
    return float(np.sum((Xc[:, 0] - Xc[:, 1] / 기울기)**2))
print(f"{'':>16}{'수직 거리^2':>14}{'세로 거리^2':>14}{'가로 거리^2':>14}")
for 이름, g in (("y~x 회귀", 기울기_yx), ("x~y 회귀", 1/기울기_xy),
                ("PCA", 기울기_pca)):
    print(f"{이름:>16}{거리합(g,'수직'):>14.6f}{거리합(g,'세로'):>14.6f}"
          f"{거리합(g,'가로'):>14.6f}")
print()
print("-> 각 방법이 자기 칸에서만 최소다. 셋은 다른 문제를 푸는 것이다.")
print("   PCA 의 수직 거리 제곱합은 lambda_2 * (m-1) =",
      f"{값[1]*(m-1):.6f} 와 같다 :", np.isclose(거리합(기울기_pca,'수직'), 값[1]*(m-1)))

print("둘째 변수의 단위만 바꿔 본다 (예: cm -> mm)")
print(f"{'배율':>8}{'제1주성분':>26}{'각도':>10}{'설명분산 1위':>14}")
for 배 in (1.0, 2.0, 10.0, 100.0):
    Y = X * np.array([1.0, 배])
    v = 제1주성분(Y, True)
    Yc = Y - Y.mean(axis=0)
    ss = np.linalg.svd(Yc, compute_uv=False)
    print(f"{배:>8.0f}{str(np.round(v, 4)):>26}"
          f"{np.degrees(np.arctan2(v[1], v[0])):>10.2f}"
          f"{ss[0]**2/np.sum(ss**2):>14.4f}")
print()
print("-> 축 하나의 눈금만 바꿨는데 제1주성분이 그쪽으로 끌려간다.")
print("   PCA 는 분산을 보는데, 분산은 단위에 딸린 양이기 때문이다.")
print("   그래서 단위가 다른 변수를 섞을 때는 표준화(상관행렬 PCA)를 한다.")

print("표준화하면 눈금에 흔들리지 않는가")
def 표준화PCA(D):
    E = D - D.mean(axis=0)
    E = E / E.std(axis=0, ddof=1)
    return np.linalg.svd(E, full_matrices=False)[2][0]
print(f"{'배율':>8}{'표준화 뒤 제1주성분':>28}{'각도':>10}")
for 배 in (1.0, 2.0, 10.0, 100.0):
    v = 표준화PCA(X * np.array([1.0, 배]))
    if v[0] < 0: v = -v
    print(f"{배:>8.0f}{str(np.round(v, 6)):>28}"
          f"{np.degrees(np.arctan2(v[1], v[0])):>10.3f}")
print()
print("-> 배율을 100배 해도 답이 그대로다. 표준화가 눈금을 지워 버리기 때문이다.")
print("   대신 '어느 변수가 원래 더 크게 흔들렸는가' 라는 정보도 함께 지워진다.")