Mathematical foundations of Singular Value Decomposition.
The Singular Value Decomposition (SVD)—originating in the independent mathematical investigations of Eugenio Beltrami (1873), Camille Jordan (1874), and James Joseph Sylvester (1889), and elevated to an algorithmic triumph by Gene Golub and William Kahan in 1965—is universally acknowledged as the absolute summit of matrix computations. In linear algebra, eigendecomposition A = VΛV⁻¹ is restricted exclusively to square, diagonalizable matrices, and can become severely ill-conditioned if eigenvectors are nearly collinear. SVD transcends all such limitations: every matrix—square or rectangular, full-rank or rank-deficient, well-conditioned or pathological—admits an SVD: A = UΣVᵀ. Here, U is an m×m orthogonal matrix containing the left singular vectors, V is an n×n orthogonal matrix containing the right singular vectors, and Σ is an m×n diagonal matrix containing the non-negative singular values sorted in descending order: σ₁ ≥ σ₂ ≥ ... ≥ σ_p ≥ 0. Geometrically, SVD decomposes any linear transformation into three elementary geometric actions: an initial rotation in the domain (Vᵀ), an axis-aligned stretching by factors σ_i (Σ), and a final rotation in the codomain (U). In computational science, SVD provides the definitive definition of matrix condition: the 2-norm condition number is given exactly by κ₂(A) = σ_max / σ_min. Furthermore, SVD resolves the universal problem of inverting non-invertible matrices via the Moore-Penrose pseudoinverse: A⁺ = V Σ⁺ Uᵀ. For any linear system Ax = b, x_{LS} = A⁺ b produces the unique solution that simultaneously minimizes the residual error ‖Ax − b‖₂ and has the minimum Euclidean norm ‖x‖₂. Finally, by the Eckart-Young-Mirsky theorem, truncating the SVD to its top k components yields the mathematically optimal low-rank matrix approximation, providing the engine for Principal Component Analysis (PCA), latent semantic indexing, and big-data denoising. In this computational laboratory, students evaluate SVD on a matrix with singular values σ = [10.0, 4.0, 0.4] (condition number κ = 10.0/0.4 = 25.000), construct the Moore-Penrose pseudoinverse, explore rank truncation, and extract the condition number κ₂(A) = 25.000 via linear regression.
Viva voce questions & analytical theory
What is the Singular Value Decomposition (SVD) and why is it called the 'rolling pin of linear algebra'?
The Singular Value Decomposition factors any real m×n matrix A into A = UΣVᵀ, where U is an m×m orthogonal matrix (left singular vectors), V is an n×n orthogonal matrix (right singular vectors), and Σ is an m×n diagonal matrix containing non-negative singular values σ₁ ≥ σ₂ ≥ ... ≥ σ_p ≥ 0 (p = min(m, n)). Geometrically, SVD reveals that every linear transformation maps the unit sphere in ℝⁿ into a hyperellipsoid in ℝᵐ, where the singular values σ_i are the lengths of the semi-axes and the left singular vectors u_i are their spatial orientations. SVD is universally applicable: unlike eigendecomposition, it exists for every matrix, square or rectangular, full-rank or singular.
What is the Moore-Penrose Pseudoinverse A⁺ and what four axioms define it?
Formulated by E. H. Moore (1920) and Arne Bjerhammar/Roger Penrose (1955), the pseudoinverse A⁺ is the unique n×m matrix satisfying the four Moore-Penrose conditions: (1) AA⁺A = A; (2) A⁺AA⁺ = A⁺; (3) (AA⁺)ᵀ = AA⁺; and (4) (A⁺A)ᵀ = A⁺A. Using SVD, the pseudoinverse is computed effortlessly as A⁺ = V Σ⁺ Uᵀ, where Σ⁺ inverts all non-zero singular values (σ_i⁺ = 1/σ_i) and leaves zeros untouched. If A is square and invertible, A⁺ = A⁻¹; if A has full column rank, A⁺ = (AᵀA)⁻¹ Aᵀ.
What is the Eckart-Young-Mirsky Theorem and how does it power data compression (PCA)?
The Eckart-Young-Mirsky theorem (1936) states that for any rank k < r = rank(A), the optimal rank-k approximation to A in both the Frobenius norm and spectral 2-norm is obtained by truncating the SVD to its top k terms: A_k = ∑_{i=1}^k σ_i u_i v_iᵀ. The approximation error is precisely ‖A − A_k‖₂ = σ_{k+1}. This theorem forms the mathematical foundation of Principal Component Analysis (PCA), latent semantic analysis in NLP, image compression, and facial recognition (Eigenfaces).
Why is solving least squares via SVD or QR vastly superior to the Normal Equations (AᵀA x = Aᵀb)?
The normal equations AᵀA x = Aᵀb require forming the product AᵀA, which squares the condition number: κ₂(AᵀA) = [κ₂(A)]². If A has condition number κ = 10⁸, AᵀA has condition number κ² = 10¹⁶, completely exhausting 64-bit IEEE 754 precision (ε_mach ≈ 10⁻¹⁶) and rendering the system insoluble! In contrast, SVD and QR solve least squares operating directly on A without forming AᵀA, preserving condition number κ₂(A) and maintaining full 8-digit accuracy.
How does Truncated SVD (TSVD) achieve regularization in ill-posed inverse problems?
In ill-posed problems (such as medical tomography deconvolution), small singular values σ_i ≈ 0 are dominated by experimental noise. In the pseudoinverse, dividing by tiny σ_i (1/σ_i >> 10¹⁰) blows up noise, destroying the reconstruction. Truncated SVD regularizes the problem by discarding all singular values below a noise threshold ε_tol (setting 1/σ_i = 0), suppressing high-frequency noise amplification while preserving the dominant signal.
What is the relationship between the singular values of A and the eigenvalues of AᵀA and AAᵀ?
Because A = UΣVᵀ, AᵀA = (VΣᵀUᵀ)(UΣVᵀ) = V(ΣᵀΣ)Vᵀ. This is the exact orthogonal eigendecomposition of the symmetric positive semi-definite matrix AᵀA! Thus, the right singular vectors v_i are eigenvectors of AᵀA, the left singular vectors u_i are eigenvectors of AAᵀ, and the singular values σ_i are the non-negative square roots of the eigenvalues of AᵀA: σ_i = √(λ_i(AᵀA)).
How does SVD directly identify the four fundamental subspaces of linear algebra?
Gilbert Strang's four fundamental subspaces are partitioned with exact precision by SVD: (1) Column space Col(A): spanned by the first r left singular vectors {u₁, ..., u_r}; (2) Left nullspace N(Aᵀ): spanned by the remaining {u_{r+1}, ..., u_m}; (3) Row space Row(A): spanned by the first r right singular vectors {v₁, ..., v_r}; and (4) Nullspace N(A): spanned by the remaining {v_{r+1}, ..., v_n}.
Algorithm precautions and stability rules
- Never form the normal equations matrix AᵀA explicitly when A is ill-conditioned; always use SVD or QR.
- Choose singular value truncation threshold ε_tol carefully (e.g. ε_tol = max(m, n) · ε_mach · σ₁) to prevent noise amplification.
- Verify that singular values are sorted in non-increasing order: σ₁ ≥ σ₂ ≥ ... ≥ σ_p ≥ 0.
- Do not invert zero or near-zero singular values; set σ_i⁺ = 0 for σ_i < ε_tol.
- Maintain 64-bit precision; Golub-Kahan bidiagonalization requires dual Householder sweeps.
- Check that U and V are strictly orthogonal: ‖UᵀU − I‖ < 10⁻¹⁴ and ‖VᵀV − I‖ < 10⁻¹⁴.
Theoretical limitations
- SVD is computationally expensive: full SVD requires approximately 4mn² + 8n³ FLOPs.
- Requires storing both orthogonal matrices U (m×m) and V (n×n) unless thin/economy SVD is used.
- Truncating singular values introduces regularization bias that trades accuracy for stability.