SVD·PCA·저랭크 근사

SVD

정의. 모든 $m \times n$ 실행렬 $\mathbf{A}$는

로 분해된다. $\mathbf{U}$($m \times m$)와 $\mathbf{V}$($n \times n$)는 직교행렬이고, $\boldsymbol{\Sigma}$($m \times n$)는 대각에 $\sigma_1 \ge \sigma_2 \ge \cdots \ge 0$을 놓고 나머지가 0인 행렬이다. $\sigma_i$를 특이값(singular value), $\mathbf{U}$의 열 $\mathbf{u}_i$를 왼쪽 특이벡터, $\mathbf{V}$의 열 $\mathbf{v}_i$를 오른쪽 특이벡터라 한다. 0이 아닌 특이값이 $r$개이면 처음 $r$개 항만 남긴 thin SVD $\mathbf{A} = \mathbf{U}_r\boldsymbol{\Sigma}_r\mathbf{V}_r^\top$로 충분하다. 열마다 쓰면

이다. 입력 공간의 직교하는 축 $\mathbf{v}_i$가 출력 공간의 직교하는 축 $\mathbf{u}_i$로 가면서 $\sigma_i$배로 늘어난다.

기하: 회전–스케일–회전. $\mathbf{A}\mathbf{x} = \mathbf{U}(\boldsymbol{\Sigma}(\mathbf{V}^\top\mathbf{x}))$를 오른쪽부터 읽는다.

1. $\mathbf{V}^\top$: 입력을 회전(또는 반사)해 축 $\mathbf{v}_i$를 좌표축에 맞춘다.
2. $\boldsymbol{\Sigma}$: 좌표축마다 $\sigma_i$배로 늘이거나 줄인다. 차원이 바뀌면 남는 축을 버리거나 0을 채운다.
3. $\mathbf{U}$: 결과를 출력 공간에서 회전해 축을 $\mathbf{u}_i$에 맞춘다.

그래서 단위원(단위구)은 반축이 $\sigma_i\mathbf{u}_i$인 타원(타원체)이 된다. 7.4장의 spectral 노름 $\lVert\mathbf{A}\rVert_2 = \sigma_1$은 "가장 많이 늘어나는 방향의 배율"이다. 아래 위젯에서 세 단계를 하나씩 재생해 본다.

고유분해와의 관계. $\mathbf{A}^\top\mathbf{A} = \mathbf{V}\boldsymbol{\Sigma}^\top\boldsymbol{\Sigma}\mathbf{V}^\top$이고 $\mathbf{A}\mathbf{A}^\top = \mathbf{U}\boldsymbol{\Sigma}\boldsymbol{\Sigma}^\top\mathbf{U}^\top$이다. 두 행렬은 대칭 PSD이므로 스펙트럴 정리가 적용되고, 위 식이 바로 그 고유분해다. 즉 오른쪽 특이벡터는 $\mathbf{A}^\top\mathbf{A}$의 고유벡터, 특이값은 그 고윳값의 제곱근이다. SVD의 존재는 이 사실에서 나온다. $\mathbf{A}^\top\mathbf{A}$를 직교 대각화해 $\mathbf{v}_i$와 $\sigma_i$를 얻고, $\sigma_i > 0$인 $i$에 대해 $\mathbf{u}_i = \mathbf{A}\mathbf{v}_i/\sigma_i$로 두면 이들이 정규직교임을 확인할 수 있다. 실제 계산에서는 $\mathbf{A}^\top\mathbf{A}$를 만들지 않는다. 7.4장에서 본 것처럼 조건수가 제곱되기 때문이다.

네 부분공간과 pseudoinverse. 특이벡터는 7.3장의 네 부분공간에 정규직교 basis를 한꺼번에 준다.

rank는 0이 아닌 특이값의 개수이고, 수치적 rank는 허용오차보다 큰 특이값의 개수다(7.3장). 7.4장의 pseudoinverse도 SVD로 한 줄에 쓴다.

