Matrix Decompositions

If you were to learn only one matrix decomposition technique for machine learning, Singular Value Decomposition (SVD) would be the unambiguous choice. It is, without exaggeration, the most versatile and widely-used matrix factorization in applied mathematics and data science.

SVD reveals the fundamental structure of any matrix—not just square matrices with special properties, but any matrix of any shape. It decomposes a matrix into three simpler matrices that expose the underlying geometry: the directions of maximum variance, the magnitudes of stretching along those directions, and how inputs and outputs relate through the transformation.

From Netflix's recommendation engine that powered its famous million-dollar prize competition to Google's latent semantic indexing for search, from image compression algorithms to noise reduction in signal processing, SVD is the mathematical workhorse operating behind countless real-world systems.

Let's begin with the precise mathematical statement. Given any matrix A ∈ ℝ^(m×n), there exist orthogonal matrices U ∈ ℝ^(m×m) and V ∈ ℝ^(n×n), and a diagonal matrix Σ ∈ ℝ^(m×n) such that:

$$A = U \Sigma V^T$$

This factorization is called the Singular Value Decomposition of A.

Let's unpack each component:

The matrix U (Left Singular Vectors):

• U is an m × m orthogonal matrix (U^T U = UU^T = I_m)
• The columns of U are called the left singular vectors of A
• These vectors form an orthonormal basis for the column space (range) of A
• They represent the output directions of the transformation

The matrix Σ (Singular Values):

• Σ is an m × n diagonal matrix (zero everywhere except the main diagonal)
• The diagonal entries σ₁, σ₂, ..., σᵣ are called the singular values of A
• By convention: σ₁ ≥ σ₂ ≥ ... ≥ σᵣ > 0 (ordered from largest to smallest)
• r is the rank of matrix A
• The singular values are always real and non-negative

The matrix V (Right Singular Vectors):

• V is an n × n orthogonal matrix (V^T V = VV^T = I_n)
• The columns of V are called the right singular vectors of A
• These vectors form an orthonormal basis for the row space of A
• They represent the input directions of the transformation

Understanding the dimensions:

For a matrix A of size m × n with rank r:

• U has size m × m (but only the first r columns are 'active')
• Σ has size m × n with r non-zero diagonal entries
• V has size n × n (but only the first r columns are 'active')

This leads to different forms of SVD:

Full SVD: U is m × m, Σ is m × n, V is n × n

Economy (Thin) SVD: U is m × r, Σ is r × r, V is n × r

The economy SVD is more storage-efficient when r << min(m, n) and is what most software packages return by default.

The power of SVD lies not just in its existence but in what it tells us about the geometry of linear transformations. To understand this, let's think about what happens when we multiply a vector x by the matrix A.

The transformation A = UΣV^T as three steps:

When we compute Ax, the SVD decomposes this into three sequential operations:

Step 1: V^T rotates (and possibly reflects) the input space The matrix V^T rotates the input vector x from the standard basis to a new coordinate system defined by the right singular vectors. This aligns the vector with the 'principal directions' of the transformation.

Step 2: Σ scales along each axis The diagonal matrix Σ stretches or compresses along each axis by the corresponding singular value σᵢ. Directions with larger singular values get amplified more; directions with σᵢ = 0 are collapsed entirely (these correspond to the null space).

Step 3: U rotates into the output space Finally, U rotates the scaled result into the output coordinate system, aligning it with the left singular vectors.

The sphere-to-ellipse transformation:

Perhaps the most illuminating visualization: consider what happens when A acts on the unit sphere. If we take all vectors x with ||x|| = 1 (a sphere in n dimensions), and apply A to each, we get an ellipsoid in m dimensions.

• The axes of the ellipsoid are aligned with the left singular vectors (columns of U)
• The lengths of the semi-axes are the singular values σ₁, σ₂, ..., σᵣ
• The orientation in input space that produces each axis is given by the corresponding right singular vectors (columns of V)

Connecting singular vectors to fundamental subspaces:

SVD provides orthonormal bases for all four fundamental subspaces of a matrix A with rank r:

This elegant connection shows why SVD is considered the 'ultimate' matrix decomposition—it reveals the complete structure of any linear transformation.

Interpretation of singular values:

• σ₁ is the maximum amount the matrix can stretch any unit vector: σ₁ = max ||Ax|| for ||x|| = 1
• The ratio σ₁/σᵣ (condition number) measures how 'ill-conditioned' the matrix is
• σᵢ represents the 'importance' of the i-th direction—larger values capture more of the matrix's action
• The sum Σσᵢ² equals ||A||²_F (the squared Frobenius norm)
• The number of non-zero singular values equals the rank of A

Understanding how SVD is computed illuminates its deep connection to eigenvalue problems. While production implementations use sophisticated numerical methods (like the Golub-Kahan algorithm), the conceptual derivation reveals essential insights.

