Least Squares, Pseudoinverses, and Low-Rank Approximation
~6 min read · Tier toggle changes depth (T1 plain · T2 standard · T3 rigour)
You should be comfortable with the Singular Value Decomposition (SVD) and its properties, as well as basic matrix algebra and the concept of orthogonal projections in inner-product spaces.
Why this matters
Real-world data is noisy, and models are often overdetermined or ill-posed.
In experimental physics, we rarely have a system of equations with a unique exact solution. Instead, we have more measurements than unknowns, or the system is sensitive to small errors. This lecture introduces the tools to find the "best" approximate solution in the sense of minimizing error, and to compress data by retaining only the most significant features. These concepts are fundamental to signal processing, image compression, and quantum state tomography.
The Geometry of Approximation
Imagine you have a set of data points that do not lie perfectly on a straight line. You want to find the line that is "closest" to the data. In vector space terms, you are trying to project a vector (your data) onto a subspace (the set of all possible lines). The best approximation is the one where the difference between the actual data and your model is perpendicular to the model space. This is the core idea of least squares: we minimize the length of the error vector. When the system is underdetermined (more unknowns than equations), we look for the solution with the smallest magnitude. The SVD provides the natural coordinate system to solve these problems efficiently and robustly.
Least Squares via SVD
Consider a linear system \(A\mathbf{x} = \mathbf{b}\), where \(A\) is an \(m \times n\) matrix with \(m > n\) (overdetermined). If \(A\) has full column rank, the columns are linearly independent, but \(\mathbf{b}\) may not lie in the column space of \(A\). We seek \(\mathbf{x}\) that minimizes \(\|\mathbf{b} - A\mathbf{x}\|^2\). Using the SVD \(A = U\Sigma V^{\dagger}\), we can transform the problem. Let \(\mathbf{y} = V^{\dagger}\mathbf{x}\) and \(\mathbf{c} = U^{\dagger}\mathbf{b}\). The norm becomes \(\|\mathbf{c} - \Sigma\mathbf{y}\|^2\). Since \(\Sigma\) is diagonal, this separates into independent scalar minimizations. For each singular value \(\sigma_i > 0\), the optimal choice is \(y_i = c_i / \sigma_i\). For \(\sigma_i = 0\), the term is independent of \(y_i\), so we set \(y_i = 0\) for the minimum norm solution. This leads to the formula \(\mathbf{x} = V\Sigma^{+}U^{\dagger}\mathbf{b}\).
The Moore-Penrose Pseudoinverse
The matrix \(A^{+} = V\Sigma^{+}U^{\dagger}\) is called the Moore-Penrose pseudoinverse. It generalizes the inverse to non-square or singular matrices. The matrix \(\Sigma^{+}\) is an \(n \times m\) matrix where the non-zero diagonal entries are the reciprocals of the singular values, and all other entries are zero. This construction ensures that \(A^{+}\) provides the unique solution to the least squares problem that also has the minimum Euclidean norm. If the system is consistent (\(\mathbf{b}\) is in the column space of \(A\)), \(A^{+}\mathbf{b}\) is the exact solution with the smallest norm. If inconsistent, it is the best approximate solution. The pseudoinverse is crucial for handling rank-deficient systems common in physics, such as determining parameters from redundant measurements.
Low-Rank Approximation
Many physical systems have structure that can be captured by a few dominant modes. The Eckart-Young theorem states that the best rank-\(k\) approximation to a matrix \(A\) in the Frobenius norm is obtained by keeping only the \(k\) largest singular values. Let \(A_k = U_k \Sigma_k V_k^{\dagger}\), where \(U_k\) and \(V_k\) contain the first \(k\) columns of \(U\) and \(V\), and \(\Sigma_k\) is the \(k \times k\) diagonal matrix of the largest singular values. The error is \(\|A - A_k\|_F = \sqrt{\sum_{i=k+1}^{p} \sigma_i^2}\). This is the mathematical basis for Principal Component Analysis (PCA) and image compression. In quantum mechanics, it relates to truncating the Hilbert space to the most relevant states, ignoring high-energy or low-probability components.
The derivation
Singular Value Decomposition
Let \(A:\mathcal{V}\to\mathcal{W}\) be a linear map between finite-dimensional inner-product spaces over \(\mathbb{C}\) (or \(\mathbb{R}\)), with \(\dim\mathcal{V}=n\) and \(\dim\mathcal{W}=m\). Then there exist an orthonormal basis \(\{v_i\}\) of \(\mathcal{V}\), an orthonormal basis \(\{u_j\}\) of \(\mathcal{W}\), and real numbers \(\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_p\ge 0\) with \(p=\min(m,n)\), such that in these bases \(A=U\Sigma V^{\dagger}\), where \(U\) and \(V\) are unitary and \(\Sigma\in\mathbb{R}^{m\times n}\) is diagonal with non-negative entries \(\Sigma_{ii}=\sigma_i\).
When singular values are very small, dividing by them in the pseudoinverse amplifies noise. In practice, one uses "truncated SVD" or Tikhonov regularization, setting \(\sigma_i^{-1} = 0\) for \(\sigma_i < \epsilon\) to stabilize the solution. This is a form of regularization that trades bias for variance. C
Worked example
Problem. A sensor array gives the measurements \(\mathbf{b} = (1, 2, 4)^{\mathrm T}\) for the model \(A\mathbf{x} = \mathbf{b}\) with \(A = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}\). The system is overdetermined and inconsistent. Find the least-squares solution and the residual norm.
- Singular values: \(A^{\mathrm T}A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) has eigenvalues 3 and 1, so \(\sigma_1 = \sqrt 3\) and \(\sigma_2 = 1\). Both are non-zero, so \(A\) has full column rank and the least-squares solution is unique.
- Normal equations \(A^{\mathrm T}A\,\mathbf{x} = A^{\mathrm T}\mathbf{b}\): \(A^{\mathrm T}\mathbf{b} = (1 + 4,\ 2 + 4)^{\mathrm T} = (5, 6)^{\mathrm T}\).
- Solve: \((A^{\mathrm T}A)^{-1} = \frac13\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}\), so \(\mathbf{x} = \frac13 (2\cdot 5 - 6,\ -5 + 2\cdot 6)^{\mathrm T} = \left(\tfrac43,\ \tfrac73\right)^{\mathrm T}\). (The SVD formula \(\mathbf{x} = V\Sigma^{+}U^{\dagger}\mathbf{b}\) gives the same vector.)
- Residual: \(A\mathbf{x} = \left(\tfrac43, \tfrac73, \tfrac{11}{3}\right)^{\mathrm T}\), so \(\mathbf{r} = \mathbf{b} - A\mathbf{x} = \left(-\tfrac13, -\tfrac13, \tfrac13\right)^{\mathrm T}\) and \(\|\mathbf{r}\| = \sqrt{3/9} = 1/\sqrt 3 \approx 0.577\). Check: \(A^{\mathrm T}\mathbf{r} = (0, 0)^{\mathrm T}\), so the residual is orthogonal to the columns of \(A\), as it must be.
Units check. If \(\mathbf{b}\) is in volts and \(A\) is dimensionless, then \(\mathbf{x}\), \(\mathbf{r}\) and \(\|\mathbf{r}\|\) are all in volts.
Three quick tests
- What is the rank of the pseudoinverse \(A^{+}\) if \(A\) has rank \(r\)?
Answer
It is also \(r\). The pseudoinverse preserves the rank of the original matrix. - If all singular values of \(A\) are equal to 1, what is \(A^{+}\)?
Answer
Then \(\Sigma^{+} = \Sigma^{\mathrm T}\), so \(A^{+} = V\Sigma^{\mathrm T}U^{\dagger} = (U\Sigma V^{\dagger})^{\dagger} = A^{\dagger}\), the conjugate transpose (\(A^{\mathrm T}\) for a real matrix). If \(A\) is square, it is unitary and \(A^{+} = A^{-1} = A^{\dagger}\). - Why is the minimum norm solution important in underdetermined systems?
Answer
It provides a unique solution among infinitely many, selecting the one with the smallest energy or magnitude, which is often physically preferred (e.g., minimum energy state).
Problem set for PU-104 →
Each derivation in this unit ends with five graded problems with full solutions.
Lecture notes drafted with a local model and checked line by line by the editor; every number in the worked example was recomputed.