$\boldsymbol{\Sigma}^+$는 $\boldsymbol{\Sigma}$를 전치하고 0이 아닌 $\sigma_i$를 $1/\sigma_i$로 바꾼 $n \times m$ 행렬이다. column space 성분은 $\sigma_i$로 나눠 되돌리고, 되돌릴 수 없는 방향(left null space)은 버린다. 그 결과는 row space 안에 놓이므로 최소 노름 해가 된다.

import numpy as np

A = np.array([[3.0, 0.0],
              [4.0, 5.0]])                               # (2, 2)
U, s, Vt = np.linalg.svd(A)                              # (2, 2), (2,), (2, 2)
print("특이값          :", np.round(s, 4))
print("U S V^T = A ?   :", np.allclose(U @ np.diag(s) @ Vt, A))
print("A v1 = s1 u1 ?  :", np.allclose(A @ Vt[0], s[0] * U[:, 0]))
print("sqrt(eig(A^T A)):", np.round(np.sqrt(np.linalg.eigvalsh(A.T @ A))[::-1], 4))

theta = np.linspace(0, 2 * np.pi, 3601)
circle = np.stack([np.cos(theta), np.sin(theta)])        # (2, 3601) 단위원
r = np.linalg.norm(A @ circle, axis=0)                   # 상의 원점 거리
print(f"단위원의 상: 최대 반지름 {r.max():.4f}, 최소 반지름 {r.min():.4f}")
A_pinv = Vt.T @ np.diag(1 / s) @ U.T                     # V S^+ U^T
print("V S^+ U^T = pinv(A) ?", np.allclose(A_pinv, np.linalg.pinv(A)))

특이값          : [6.7082 2.2361]
U S V^T = A ?   : True
A v1 = s1 u1 ?  : True
sqrt(eig(A^T A)): [6.7082 2.2361]
단위원의 상: 최대 반지름 6.7082, 최소 반지름 2.2361
V S^+ U^T = pinv(A) ? True

NumPy의 svd는 $\mathbf{V}$가 아니라 $\mathbf{V}^\top$(Vt)를 돌려준다. 그래서 오른쪽 특이벡터 $\mathbf{v}_1$은 Vt[0], 즉 행이다. 자주 틀리는 부분이다. 단위원 위 3601개 점의 상에서 원점까지의 최대·최소 거리가 두 특이값 $3\sqrt{5}$와 $\sqrt{5}$와 일치한다.

저랭크 근사

왜 필요한가. 7.3장에서 데이터와 가중치가 근사적으로 저랭크인 경우가 많다고 했다. 그렇다면 주어진 rank $k$에서 원래 행렬에 가장 가까운 행렬은 무엇이고, 그 근사의 오차는 얼마일까. 이 질문에 정확한 답이 있다는 것이 SVD를 실용적으로 만든다.

rank-1 행렬의 합. SVD를 열과 행의 곱으로 펼치면

이다. 각 항 $\mathbf{u}_i\mathbf{v}_i^\top$는 rank 1인 행렬이고, 특이값이 그 항의 크기다. 앞의 $k$개만 남긴 것을 $\mathbf{A}_k = \sum_{i=1}^k \sigma_i\mathbf{u}_i\mathbf{v}_i^\top$(truncated SVD)라 한다.

정리(Eckart–Young). spectral 노름에 대한 부분까지 포함해 Eckart–Young–Mirsky 정리라고도 한다. rank가 $k$ 이하인 모든 행렬 $\mathbf{B}$에 대해

즉 truncated SVD는 Frobenius 노름과 spectral 노름 모두에서 최적의 rank-$k$ 근사이고, 그 오차는 버린 특이값으로 정확히 계산된다. 탐욕적으로 큰 항부터 고르는 것이 전역 최적이라는 점이 핵심이다.

에너지 비율. $\lVert\mathbf{A}\rVert_F^2 = \sum_i \sigma_i^2$이므로, rank-$k$ 근사가 담는 비율