Deriving SVD from eigenvalues:

The key insight is that if we have A = UΣV^T, then:

$$A^T A = V \Sigma^T U^T U \Sigma V^T = V \Sigma^T \Sigma V^T = V \Lambda_V V^T$$

$$A A^T = U \Sigma V^T V \Sigma^T U^T = U \Sigma \Sigma^T U^T = U \Lambda_U U^T$$

where Λ_V = Σ^TΣ and Λ_U = ΣΣ^T are diagonal matrices containing σ₁², σ₂², ...

This reveals that:

• V contains the eigenvectors of A^TA (an n × n symmetric matrix)
• U contains the eigenvectors of AA^T (an m × m symmetric matrix)
• The singular values σᵢ are the square roots of the eigenvalues of A^TA (or equivalently, AA^T)

Since A^TA and AA^T are symmetric positive semi-definite, their eigenvalues are real and non-negative. This guarantees that singular values are always real and non-negative.

The computational algorithm (conceptual):

1. Compute B = A^TA (if n ≤ m) or B = AA^T (if m < n)
2. Find the eigenvalues λ₁ ≥ λ₂ ≥ ... ≥ λᵣ > 0 and eigenvectors of B
3. Set σᵢ = √λᵢ for the singular values
4. The eigenvectors of A^TA give V; the eigenvectors of AA^T give U
5. Alternatively: once V is known, compute uᵢ = (1/σᵢ) A vᵢ to get U

Computational complexity:

For an m × n matrix (assume m ≥ n):

• Full SVD: O(mn² + n³) operations
• Truncated SVD (k singular values): O(mnk) using iterative methods
• Memory: O(mn) for the input, plus O(mk + nk + k) for truncated output

This makes SVD relatively expensive for large matrices, motivating the use of approximate and randomized algorithms for big data applications.

Uniqueness of SVD:

SVD is essentially unique with the following caveats:

• If all singular values are distinct, the singular vectors are unique up to sign
• If some singular values are equal (multiplicity > 1), the corresponding singular vectors are only unique up to rotation within their eigenspace
• The singular values themselves are always uniquely determined

One of the most important applications of SVD is finding the best low-rank approximation to a matrix. This is fundamental to dimensionality reduction, data compression, noise reduction, and matrix completion.

The problem statement:

Given a matrix A of rank r, find the rank-k matrix Aₖ (where k < r) that is closest to A. 'Closest' means minimizing either:

• The Frobenius norm: ||A - Aₖ||_F = √(Σᵢⱼ (A - Aₖ)ᵢⱼ²)
• The spectral norm: ||A - Aₖ||₂ = largest singular value of (A - Aₖ)

The Eckart-Young-Mirsky Theorem (1936):

The optimal rank-k approximation Aₖ is obtained by truncating the SVD:

$$A_k = \sum_{i=1}^{k} \sigma_i u_i v_i^T = U_k \Sigma_k V_k^T$$

where Uₖ contains the first k columns of U, Σₖ is the k × k diagonal matrix of the top k singular values, and Vₖ contains the first k columns of V.

The approximation errors are:

• ||A - Aₖ||₂ = σₖ₊₁ (the next singular value)
• ||A - Aₖ||_F = √(σₖ₊₁² + σₖ₊₂² + ... + σᵣ²)

This is a profound result: among all rank-k matrices, the truncated SVD is provably optimal in both common matrix norms!

Practical implications:

1. Data Compression: Instead of storing m × n numbers, store:

• Uₖ: m × k numbers
• σ₁...σₖ: k numbers
• Vₖ: n × k numbers

Total: k(m + n + 1) numbers. If k << min(m,n), this is massive compression.

2. Noise Reduction: If data = signal + noise, and the signal is effectively low-rank while noise is full-rank, truncating small singular values removes noise while preserving signal.

3. Choosing k:

• Explained variance ratio: keep enough singular values to explain a target percentage (e.g., 95%) of the total variance (sum of σᵢ²)
• Elbow method: plot singular values and look for a sharp drop
• Cross-validation: choose k that minimizes held-out reconstruction error

Outer product form of SVD:

The SVD can be written as a sum of rank-1 matrices: $$A = \sum_{i=1}^{r} \sigma_i u_i v_i^T$$

Each term σᵢuᵢvᵢᵀ is a rank-1 'building block'. The truncated SVD simply takes the first k of these blocks.

SVD and eigendecomposition are related but distinct factorizations. Understanding when to use each is crucial for ML practitioners.

Eigendecomposition recap:

For a square matrix A ∈ ℝⁿˣⁿ, the eigendecomposition is: $$A = Q \Lambda Q^{-1}$$

where Λ is diagonal containing eigenvalues, and Q contains eigenvectors.

