SVD — 4TB Memory Crash on Sparse Matrices

SVD is simultaneously the most powerful and most practically useful decomposition in linear algebra. Unlike eigendecomposition (restricted to square diagonalisable matrices), SVD works on any m×n matrix. It has a geometric interpretation (rotation-scale-rotation), a dimensionality reduction application (truncated SVD = PCA), and a computational application (pseudoinverse, rank, condition number).

Netflix's original collaborative filtering used SVD to decompose the user-movie rating matrix. Image compression, latent semantic analysis (LSA) for text, and principal component analysis all reduce to truncated SVD. NumPy's np.linalg.svd and scipy's truncated_svd (for large sparse matrices) implement it.

SVD Decomposition

A = U Σ V^T where U (m×m) and V (n×n) are orthogonal matrices, Σ (m×n) is diagonal with non-negative singular values σ₁ ≥ σ₂ ≥ ... ≥ 0.

import numpy as np

A = np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], dtype=float)
U, s, Vt = np.linalg.svd(A, full_matrices=False)

print(f'A shape: {A.shape}')
print(f'U: {U.shape}, s: {s.shape}, Vt: {Vt.shape}')
print(f'Singular values: {s}')
print(f'Rank (non-zero singular values): {np.sum(s > 1e-10)}')

# Reconstruct A from SVD
A_reconstructed = U @ np.diag(s) @ Vt
print(f'Reconstruction error: {np.max(np.abs(A - A_reconstructed)):.2e}')

Low-Rank Approximation — Truncated SVD

Keep only the k largest singular values — the rank-k approximation is the closest rank-k matrix to A (in Frobenius norm).

def truncated_svd(A, k):
    """Best rank-k approximation of A."""
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    return U[:,:k] @ np.diag(s[:k]) @ Vt[:k,:]

# Image compression example
import numpy as np
np.random.seed(42)
image = np.random.rand(100, 100)  # 100x100 'image'
for k in [1, 5, 10, 50]:
    approx = truncated_svd(image, k)
    error  = np.linalg.norm(image - approx, 'fro') / np.linalg.norm(image, 'fro')
    compression = k * (100 + 100 + 1) / (100*100)
    print(f'k={k:3}: error={error:.3f}, storage={compression:.2%}')

Geometric Intuition of SVD: Rotation, Scale, Rotation

Any linear transformation (matrix) applied to a vector can be broken down into a rotation (V^T), a scaling (Σ) along the new axes, and another rotation (U). The columns of V are the right singular vectors (input directions), columns of U are the left singular vectors (output directions). In 2D, a unit circle maps to an ellipse whose semi-axes lengths are the singular values.

import numpy as np
import matplotlib.pyplot as plt

A = np.array([[3, 1], [1, 2]])
U, s, Vt = np.linalg.svd(A)

# Circle of points
angles = np.linspace(0, 2*np.pi, 100)
circle = np.array([np.cos(angles), np.sin(angles)])

# Transform circle
ellipse = A @ circle

plt.figure(figsize=(6,6))
plt.plot(circle[0], circle[1], 'b--', label='Unit circle')
plt.plot(ellipse[0], ellipse[1], 'r-', label='Transformed')
# Plot axes
plt.arrow(0,0, s[0]*U[0,0], s[0]*U[1,0], head_width=0.2, color='g')
plt.arrow(0,0, s[1]*U[0,1], s[1]*U[1,1], head_width=0.2, color='g')
plt.axis('equal')
plt.legend()
plt.show()

Numerical Computation: Algorithms and Stability

Computing SVD is numerically stable but expensive. Full SVD uses QR decomposition (Golub-Kahan algorithm) in O(mn·min(m,n)). For large matrices, use randomized SVD: project the matrix onto a random low-dimensional subspace, then compute QR and small SVD. Power iteration can improve accuracy. scipy.linalg.svd is the workhorse; for sparse matrices use scipy.sparse.linalg.svds (ARPACK).

import numpy as np
from sklearn.utils.extmath import randomized_svd

A = np.random.rand(1000, 100)
k = 10
U, s, Vt = randomized_svd(A, n_components=k, n_iter='auto', random_state=42)
print(f'Computed top {k} singular values in ~0.1s vs 2s for full SVD')

Applications: PCA, Pseudoinverse, and Recommender Systems

PCA: center the data, compute SVD of centered matrix, left singular vectors multiplied by singular values give principal components. Pseudoinverse A⁺ = V Σ⁺ U^T (reciprocal of non-zero singular values). Used for solving least-squares and underdetermined systems. In recommender systems: user-item matrix R ≈ U_k Σ_k V_k^T; predict missing ratings via inner product. This is the FunkSVD / SVD++ family.

import numpy as np
X = np.random.randn(100, 5)
X_centered = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(X_centered, full_matrices=False)
# Principal components are columns of Vt.T
components = Vt.T
# Project onto first 2 PCs
X_pca = X_centered @ components[:,:2]
print(f'Explained variance ratio: {s[:2]**2 / np.sum(s**2)}')