를 흔히 "에너지" 비율이라 부른다. 데이터를 중심화한 경우에는 이 비율이 PCA의 설명 분산 비율과 같다(다음 절).

공간 60점 × 시각 40개의 장(field)을 만들어 확인한다. 세 개의 공간 패턴이 각자 다른 시간 변화를 하고, 작은 잡음이 섞여 있다.

rng = np.random.default_rng(0)
x, t = np.linspace(0, 1, 60), np.linspace(0, 1, 40)
X = (np.outer(np.sin(np.pi * x), np.exp(-t))
     + 0.5 * np.outer(np.sin(2 * np.pi * x), np.cos(3 * t))
     + 0.2 * np.outer(np.sin(3 * np.pi * x), t)
     + 0.01 * rng.normal(size=(60, 40)))                 # (60, 40) 세 패턴 + 잡음
U, s, Vt = np.linalg.svd(X, full_matrices=False)
print("특이값 앞 6개:", np.round(s[:6], 3))
for k in [1, 2, 3, 4]:
    Xk = U[:, :k] @ np.diag(s[:k]) @ Vt[:k]              # (60, 40) rank-k 근사
    err = np.linalg.norm(X - Xk)
    G = np.linalg.qr(rng.normal(size=(60, k)))[0]        # 무작위 k차원 부분공간
    err_rand = np.linalg.norm(X - G @ (G.T @ X))
    energy = (s[:k] ** 2).sum() / (s ** 2).sum()
    print(f"k = {k}  오차 {err:.4f}  sqrt(꼬리 s^2 합) {np.sqrt((s[k:] ** 2).sum()):.4f}  "
          f"무작위 부분공간 {err_rand:.4f}  에너지 {energy:.4%}")

특이값 앞 6개: [23.271 11.518  0.42   0.137  0.124  0.123]
k = 1  오차 11.5352  sqrt(꼬리 s^2 합) 11.5352  무작위 부분공간 25.8801  에너지 80.2759%
k = 2  오차 0.6251  sqrt(꼬리 s^2 합) 0.6251  무작위 부분공간 24.3924  에너지 99.9421%
k = 3  오차 0.4626  sqrt(꼬리 s^2 합) 0.4626  무작위 부분공간 25.2986  에너지 99.9683%
k = 4  오차 0.4418  sqrt(꼬리 s^2 합) 0.4418  무작위 부분공간 25.1137  에너지 99.9711%

근사 오차가 버린 특이값으로 계산한 값과 넷째 자리까지 같다. 같은 rank의 무작위 부분공간에 사영하면 오차가 훨씬 크다. 특이값은 셋째(0.42)까지 뚜렷하고 넷째부터는 0.12–0.14 근처로 평평하다. 이 평평한 바닥이 잡음이다. 셋째 패턴은 에너지로는 0.03%에 불과해 "에너지 99.9%" 같은 기준으로는 버려진다. 에너지 비율만으로 $k$를 고르면 작지만 실제로 있는 구조를 놓칠 수 있으므로, 특이값의 감소 양상을 함께 본다.

PCA

왜 필요한가. 특징이 수십·수백 개인 데이터에서 "변동의 대부분이 몇 개의 방향에 몰려 있는가", "그 방향은 무엇인가"를 묻는 것은 탐색적 분석의 출발점이다. 주성분 분석(PCA, principal component analysis)은 이 질문에 답한다. 7.5장에서 공분산 행렬의 고유벡터가 주축이라는 것을 보았다. 이 절은 같은 답을 데이터 행렬의 SVD로 구한다.

정의. 데이터 행렬 $\mathbf{X} \in \mathbb{R}^{N \times D}$(행 = 샘플)에서 열 평균을 뺀 중심화 행렬 $\mathbf{X}_c$를 만든다. 그 thin SVD가 $\mathbf{X}_c = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top$이면, 표본 공분산은

이다. 이것이 $\mathbf{C}$의 고유분해이므로 다음이 성립한다.

