LINEAR ALGEBRA / 8. SVD & DECOMPOSITIONS

SVD & Matrix Decompositions

The most useful matrix factorization in all of data science

EXPLANATION

Singular Value Decomposition (SVD) decomposes any matrix A into:
A = U Σ Vᵀ

Where:
• U (m×m) → left singular vectors (orthonormal), column space
• Σ (m×n) → diagonal matrix of singular values σ₁ ≥ σ₂ ≥ ... ≥ 0
• Vᵀ (n×n) → right singular vectors (orthonormal), row space

Geometric interpretation: any linear transformation can be broken into:
1. Rotation (Vᵀ)
2. Scaling (Σ)
3. Another rotation (U)

SVD applications in ML:
• PCA: eigenvectors of XᵀX = right singular vectors of X
• Low-rank approximation: keep only top-k singular values
• Recommender systems: matrix factorization
• Solving overdetermined systems: pseudoinverse A⁺ = V Σ⁺ Uᵀ
• Measuring matrix stability: condition number = σ_max/σ_min
• Compressing neural network weights: LoRA fine-tuning uses low-rank decomposition!

Truncated SVD (rank-k approximation):
A ≈ UₖΣₖVₖᵀ (keep only k largest singular values)
This gives the best rank-k approximation in terms of Frobenius norm.

DIAGRAM

SVD: A(m×n) = U(m×m) Σ(m×n) Vᵀ(n×n)

  For A(3×2):
  U(3×3) · Σ(3×2) · Vᵀ(2×2)
  [rotates] [scales] [rotates]

  Σ = [[σ₁, 0 ],
       [0,  σ₂],   σ₁ ≥ σ₂ ≥ 0
       [0,  0 ]]

  Low-rank approximation (rank-1):
  A ≈ σ₁ · u₁ · v₁ᵀ   (outer product)
  Each additional term adds more detail

  LoRA (Low-Rank Adaptation) for LLMs:
  ΔW = A · B  where A:(d×r), B:(r×d), r << d
  Only fine-tune the low-rank factors A and B
  Instead of updating the full weight matrix W

CODE

PYTHON

1import numpy as np
2from sklearn.decomposition import TruncatedSVD
3
4# ── SVD decomposition ─────────────────────────────────────────────
5A = np.array([[1., 2., 3.],
6              [4., 5., 6.],
7              [7., 8., 9.],
8              [10.,11.,12.]])   # (4, 3)
9
10U, S, Vt = np.linalg.svd(A, full_matrices=True)
11print(f"SVD shapes: U={U.shape}, S={S.shape}, Vt={Vt.shape}")
12print(f"Singular values: {S.round(4)}")
13
14# Reconstruct A from SVD
15A_reconstructed = U[:, :len(S)] @ np.diag(S) @ Vt[:len(S), :]
16print(f"Reconstruction error: {np.max(abs(A - A_reconstructed)):.2e}")
17
18# ── Low-rank approximation ────────────────────────────────────────
19np.random.seed(42)
20A_data = np.random.randn(50, 30)   # 50 samples, 30 features
21
22U, S, Vt = np.linalg.svd(A_data, full_matrices=False)
23
24print("
25Low-rank approximations:")
26for k in [1, 3, 5, 10, 30]:
27    A_k = U[:,:k] @ np.diag(S[:k]) @ Vt[:k,:]
28    error = np.linalg.norm(A_data - A_k, "fro")
29    variance = (S[:k]**2).sum() / (S**2).sum()
30    print(f"  rank-{k:2d}: error={error:.3f}, variance explained={variance:.3f}")
31
32# ── SVD for PCA ───────────────────────────────────────────────────
33# PCA via SVD (numerically more stable than eigendecomposition)
34np.random.seed(42)
35X = np.random.randn(100, 5)
36X = X - X.mean(axis=0)   # center
37
38U, S, Vt = np.linalg.svd(X, full_matrices=False)
39# Principal components = rows of Vt
40# Singular values² / (n-1) = eigenvalues of covariance matrix
41print(f"
42PCA via SVD:")
43print(f"  Principal component 1: {Vt[0].round(3)}")
44var_explained = S**2 / (S**2).sum()
45for i, v in enumerate(var_explained):
46    print(f"  PC{i+1}: {v*100:.1f}% variance")
47
48# ── Pseudoinverse (Moore-Penrose) ─────────────────────────────────
49# A⁺ = V Σ⁺ Uᵀ  — solves overdetermined systems in least-squares sense
50A_rect = np.random.randn(10, 4)   # more equations than unknowns
51b = np.random.randn(10)
52x_lstsq = np.linalg.pinv(A_rect) @ b   # least-squares solution
53x_np    = np.linalg.lstsq(A_rect, b, rcond=None)[0]
54print(f"
55Pseudoinverse vs lstsq match: {np.allclose(x_lstsq, x_np)}")

← PREV7. Dot Product & Norms NEXT →9. PCA from Scratch