Singular Value Decomposition

The Singular Value Decomposition states that every matrix $A \in \mathbb{R}^{m\times n}$ can be written as

where $U \in \mathbb{R}^{m\times m}$ is orthogonal ($U^\top U = I$), $V \in \mathbb{R}^{n\times n}$ is orthogonal ($V^\top V = I$), and $\Sigma \in \mathbb{R}^{m\times n}$ is diagonal with nonnegative entries $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_{\min(m,n)} \geq 0$ on its main diagonal, called singular values. The columns of $U$ are left singular vectors, and the columns of $V$ are right singular vectors. Geometrically, $A$ acts on a vector by first rotating it (via $V^\top$), then stretching it along coordinate axes by the factors $\sigma_i$ (via $\Sigma$), then rotating again (via $U$). This decomposition always exists — no assumption of square shape, invertibility, or diagonalizability is needed — which is precisely why SVD is more broadly applicable than eigendecomposition. Eigendecomposition $A = P D P^{-1}$ requires $A$ to be square and diagonalizable, and $P$ need not be orthogonal. The SVD sidesteps all of this by working with two different orthonormal bases (one for the domain, one for the codomain) instead of one shared basis. The connection to eigenvalues is direct. Since $U, V$ are orthogonal,

which is an eigendecomposition of the symmetric positive semi-definite matrix $A^\top A$: the columns of $V$ are its eigenvectors, and the diagonal entries of $\Sigma^\top \Sigma$ (i.e., $\sigma_i^2$) are its eigenvalues. Symmetrically, $A A^\top = U (\Sigma \Sigma^\top) U^\top$ gives $U$ as the eigenvectors of $AA^\top$. This is exactly how SVD connects to Principal Component Analysis: if $A$ is a mean-centered $n \times p$ data matrix (rows are samples), then $A^\top A / (n-1)$ is the sample covariance matrix, whose eigenvectors are the principal component directions — and those eigenvectors are precisely the right singular vectors of $A$.

Step 1 — reduce to an eigenproblem. Conceptually (not how numerical libraries do it internally, but mathematically equivalent), form $A^\top A$, an $n \times n$ symmetric positive semi-definite matrix. Its eigenvalues $\lambda_1 \geq \dots \geq \lambda_n \geq 0$ and orthonormal eigenvectors $v_1, \dots, v_n$ give the right singular vectors directly, with $\sigma_i = \sqrt{\lambda_i}$.

Step 2 — recover the left singular vectors. For each nonzero singular value, set

These $u_i$ are automatically orthonormal because $u_i^\top u_j = \frac{1}{\sigma_i \sigma_j} v_i^\top A^\top A v_j = \frac{\lambda_j}{\sigma_i \sigma_j} v_i^\top v_j$, which is $1$ when $i=j$ and $0$ otherwise (using $v_i^\top v_j = \delta_{ij}$). If $m > n$ (or there are zero singular values), extend $\{u_i\}$ to a full orthonormal basis of $\mathbb{R}^m$ via Gram–Schmidt on any complementary vectors; these extra columns multiply zero rows of $\Sigma$ and do not affect $A$.

Step 3 — assemble. Stack the $v_i$ as columns of $V$, the $u_i$ as columns of $U$, and place $\sigma_i$ on the diagonal of $\Sigma$ in decreasing order. Then $A = U\Sigma V^\top$ holds by construction, since applying both sides to each $v_i$ gives $A v_i = \sigma_i u_i$ on both sides, and the $v_i$ form a basis. In practice, numerical libraries never explicitly form $A^\top A$ (this squares the condition number and loses precision); instead they use bidiagonalization followed by an implicit QR or divide-and-conquer iteration directly on $A$. The pseudocode below presents the conceptual algorithm (correct and pedagogically transparent) alongside the key downstream use — truncated SVD for low-rank approximation. ```text function SVD(A): # A is m x n

# Conceptual computation (production code uses bidiagonal QR / divide&conquer) C := transpose(A) @ A # n x n, symmetric PSD (eigenvalues λ, eigenvectors v) := EIGENDECOMPOSE(C) sort λ descending; permute v to match for i in 1..n: σ[i] := sqrt(max(λ[i], 0)) # numerical floor at 0 for i in 1..rank: # rank = number of nonzero σ[i] u[i] := (A @ v[i]) / σ[i] extend {u[i]} to an orthonormal basis of R^m via Gram-Schmidt (fills U) V := columns v[1..n] U := columns u[1..m] Σ := m x n diagonal matrix with entries σ[1..min(m,n)] return U, Σ, V function TRUNCATED_SVD(A, k): U, Σ, V := SVD(A) U_k := first k columns of U Σ_k := top-left k x k block of Σ V_k := first k columns of V return U_k, Σ_k, V_k # A_k = U_k Σ_k V_k^T is the best rank-k approx function LOW_RANK_RECONSTRUCT(U_k, Σ_k, V_k): return U_k @ Σ_k @ transpose(V_k) # explicit m x n matrix, rank <= k ``` The truncated form is the workhorse in applications: instead of storing $A_k$ explicitly, one stores the three thin factors $U_k$ ($m\times k$), $\Sigma_k$ ($k\times k$), $V_k$ ($n \times k$), which is far cheaper when $k \ll \min(m,n)$.