• 주성분 방향(loading)은 오른쪽 특이벡터 $\mathbf{v}_i$이다.
• $i$번째 방향의 분산은 $\sigma_i^2/(N - 1)$이다.
• 각 샘플의 주성분 점수(score)는 $\mathbf{X}_c\mathbf{V} = \mathbf{U}\boldsymbol{\Sigma}$이다.

두 관점이 같은 답을 준다. PCA는 두 가지로 정의할 수 있다.

• 분산 최대화. 사영한 데이터의 분산 $\mathbf{w}^\top\mathbf{C}\mathbf{w}$를 최대화하는 단위 벡터 $\mathbf{w}$를 찾는다. Rayleigh quotient(7.5장)에 의해 답은 최대 고유벡터 $\mathbf{v}_1$이다.
• 재구성 오차 최소화. $k$차원 부분공간에 사영했다가 되돌렸을 때의 오차 $\lVert\mathbf{X}_c - \mathbf{X}_c\mathbf{W}\mathbf{W}^\top\rVert_F^2$를 최소화하는 정규직교 $\mathbf{W}$를 찾는다. Eckart–Young에 의해 답은 $\mathbf{V}_k$다.

두 관점이 일치하는 이유는 피타고라스 정리다. 각 점의 제곱 거리는 "부분공간 안의 성분 제곱 + 직교 성분 제곱"으로 나뉘고 그 합은 고정되어 있으므로, 사영된 분산을 최대화하는 것과 남는 오차를 최소화하는 것은 같은 일이다.

rng = np.random.default_rng(1)
N = 500
z = rng.normal(size=(N, 2))                               # (N, 2) 잠재 요인
W = np.array([[2.0, 0.0, 1.0], [0.0, 1.0, 1.0]])          # (2, 3)
D = z @ W + 0.1 * rng.normal(size=(N, 3)) + [1.0, 10.0, -2.0]   # (N, 3)

Dc = D - D.mean(axis=0)                                   # 중심화
U, s, Vt = np.linalg.svd(Dc, full_matrices=False)         # (N, 3), (3,), (3, 3)
var_svd = s ** 2 / (N - 1)
w, Q = np.linalg.eigh(np.cov(D, rowvar=False))            # 오름차순
print("SVD 분산 :", np.round(var_svd, 4))
print("eigh 분산:", np.round(w[::-1], 4))
print("주성분 방향 |cos|:", np.round(np.abs(np.sum(Vt * Q[:, ::-1].T, axis=1)), 6))
print("설명 분산 비율:", np.round(var_svd / var_svd.sum(), 4))
Z = Dc @ Vt.T / np.sqrt(var_svd)                          # (N, 3) whitening
print("whitening 후 공분산 = I ?", np.allclose(np.cov(Z, rowvar=False), np.eye(3)))

SVD 분산 : [4.9071 1.7522 0.0097]
eigh 분산: [4.9071 1.7522 0.0097]
주성분 방향 |cos|: [1. 1. 1.]
설명 분산 비율: [0.7358 0.2627 0.0014]
whitening 후 공분산 = I ? True

잠재 요인 2개로 만든 3차원 데이터에서, 두 주성분이 분산의 약 99.85%를 설명하고 셋째 성분은 잡음(분산 $0.1^2 = 0.01$ 근처)만 남는다. SVD와 공분산 고유분해는 같은 분산과 방향을 준다(고유벡터의 부호는 다를 수 있어 $|\cos|$로 비교했다). SVD를 쓰는 이유는 $\mathbf{X}_c^\top\mathbf{X}_c$를 만들지 않아 조건수를 제곱하지 않기 때문이다(7.4장).

