The Singular Value Decomposition (SVD) is the pinnacle of matrix factorizations in computational linear algebra. Unlike the eigendecomposition, which is defined only for square matrices and may fail to exist for defective operators, the SVD exists unconditionally for every matrix of arbitrary rectangular dimensions \(\bA \in \fR^{m \times n}\).
Geometrically, the SVD states that the action of any linear transformation on the unit sphere in \(\fR^n\) maps it into a hyperellipsoid in \(\fR^m\). The right singular vectors \(\bv_i\) define the principal orthogonal directions in the domain, the singular values \(\sigma_i\) specify the scaling lengths along each semi-axis, and the left singular vectors \(\bu_i\) define the corresponding orthogonal directions in the codomain.
The SVD Theorem and Subspace Characterization
Let \(\bA \in \fR^{m \times n}\) with \(m \geq n\). There exist orthogonal matrices \(\bU \in \fR^{m \times m}\) and \(\bV \in \fR^{n \times n}\), and a diagonal matrix \(\bsigma \in \fR^{m \times n}\), such that \[
\begin{align}
\bA = \bU \bsigma \bV^T,
\end{align}
\] where the diagonal entries of \(\bsigma = \operatorname{diag}(\sigma_1, \sigma_2, ..., \sigma_n)\) are real, non-negative, and sorted in non-increasing order: \[
\begin{align}
\sigma_1 \geq \sigma_2 \geq ... \geq \sigma_r > \sigma_{r+1} = ... = \sigma_n = 0.
\end{align}
\] The number of strictly positive singular values \(r = \operatorname{rank}(\bA)\) equals the rank of \(\bA\). In dyadic (outer product) form, \[
\begin{align}
\bA = \sum_{i=1}^r \sigma_i \bu_i \bv_i^T.
\end{align}
\]
Let \(\bA \in \fR^{m \times n}\) have rank \(r \leq n \leq m\). The thin (economy) SVD is given by \[
\begin{align}
\bA = \bU_1 \bsigma_1 \bV_1^T,
\end{align}
\] where \(\bU_1 \in \fR^{m \times r}\) and \(\bV_1 \in \fR^{n \times r}\) have orthonormal columns, and \(\bsigma_1 = \operatorname{diag}(\sigma_1, ..., \sigma_r) \in \fR^{r \times r}\) is strictly positive. The columns of the complete singular vectors provide optimal orthonormal bases for the four fundamental subspaces:
Column Space (Range): \(\mathcal{R}(\bA) = \operatorname{span}\{\bu_1, ..., \bu_r\} \subset \fR^m\).
Left Null Space: \(\mathcal{N}(\bA^T) = \operatorname{span}\{\bu_{r+1}, ..., \bu_m\} \subset \fR^m\).
Row Space: \(\mathcal{R}(\bA^T) = \operatorname{span}\{\bv_1, ..., \bv_r\} \subset \fR^n\).
Null Space (Kernel): \(\mathcal{N}(\bA) = \operatorname{span}\{\bv_{r+1}, ..., \bv_n\} \subset \fR^n\).
From \(\bA = \bU\bsigma\bV^T\), we have \(\bA\bv_i = \sigma_i \bu_i\) for \(i = 1, ..., r\), and \(\bA\bv_i = \bzero\) for \(i = r+1, ..., n\). Any vector \(\bx \in \fR^n\) can be expanded in the orthonormal basis \(\{\bv_1, ..., \bv_n\}\) as \(\bx = \sum_{i=1}^n (\bv_i^T\bx)\bv_i\). Applying \(\bA\) gives \[
\begin{align}
\bA\bx = \sum_{i=1}^n (\bv_i^T\bx) \bA\bv_i = \sum_{i=1}^r \sigma_i (\bv_i^T\bx) \bu_i.
\end{align}
\] This expresses \(\bA\bx\) as a linear combination of \(\{\bu_1, ..., \bu_r\}\), proving that \(\mathcal{R}(\bA) = \operatorname{span}\{\bu_1, ..., \bu_r\}\). Furthermore, \(\bA\bx = \bzero\) if and only if \(\bv_i^T\bx = 0\) for all \(i = 1, ..., r\), which proves that \(\mathcal{N}(\bA) = \operatorname{span}\{\bv_{r+1}, ..., \bv_n\}\). Repeating the argument for \(\bA^T = \bV\bsigma^T\bU^T\) completes the proof for the row space and left null space.
Low-Rank Matrix Approximation
One of the most important practical applications of the SVD is optimal low-rank matrix approximation, which forms the computational basis for Principal Component Analysis, latent semantic indexing, and data compression.
Let \(\bA \in \fR^{m \times n}\) have singular value decomposition \(\bA = \sum_{i=1}^r \sigma_i \bu_i \bv_i^T\). For any integer \(k < r\), let the truncated SVD be defined by \[
\begin{align}
\bA_k = \sum_{i=1}^k \sigma_i \bu_i \bv_i^T.
\end{align}
\] Then \(\bA_k\) is the best rank-\(k\) approximation to \(\bA\) in both the spectral norm and the Frobenius norm: \[
\begin{align}
\min_{\operatorname{rank}(\bB) \leq k} \|\bA - \bB\|_2 &= \|\bA - \bA_k\|_2 = \sigma_{k+1}, \\
\min_{\operatorname{rank}(\bB) \leq k} \|\bA - \bB\|_F &= \|\bA - \bA_k\|_F = \sqrt{\sum_{j=k+1}^r \sigma_j^2}.
\end{align}
\]
(Numerical rank and noise filtering) In practical data analysis, singular values rarely equal zero exactly due to measurement noise and roundoff. The numerical rank of \(\bA\) is defined as the number of singular values exceeding a tolerance threshold \(\tau\): \[
\begin{align}
r_{\text{num}} = \max \{ i : \sigma_i > \tau \}, \qquad \tau = \varepsilon_{\text{mach}} \max(m,n) \sigma_1.
\end{align}
\] Truncating the singular value expansion at \(k = r_{\text{num}}\) removes noise components and extracts the dominant low-dimensional dynamical modes.
The Moore-Penrose Pseudoinverse
Let \(\bA \in \fR^{m \times n}\) have thin SVD \(\bA = \bU_1 \bsigma_1 \bV_1^T\). The Moore-Penrose pseudoinverse of \(\bA\) is defined by \[
\begin{align}
\bA^\dagger = \bV_1 \bsigma_1^{-1} \bU_1^T = \sum_{i=1}^r \frac{1}{\sigma_i} \bv_i \bu_i^T.
\end{align}
\] For any right-hand side vector \(\bb \in \fR^m\), the vector \(\hat{\bx} = \bA^\dagger\bb\) is the unique solution that minimizes \(\|\bA\bx - \bb\|_2\), and has minimal Euclidean norm \(\|\hat{\bx}\|_2\) among all least-squares minimizers.
(Step-by-step SVD calculation) Consider the rank-1 matrix \[
\begin{align}
\bA = \begin{pmatrix} 3 & 3 \\ 4 & 4 \end{pmatrix}.
\end{align}
\] We compute \(\bA^T\bA\): \[
\begin{align}
\bA^T\bA = \begin{pmatrix} 3 & 4 \\ 3 & 4 \end{pmatrix} \begin{pmatrix} 3 & 3 \\ 4 & 4 \end{pmatrix} = \begin{pmatrix} 25 & 25 \\ 25 & 25 \end{pmatrix}.
\end{align}
\] The eigenvalues of \(\bA^T\bA\) are \(\lambda_1 = 50\) and \(\lambda_2 = 0\). The singular values are \(\sigma_1 = \sqrt{50} = 5\sqrt{2}\) and \(\sigma_2 = 0\). The normalized eigenvector corresponding to \(\lambda_1 = 50\) is \(\bv_1 = \frac{1}{\sqrt{2}}(1, 1)^T\). The first left singular vector is \[
\begin{align}
\bu_1 = \frac{1}{\sigma_1} \bA\bv_1 = \frac{1}{5\sqrt{2}} \begin{pmatrix} 3 & 3 \\ 4 & 4 \end{pmatrix} \begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix} = \frac{1}{10} \begin{pmatrix} 6 \\ 8 \end{pmatrix} = \begin{pmatrix} 3/5 \\ 4/5 \end{pmatrix}.
\end{align}
\] The thin SVD is \(\bA = \sigma_1 \bu_1 \bv_1^T = (5\sqrt{2}) \begin{pmatrix} 3/5 \\ 4/5 \end{pmatrix} \begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\).
Compute the singular values and singular vectors by hand for \(\bA = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}\).
Let \(\bA \in \fR^{n \times n}\) be invertible. Prove that \(\|\bA\|_2 = \sigma_{\max}\) and \(\|\bA^{-1}\|_2 = 1/\sigma_{\min}\), confirming that \(\kappa_2(\bA) = \sigma_{\max}/\sigma_{\min}\).
Prove that the pseudoinverse satisfies the four Moore-Penrose conditions: (1) \(\bA\bA^\dagger\bA = \bA\), (2) \(\bA^\dagger\bA\bA^\dagger = \bA^\dagger\), (3) \((\bA\bA^\dagger)^T = \bA\bA^\dagger\), and (4) \((\bA^\dagger\bA)^T = \bA^\dagger\bA\).
Load an image as a matrix in Python. Compute its SVD and plot the relative Frobenius approximation error as a function of the truncation rank \(k\).