Time complexity. Computing a full SVD of a dense $m \times n$ matrix (assume $m \geq n$ without loss of generality) via the standard approach — Householder bidiagonalization followed by an implicit-shift QR or divide-and-conquer algorithm on the bidiagonal form — costs

for the singular values and $V$ alone, and $O(mn^2 + n^3)$ (dominated by the same term when $m \geq n$) if $U$ is also required in full. More generally this is stated as $O(mn \cdot \min(m,n))$, matching $O(mn^2)$ when $n \leq m$. This dominates because bidiagonalization is a sequence of $O(n)$ Householder reflections, each costing $O(mn)$, and the subsequent bidiagonal SVD iteration converges in $O(n)$ sweeps each costing $O(n)$ per superdiagonal entry, giving lower-order terms. Forming $A^\top A$ explicitly (the conceptual algorithm above) costs $O(mn^2)$ just for the product and then $O(n^3)$ for its eigendecomposition — asymptotically comparable but numerically inferior, since squaring $A$ squares its condition number and can destroy precision in small singular values. Production libraries (LAPACK's `gesdd`/`gesvd`) avoid this by working with $A$ directly.

Space complexity. Storing the full factorization needs $O(m^2 + mn + n^2)$ in the worst case (full $U$ and $V$), but truncated SVD needs only $O(mk + k + nk) = O(k(m+n))$ for rank $k$, which is the practically relevant figure since $k \ll \min(m,n)$ in most applications.

Randomized/truncated speedups. When only the top $k$ singular vectors are needed and $A$ is large or sparse, randomized SVD algorithms (Halko–Martinsson–Tropp) compute an approximate truncated SVD in $O(mn\log k + (m+n)k^2)$ time by projecting onto a random $k$-dimensional subspace, orthogonalizing, and solving a small dense SVD on the projected matrix. This is the practical choice for matrices with millions of rows/columns, such as recommender-system rating matrices.

Correctness — existence. The existence of the SVD follows from the eigendecomposition of $A^\top A$: it is symmetric and positive semi-definite (since $x^\top A^\top A x = \|Ax\|^2 \geq 0$), so the spectral theorem guarantees real nonnegative eigenvalues and an orthonormal eigenbasis. The construction in Step 2 above (defining $u_i = Av_i/\sigma_i$) is then verified to yield orthonormal vectors satisfying $Av_i = \sigma_i u_i$ for every $i$, which is exactly the statement $AV = U\Sigma$, i.e. $A = U\Sigma V^\top$ since $V$ is orthogonal ($V^{-1}=V^\top$).

Correctness — optimality (Eckart–Young theorem). The deepest correctness claim is that truncating the SVD gives the provably best low-rank approximation. Formally, for any matrix $B$ with $\text{rank}(B) \leq k$,

where $A_k = U_k \Sigma_k V_k^\top$ is the truncated SVD. Moreover the exact approximation error is known in closed form:

Proof sketch. Write $A$ in its SVD basis; since Frobenius norm is invariant under orthogonal transformations, $\|A - A_k\|_F^2 = \|\Sigma - \Sigma_k\|_F^2 = \sum_{i>k} \sigma_i^2$ directly, giving the error formula for the truncated SVD. For optimality, suppose $B$ has rank $\leq k$; its null space has dimension $\geq n-k$, so it intersects the span of the top $k+1$ right singular vectors $\{v_1,\dots,v_{k+1}\}$ (a $(k+1)$-dimensional space) in a nonzero vector $z$ (by a dimension-counting argument: two subspaces of dimensions $\geq n-k$ and $k+1$ inside $\mathbb{R}^n$ must intersect nontrivially). Normalize $\|z\|=1$; then $Bz = 0$, so

using that $\sigma_i \geq \sigma_{k+1}$ for $i \leq k+1$ and $\sum (v_i^\top z)^2 = \|z\|^2 = 1$. Hence no rank-$k$ matrix beats the truncated SVD in spectral norm, and a parallel argument (or the Frobenius-norm case of Mirsky's theorem) extends this to Frobenius norm. This is the mathematical fact that makes SVD the gold standard for compression, denoising, and dimensionality reduction: it is not merely a good heuristic but a certified optimum under two of the most natural matrix norms. $\blacksquare$

Tradeoffs. SVD assumes linear (subspace) structure; when data lies on a curved manifold, a linear rank-$k$ approximation can badly misrepresent geometry — this is exactly the gap that Manifold Learning techniques address. SVD is also sensitive to scale: features measured in different units will have their singular values dominated by scale rather than information content, so standardization (as in PCA) matters. Numerically, $\sigma_i$ near zero are ill-determined and their corresponding $u_i, v_i$ are unstable to perturbation; truncation before this regime, guided by the spectrum, is standard practice.

Consider the small matrix

Step 1: form $A^\top A$.

Step 2: eigenvalues of $A^\top A$. The characteristic equation is $(25-\lambda)^2 - 400 = 0$, so $25-\lambda = \pm 20$, giving $\lambda_1 = 45$, $\lambda_2 = 5$. Both positive, confirming positive semi-definiteness.

Step 3: singular values.

Step 4: eigenvectors of $A^\top A$ (right singular vectors). For $\lambda_1 = 45$: solve $(25-45)x_1 + 20 x_2 = 0 \Rightarrow -20 x_1 + 20 x_2 = 0 \Rightarrow x_1 = x_2$. Normalized: $v_1 = \frac{1}{\sqrt2}(1,1)^\top$. For $\lambda_2 = 5$: $(25-5)x_1 + 20x_2 = 0 \Rightarrow 20x_1 + 20x_2 = 0 \Rightarrow x_1 = -x_2$. Normalized: $v_2 = \frac{1}{\sqrt2}(1,-1)^\top$.

Step 5: left singular vectors via $u_i = Av_i/\sigma_i$.

Check norm: $\frac{1}{10}(1+9) = 1$. ✓.

Check norm: $\frac{1}{10}(9+1)=1$. ✓. Check orthogonality $u_1^\top u_2 = \frac{1}{10}(1\cdot3 + 3\cdot(-1)) = 0$. ✓.

Step 6: assemble and verify.

A direct multiplication $U\Sigma V^\top$ reconstructs $A$ exactly (verified entrywise: e.g., entry $(1,1)$ is $\frac{1}{\sqrt{10}}\cdot 3\sqrt5 \cdot \frac{1}{\sqrt2} + \frac{3}{\sqrt{10}}\cdot\sqrt5\cdot\frac{1}{\sqrt2} = \frac{3\sqrt5+3\sqrt5}{\sqrt{20}} = \frac{6\sqrt5}{2\sqrt5} = 3$, matching $A_{11}=3$).

Rank-1 truncation. Keeping only $\sigma_1$:

The Eckart–Young theorem guarantees this is the closest rank-1 matrix to $A$ in Frobenius norm, with error exactly $\sqrt{\sigma_2^2} = \sqrt5 \approx 2.236$, which one can verify directly: $\|A - A_1\|_F^2 = (3-1.5)^2+(0-1.5)^2+(4-4.5)^2+(5-4.5)^2 = 2.25+2.25+0.25+0.25 = 5$, and indeed $\sigma_2^2 = 5$. ✓ This confirms the closed-form error formula on a concrete numeric case, not just symbolically.

The SVD's development spans nearly a century of differential geometry and numerical analysis, culminating in the algorithm that made it computationally practical for real matrices at scale.

Recommender systems. Matrix factorization methods for collaborative filtering — famously used in the Netflix Prize competition — are directly inspired by truncated SVD: user and item latent factor matrices approximate a sparse ratings matrix as a low-rank product, letting the system predict unseen ratings.

Latent Semantic Analysis (LSA) / information retrieval. Applying SVD to term-document matrices reveals latent semantic dimensions, letting search engines and document-similarity systems match on meaning rather than exact keyword overlap.

Image and signal compression. Treating an image as a matrix of pixel intensities, truncating its SVD to a small rank $k$ reconstructs a visually similar image using far fewer numbers than the original ($k(m+n+1)$ versus $mn$), illustrating the Eckart–Young guarantee directly.

Noise reduction and data denoising. Because noise tends to spread energy across many small singular values while true signal concentrates in the top few, truncating the SVD often removes noise while preserving structure — used in genomics (gene expression matrices), sensor arrays, and finance (covariance matrix denoising).

Numerical linear algebra infrastructure. The Moore–Penrose pseudo-inverse, used to solve least-squares problems and ill-conditioned linear systems, is defined directly via the SVD: $A^+ = V \Sigma^+ U^\top$, where $\Sigma^+$ inverts the nonzero singular values. This underlies ridge regression solvers, `numpy.linalg.pinv`, and robust regression pipelines throughout scientific computing.

Foundation for PCA and beyond. Every standard PCA implementation (`sklearn.decomposition.PCA`, R's `prcomp`) computes its result via SVD of the centered data matrix rather than by eigendecomposing the covariance matrix directly, for the same numerical-stability reasons discussed above. This also sets the stage for Manifold Learning methods, which relax SVD's linear-subspace assumption to model curved, nonlinear data geometry.

Ask the tutor: How does randomized SVD achieve near-optimal low-rank approximation without computing the full decomposition?