whitening. 점수를 각 성분의 표준편차로 나누면($\mathbf{Z} = \mathbf{X}_c\mathbf{V}\boldsymbol{\Lambda}^{-1/2}$, $\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^2/(N-1)$) 공분산이 단위행렬이 된다. 특징 사이의 상관을 없애고 분산을 1로 맞추는 전처리다. 주의할 점은 분산이 아주 작은 방향(위의 셋째 성분처럼 잡음뿐인 방향)도 1로 늘려 잡음을 크게 키운다는 것이다. 보통 작은 성분은 버리거나 분모에 작은 상수를 더한다.

중심화. PCA의 모든 식은 $\mathbf{X}_c$, 즉 평균을 뺀 데이터에 대한 것이다. 평균을 빼지 않으면 무슨 일이 생기는지 같은 데이터로 본다.

_, s_raw, Vt_raw = np.linalg.svd(D, full_matrices=False)  # 평균을 빼지 않음
mean_dir = D.mean(axis=0) / np.linalg.norm(D.mean(axis=0))   # (3,)
print("|cos(중심화 없는 제1방향, 평균 방향)| =", round(abs(Vt_raw[0] @ mean_dir), 4))
print("|cos(중심화한 제1방향,   평균 방향)| =", round(abs(Vt[0] @ mean_dir), 4))
print("중심화 없이 '분산' 비율:", np.round(s_raw ** 2 / (s_raw ** 2).sum(), 4))

|cos(중심화 없는 제1방향, 평균 방향)| = 1.0
|cos(중심화한 제1방향,   평균 방향)| = 0.0657
중심화 없이 '분산' 비율: [0.9459 0.0437 0.0104]

중심화하지 않은 첫 방향은 데이터의 변동이 아니라 원점에서 데이터 구름의 중심을 가리키는 방향이다. 그리고 그 방향이 "분산"의 95%를 설명하는 것처럼 보인다. 실제 변동의 주축과는 거의 직교한다($|\cos| = 0.066$). 원점에서 먼 데이터일수록 이 착시가 심하다.

계산

비용. $m \times n$ 행렬($m \ge n$)의 전체 SVD는 대략 $mn^2$에 비례하는 연산이 필요하다. 행렬이 수만 × 수만이면 계산도 메모리도 감당하기 어렵다. 그런데 대부분의 응용은 앞의 $k$개 특이값과 특이벡터만 필요하다.

truncated SVD. 앞의 $k$개만 반복법으로 구하는 방법이 있다. 7.5장의 거듭제곱법을 블록으로 확장하거나 Krylov 부분공간 방법(Lanczos 등)을 쓴다. 이런 방법은 행렬–벡터 곱만 필요하므로 희소 행렬이나 명시적으로 저장하지 않은 연산자에도 적용된다. 라이브러리 예로 scipy.sparse.linalg.svds가 있다(집필 시점 예시. 인자와 기본 알고리즘은 사용하는 버전의 문서를 확인한다).

randomized SVD. 무작위 벡터 몇 개를 행렬에 곱하면, 그 결과는 높은 확률로 앞의 특이벡터 방향을 잘 잡아낸다. 이 관찰을 이용하는 절차는 다음과 같다.

1. $n \times (k + p)$ 무작위 행렬 $\boldsymbol{\Omega}$를 만든다($p$는 5–10 정도의 여유분, oversampling).
2. $\mathbf{Y} = \mathbf{A}\boldsymbol{\Omega}$를 계산한다. 특이값이 천천히 감소하면 $\mathbf{Y} = (\mathbf{A}\mathbf{A}^\top)^q\mathbf{A}\boldsymbol{\Omega}$로 거듭제곱을 $q$번 섞어 큰 특이값 방향을 강조한다.
3. $\mathbf{Y}$의 QR 분해로 정규직교 basis $\mathbf{Q}$를 얻는다. $\mathbf{A} \approx \mathbf{Q}\mathbf{Q}^\top\mathbf{A}$이다.
4. 작은 행렬 $\mathbf{B} = \mathbf{Q}^\top\mathbf{A}$($(k + p) \times n$)의 SVD $\mathbf{B} = \tilde{\mathbf{U}}\boldsymbol{\Sigma}\mathbf{V}^\top$를 구하고, $\mathbf{U} = \mathbf{Q}\tilde{\mathbf{U}}$로 되돌린다.