For symmetric matrices A = Aᵀ, this simplifies to: $$A = Q \Lambda Q^T$$

where Q is orthogonal (Q⁻¹ = Qᵀ).

Key differences:

Relationship for symmetric matrices:

For a symmetric positive semi-definite matrix A:

• Eigenvalues = Singular values (both non-negative)
• Eigenvectors = Left singular vectors = Right singular vectors
• Both decompositions are essentially the same

For a symmetric matrix with negative eigenvalues:

• Singular values = |eigenvalues|
• The signs are 'absorbed' into the singular vectors

When to use SVD:

1. Rectangular matrices — eigendecomposition doesn't apply
2. Low-rank approximation — SVD gives optimal rank-k approximation
3. Solving least squares — SVD handles rank-deficient cases gracefully
4. Computing pseudoinverse — A⁺ = VΣ⁺Uᵀ
5. Principal Component Analysis — SVD of data matrix is more numerically stable than eigendecomposition of covariance matrix
6. Matrix completion/recommendation systems — natural fit for user-item matrices

When to use eigendecomposition:

1. Understanding dynamics — eigenvalues determine system behavior (stability, oscillation, decay)
2. Graph problems — graph Laplacian eigenvalues reveal connectivity, spectral clustering
3. Markov chains — stationary distribution is the eigenvector for eigenvalue 1
4. Quadratic forms — diagonalization simplifies xᵀAx
5. Matrix powers — Aⁿ = QΛⁿQ⁻¹ is efficient for large n
6. Matrix exponential — exp(A) = Q exp(Λ) Q⁻¹

SVD is not just a theoretical tool—it's the computational backbone of many ML systems. Let's explore the major applications with practical details.

1. Principal Component Analysis (PCA):

PCA finds the orthogonal directions of maximum variance in data. Given centered data matrix X (n samples × d features), compute SVD: X = UΣVᵀ

• Principal components (directions): columns of V
• Projections onto k dimensions: XVₖ = UₖΣₖ
• Explained variance ratios: σᵢ² / Σⱼσⱼ²

2. Latent Semantic Analysis (LSA):

For text analysis, create a term-document matrix A where Aᵢⱼ = frequency of term i in document j. SVD reveals:

• Latent topics as singular vectors
• Topic strengths as singular values
• Documents and terms mapped to shared latent space

3. Recommendation Systems:

The Netflix Prize made matrix factorization famous. For user-item rating matrix R:

• Factor into user features U and item features V: R ≈ UVᵀ
• SVD provides initialization or direct solution
• Predict missing ratings by reconstructing from low-rank factors

4. Image Compression:

Treat an image as a matrix (or separately for each color channel). Truncated SVD gives lossy compression: storing k(m+n+1) values instead of mn, with provably minimal reconstruction error.

5. Pseudoinverse and Least Squares:

The Moore-Penrose pseudoinverse A⁺ is elegantly computed via SVD: $$A^+ = V \Sigma^+ U^T$$

where Σ⁺ has 1/σᵢ on the diagonal for non-zero singular values (and 0 elsewhere).

This solves least squares problems even when A is rank-deficient: $$x = A^+ b$$

6. Noise Reduction:

If measurements are corrupted by noise, and the true signal is approximately low-rank:

1. Compute SVD of noisy data
2. Threshold or truncate small singular values
3. Reconstruct from remaining components

The small singular values typically correspond to noise, while large ones capture the signal.

7. Matrix Completion:

Given a matrix with missing entries (e.g., user-movie ratings), find a low-rank matrix that matches the observed entries. Nuclear norm minimization (sum of singular values) is a convex relaxation that often recovers the true matrix. SVD is central to algorithms like Soft-Impute.

Computing full SVD on large matrices is expensive. Understanding computational trade-offs is essential for practical ML applications.

Complexity Analysis:

Truncated SVD:

When you only need the top k singular values/vectors (typical in ML), use iterative methods that never compute the full decomposition:

• scipy.sparse.linalg.svds for sparse matrices
• sklearn.decomposition.TruncatedSVD for general use

Randomized SVD:

For very large matrices, randomized algorithms provide excellent approximations much faster:

1. Create a random projection matrix Ω (n × k or slightly larger)
2. Form Y = AΩ (captures the range of A)
3. Orthonormalize Y to get Q
4. Form the smaller matrix B = QᵀA
5. Compute SVD of B
6. Recover U from the SVD of B

This trades exact computation for probabilistic guarantees with much lower cost.

We've covered SVD in depth—from its formal definition through geometric interpretation, computation, and practical applications. Let's consolidate the key insights.

What's next:

Having mastered SVD, we'll explore Eigendecomposition in depth—understanding when and why to decompose square matrices by their eigenvalues and eigenvectors. While SVD is the general workhorse, eigendecomposition is essential for understanding matrix dynamics, spectral graph theory, and many probabilistic models.