Singular Value Decomposition (SVD)

Eigen-decomposition works only for square matrices, and only some of them. The singular value decomposition works for every matrix, of any shape:

A = U Σ Vᵀ

For an m × n matrix A of rank r:

• U is an m × m orthogonal matrix; its columns are the left singular vectors.
• V is an n × n orthogonal matrix; its columns are the right singular vectors.
• Σ is an m × n "diagonal" matrix whose diagonal entries σ₁ ≥ σ₂ ≥ … ≥ σᵣ > 0 are the singular values; all other entries are 0.

Geometrically, every linear map is a rotation (Vᵀ), followed by a stretch along the axes (Σ), followed by another rotation (U).

Connection with eigenvalues

From A = UΣVᵀ:

• AᵀA = V (ΣᵀΣ) Vᵀ: the columns of V are eigenvectors of AᵀA, with eigenvalues σᵢ².
• AAᵀ = U (ΣΣᵀ) Uᵀ: the columns of U are eigenvectors of AAᵀ, with the same non-zero eigenvalues σᵢ².

So the singular values are the square roots of the eigenvalues of AᵀA (which are never negative, because AᵀA is positive semidefinite). And Avᵢ = σᵢuᵢ: each right singular vector is mapped to a stretched left singular vector.

For a symmetric positive semidefinite matrix, the SVD and the eigen-decomposition coincide. For a general symmetric matrix, the singular values are the absolute values of the eigenvalues.

Computing an SVD by hand

Example. A = [[3, 0], [4, 5]].

1. AᵀA = [[3, 4], [0, 5]] · [[3, 0], [4, 5]] = [[25, 20], [20, 25]].
2. Its eigenvalues: trace 50, determinant 625 − 400 = 225, so λ² − 50λ + 225 = 0 gives λ = 45 and λ = 5.
3. The singular values are σ₁ = √45 = 3√5 ≈ 6.708 and σ₂ = √5 ≈ 2.236.
4. The eigenvectors of AᵀA give V: v₁ = (1, 1)/√2 for 45 and v₂ = (1, −1)/√2 for 5.
5. Then uᵢ = Avᵢ / σᵢ: u₁ = (3, 9)/(√2 · 3√5) = (1, 3)/√10 and u₂ = (3, −1)/(√2 · √5) = (3, −1)/√10.

Checks: σ₁σ₂ = √225 = 15 = |det A|, and σ₁² + σ₂² = 50 = the sum of the squares of all entries of A (9 + 0 + 16 + 25).

What the singular values tell you

• rank(A) = the number of non-zero singular values.
• The Frobenius norm ‖A‖_F = √(Σ aᵢⱼ²) = √(σ₁² + … + σᵣ²).
• The spectral norm ‖A‖₂ = max ‖Ax‖ over unit vectors x = σ₁.
• For a square matrix, |det(A)| = σ₁σ₂⋯σₙ.
• The condition number σ₁/σₙ measures how sensitive Ax = b is to errors; a huge value means A is nearly singular.

Low-rank approximation

Writing the SVD as a sum of rank-1 pieces:

A = σ₁u₁v₁ᵀ + σ₂u₂v₂ᵀ + … + σᵣuᵣvᵣᵀ

Keeping only the first k terms gives Aₖ. The Eckart-Young theorem: Aₖ is the best rank-k approximation of A, in both the Frobenius and spectral norms, and the errors are

‖A − Aₖ‖₂ = σₖ₊₁ and ‖A − Aₖ‖_F = √(σₖ₊₁² + … + σᵣ²).

This is the basis of image compression, recommender systems (latent factors), noise reduction, and PCA: the principal components of a centred data matrix X are its right singular vectors, and the variance along the i-th component is σᵢ²/(n − 1).

The pseudo-inverse

The Moore-Penrose pseudo-inverse is A⁺ = V Σ⁺ Uᵀ, where Σ⁺ inverts each non-zero singular value and transposes the shape. For any b, x = A⁺b is the least-squares solution of minimum norm of Ax ≈ b. When A has independent columns, A⁺ = (AᵀA)⁻¹Aᵀ, exactly the matrix in the normal equations of linear regression.

import math

A = [[3, 0], [4, 5]]
AtA = [[sum(A[k][i] * A[k][j] for k in range(2)) for j in range(2)] for i in range(2)]
print("A^T A =", AtA)
tr, det = AtA[0][0] + AtA[1][1], AtA[0][0] * AtA[1][1] - AtA[0][1] * AtA[1][0]
l1 = (tr + math.sqrt(tr * tr - 4 * det)) / 2
l2 = (tr - math.sqrt(tr * tr - 4 * det)) / 2
s1, s2 = math.sqrt(l1), math.sqrt(l2)
print(f"eigenvalues of A^T A: {l1:g}, {l2:g}; singular values: {s1:.4f}, {s2:.4f}")

r = 1 / math.sqrt(2)
V = [[r, r], [r, -r]]                    # columns v1, v2
mul = lambda M, v: [sum(M[i][k] * v[k] for k in range(2)) for i in range(2)]
u1 = [x / s1 for x in mul(A, [V[0][0], V[1][0]])]
u2 = [x / s2 for x in mul(A, [V[0][1], V[1][1]])]
print(f"u1 = ({u1[0]:.4f}, {u1[1]:.4f}), u2 = ({u2[0]:.4f}, {u2[1]:.4f}), orthogonal: {abs(u1[0] * u2[0] + u1[1] * u2[1]) < 1e-12}")

# Rebuild A from the two rank-1 pieces, and the best rank-1 approximation
A1 = [[s1 * u1[i] * V[j][0] for j in range(2)] for i in range(2)]
A2 = [[A1[i][j] + s2 * u2[i] * V[j][1] for j in range(2)] for i in range(2)]
print("sigma1 u1 v1^T =", [[round(x, 4) + 0.0 for x in row] for row in A1])
print("sum of both pieces =", [[round(x, 4) + 0.0 for x in row] for row in A2])
err = math.sqrt(sum((A[i][j] - A1[i][j]) ** 2 for i in range(2) for j in range(2)))
print(f"Frobenius error of the rank-1 approximation = {err:.4f} = sigma2")

A^T A = [[25, 20], [20, 25]]
eigenvalues of A^T A: 45, 5; singular values: 6.7082, 2.2361
u1 = (0.3162, 0.9487), u2 = (0.9487, -0.3162), orthogonal: True
sigma1 u1 v1^T = [[1.5, 1.5], [4.5, 4.5]]
sum of both pieces = [[3.0, 0.0], [4.0, 5.0]]
Frobenius error of the rank-1 approximation = 2.2361 = sigma2

Practice questions