Computing a Singular Value Decomposition

Factor a small matrix into rotation, scaling, and rotation, and read off rank, range, and the best low-rank approximation.

What the SVD is

Any real matrix A (m by n) factors as A = U S V^T, where U and V are orthogonal and S is diagonal with non-negative entries, the singular values, in decreasing order. Geometrically A maps the unit sphere to an ellipsoid; the singular values are the semi-axis lengths and U gives their directions.

By hand via A^T A

The singular values are the square roots of the eigenvalues of A^T A; the right singular vectors V are its eigenvectors; then U columns come from u_i = A v_i / s_i. For large matrices this is numerically poor, but it shows where the pieces come from.

import numpy as np
A=np.array([[3.0,0.0],[4.0,5.0]])
U,s,Vt=np.linalg.svd(A)
print('singular values:',np.round(s,4))
print('reconstruct:',np.round(U@np.diag(s)@Vt,4))
# rank-1 best approximation
A1=s[0]*np.outer(U[:,0],Vt[0])
print('rank-1:',np.round(A1,4))

What it tells you

The number of nonzero singular values is the rank. Their ratio (largest over smallest) is the condition number, warning of near-singularity. The columns of U for nonzero singular values span the range of A. Small singular values flag directions the matrix nearly collapses.

Low-rank approximation

Keeping only the top k singular triplets gives the best rank-k approximation in both the spectral and Frobenius norms (the Eckart-Young theorem). This single fact powers image compression, principal component analysis, latent-semantic indexing, and noise reduction - throw away the small singular values and you keep the essential structure.