rng = np.random.default_rng(2)
m, n, k, p = 1000, 200, 5, 10
Ua = np.linalg.qr(rng.normal(size=(m, n)))[0]             # (m, n)
Va = np.linalg.qr(rng.normal(size=(n, n)))[0]             # (n, n)
sv = 0.7 ** np.arange(n)                                  # 빠르게 감소하는 특이값
M = (Ua * sv) @ Va.T                                      # (m, n)

Omega = rng.normal(size=(n, k + p))                       # (n, k + p) 무작위 시험 벡터
Y = M @ (M.T @ (M @ Omega))                               # (m, k + p) 거듭제곱 1회
Qr = np.linalg.qr(Y)[0]                                   # (m, k + p) 상의 basis
s_small = np.linalg.svd(Qr.T @ M, compute_uv=False)       # (k + p, n) 작은 행렬의 SVD
print("참 특이값       :", np.round(sv[:k], 6))
print("randomized SVD  :", np.round(s_small[:k], 6))
print(f"최대 상대 오차  : {np.max(np.abs(s_small[:k] - sv[:k]) / sv[:k]):.1e}")

참 특이값       : [1.     0.7    0.49   0.343  0.2401]
randomized SVD  : [1.     0.7    0.49   0.343  0.2401]
최대 상대 오차  : 3.0e-11

큰 SVD는 $15 \times 200$ 행렬의 SVD로 바뀌었다. 이 예는 특이값이 $0.7^i$로 빠르게 감소해 정확도가 매우 높다. 특이값이 천천히 감소하는 행렬에서는 $k + 1$번째 이후 성분이 섞여 들어와 오차가 커지므로, 거듭제곱 횟수 $q$나 여유분 $p$를 늘리고, 정확한 SVD와 비교해 정확도를 확인하는 습관이 필요하다.

활용

SVD와 PCA는 "어떤 방향이 중요하고, 어떤 방향은 버려도 되는가"를 정하는 도구로 거의 모든 분야에서 쓰인다.

• 차원 축소. 주성분 점수 $k$개로 데이터를 요약해 시각화하거나 다음 모델의 입력으로 쓴다. 단, PCA가 고르는 방향은 분산이 큰 방향이지 예측에 유용한 방향이 아니다. 분산이 작은 방향에 타깃과 관련된 신호가 있으면 PCA로 줄이는 순간 사라진다.
• 노이즈 제거. 신호가 저랭크이고 잡음이 모든 방향에 고르게 퍼져 있으면, 작은 특이값을 버리는 것이 잡음을 줄인다. 앞의 장(field) 예에서 rank 3 근사는 잡음 바닥을 버리고 세 패턴만 남긴다. 몇 개를 남길지는 특이값의 바닥 수준과 평가 데이터로 정한다(실습).
• 모델 압축. 가중치 $\mathbf{W} \in \mathbb{R}^{m \times n}$을 $\mathbf{U}_k\boldsymbol{\Sigma}_k\mathbf{V}_k^\top$로 바꾸면 파라미터가 $mn$에서 $k(m + n)$으로 준다. 층 하나를 폭 $k$의 두 층으로 바꾸는 것과 같다. Eckart–Young은 가중치 공간에서의 최적성을 말할 뿐, 모델 출력의 정확도가 가장 덜 떨어진다는 보장은 아니다. 압축 뒤에는 반드시 과제 성능을 다시 잰다.
• 표현 분석. 신경망의 중간 표현(embedding)을 샘플 × 차원 행렬로 모아 특이값 스펙트럼을 보면, 표현이 실제로 몇 개의 방향을 쓰는지(유효 차원) 알 수 있다. 두 모델이나 두 층의 표현을 각각의 주성분 부분공간으로 비교하는 방법도 이 위에 세워진다(18.8장).