Singular Value Decomposition (SVD)

The Singular Value Decomposition (SVD) is arguably the most important result in applied linear algebra. While eigendecomposition works only for special square matrices, SVD works for any matrix—tall, wide, square, singular, or non-singular. It provides a way to “see” the underlying structure of data by decomposing it into its most significant components.

The Geometric Idea

Every matrix $A$ represents a linear map. SVD says that any such map can be broken down into three simple steps:

1. A rotation in the input space ($V$).
2. A scaling along the principal axes ($\Sigma$).
3. A rotation in the output space ($U$).

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

• $V$: Columns are “right singular vectors.” They define an orthonormal basis in the input space.
• $\Sigma$: A diagonal matrix of “singular values” $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_n \ge 0$. These tell you the “strength” or “gain” of the matrix in each direction.
• $U$: Columns are “left singular vectors.” They define an orthonormal basis in the output space.

Data Compression: The Best Low-Rank Approximation

The real magic of SVD is the Eckart-Young Theorem. It states that if you want the best possible “summary” of a matrix $A$ using only $k$ dimensions (where $k$ is less than the rank of $A$), the answer is to keep only the $k$ largest singular values and their corresponding vectors.

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

This is how image compression and noise reduction work. By throwing away small singular values, we lose “noise” or “unimportant detail” but keep the overall structure.

import numpy as np

# Create a 'rank-2' matrix with some noise
# This matrix has 100 rows (data points) and 5 columns (features)
np.random.seed(42)
base_data = np.random.randn(100, 2) @ np.random.randn(2, 5)
noise = 0.1 * np.random.randn(100, 5)
A = base_data + noise

# Perform SVD
U, s, Vt = np.linalg.svd(A, full_matrices=False)

print("Singular Values:", s)

# Construct Rank-2 approximation
k = 2
Sk = np.diag(s[:k])
Ak = U[:, :k] @ Sk @ Vt[:k, :]

# Calculate Error
full_norm = np.linalg.norm(A, 'fro')
error_norm = np.linalg.norm(A - Ak, 'fro')

print(f"\nRelative Error of Rank-{k} approximation: {error_norm/full_norm:.4f}")
print(f"Original Rank: {np.linalg.matrix_rank(A)}")
print(f"Approximation Rank: {np.linalg.matrix_rank(Ak)}")

The Pseudoinverse: Solving the Solvable and Unsolvable

When a matrix $A$ is not invertible (e.g., it is not square or is singular), we can still “solve” $Ax = b$ using the Moore-Penrose Pseudoinverse $A^\dagger$.

Using SVD, the pseudoinverse is trivial to compute: $A^\dagger = V \Sigma^\dagger U^T$ where $\Sigma^\dagger$ is formed by transposing $\Sigma$ and replacing every non-zero $\sigma_i$ with $1/\sigma_i$.

The solution $x = A^\dagger b$ is the “best” solution in two senses:

1. It minimizes the error $\|Ax - b\|^2$ (Least Squares).
2. If there are many such solutions, it picks the one with the smallest length $\|x\|$.

Exercises