15 Overdetermined Systems and Least Squares
Approximating the solution to an overdetermined system where \(m > n\) and no exact solution exists.
When \(\bb\) is not in the column space of \(\bA\), the equation \(\bA\bx=\bb\) has no exact solution. Least squares replaces the impossible goal make the residual zero'' with the geometric goalmake the residual as short as possible.’’
15.1 Least Squares Problem
\(\bA\bx = \bb\) with \(\bA \in \fR^{m \times n}, m > n\). Usually \(\bb \notin \mathcal{R}(\bA)\).
Find \(\hat{\bx}\) that minimizes the Euclidean residual: \[ \begin{align} \hat{\bx} = \arg\min_{\bx \in \fR^n} \|\bA\bx - \bb\|_2^2. \end{align} \] Geometry: \(\bA\hat{\bx}\) is the orthogonal projection of \(\bb\) onto \(\mathcal{R}(\bA)\).
For the least squares problem in the result above, the result above applied to \(S=\mathcal{R}(\bA)\) implies that the solution \(\hat{\bx}\) is characterized by \[ \begin{align} \bA\hat{\bx} = \operatorname{proj}_{\mathcal{R}(\bA)}(\bb), \qquad \br = \bb-\bA\hat{\bx} \perp \mathcal{R}(\bA). \end{align} \] Equivalently, \(\bA^T\br=\bzero\).
(Fitting a line) For data \((t_i,y_i)\), fitting \(y\approx c_0+c_1t\) gives \[ \begin{align} \bA = \begin{pmatrix} 1&t_1\\ 1&t_2\\ \vdots&\vdots\\ 1&t_m \end{pmatrix}, \qquad \bx=\begin{pmatrix}c_0\\c_1\end{pmatrix}, \qquad \bb=\begin{pmatrix}y_1\\y_2\\ \vdots\\ y_m\end{pmatrix}. \end{align} \] The residuals are the vertical errors between the data and the fitted line.
15.2 Normal Equations
Using the orthogonality condition in the result above, the minimizer satisfies \(\bA^T\bA\hat{\bx} = \bA^T\bb\).
Condition: \(\bA^T\br = \bzero\) (Residual must be normal to the column space).
Unique Solution: \(\hat{\bx} = (\bA^T\bA)^{-1}\bA^T\bb\) if \(\bA\) has full column rank.
Starting from the result above, derive the normal equations in the result above.
Show that \(\bA^T\br=\bzero\) is equivalent to \(\br\in\mathcal{N}(\bA^T)\).
Use the result above to explain why \(\br\in\mathcal{N}(\bA^T)\) means \(\bA\hat{\bx}\) is the projection of \(\bb\) onto \(\mathcal{R}(\bA)\).
Fit a line \(y = mx + b\) to points \((1,2), (2,3), (3,5), (4,4)\) via normal equations.
Prove \(\kappa_2(\bA^T\bA) = \kappa_2(\bA)^2\) using singular values. If \(\kappa_2(\bA)=10^6\), estimate the number of digits that may be lost in double precision.
Explain why the normal equations are useful for derivation but should not be the default numerical algorithm. Your answer should mention \(\bA^T\br=\bzero\), conditioning, QR, and SVD.
15.3 QR-Based Least Squares
The numerically superior approach.
If \(\bA = \bQ\bR\) (thin QR), the least squares solution satisfies: \[ \begin{align} \bR\hat{\bx} = \bQ^T\bb. \end{align} \] This avoids forming \(\bA^T\bA\), the instability diagnosed in the result above.
The goal is to minimize \(\|\bA\bx - \bb\|_2^2\). Let \(\bA = \bQ\bR\) be the thin QR factorization. Using the orthogonal projector \(\bP = \bQ\bQ^T\) and its complement \(\bI - \bP\): \[ \begin{align} \|\bA\bx - \bb\|_2^2 &= \|\bQ\bR\bx - \bb\|_2^2 \\ &= \|\bQ\bR\bx - (\bQ\bQ^T\bb + (\bI - \bQ\bQ^T)\bb)\|_2^2 \\ &= \|(\bQ\bR\bx - \bQ\bQ^T\bb) - (\bI - \bQ\bQ^T)\bb\|_2^2. \end{align} \] Since \(\bQ(\bR\bx - \bQ^T\bb) \in \mathcal{R}(\bQ)\) and \((\bI - \bQ\bQ^T)\bb \in \mathcal{R}(\bQ)^\perp\), the Pythagorean theorem gives: \[ \begin{align} \|\bA\bx - \bb\|_2^2 = \|\bQ(\bR\bx - \bQ^T\bb)\|_2^2 + \|(\bI - \bQ\bQ^T)\bb\|_2^2. \end{align} \] The first term is minimized (set to zero) when \(\bR\bx = \bQ^T\bb\), and the second term is independent of \(\bx\). Thus \(\hat{\bx} = \bR^{-1}\bQ^T\bb\) is the unique minimizer.
(Least squares solver hierarchy)
Normal equations: Fast and simple, but least stable; condition number is squared.
QR: Standard method for full-rank least squares; stable and efficient.
SVD: Most robust; handles rank deficiency and reveals numerical rank.
15.4 Rank Deficiency and Regularization
If \(\bA\) does not have full column rank, there may be infinitely many least squares minimizers. The SVD selects the minimizer with smallest Euclidean norm. Regularization modifies the problem to penalize large coefficients and improve conditioning.
For any \(\bA\), the vector \[ \begin{align} \hat{\bx}=\bA^\dagger\bb \end{align} \] is the least squares minimizer with smallest Euclidean norm.
Write \(\bA=\bU\bsigma\bV^T\).
Transform the least squares problem into SVD coordinates.
Identify which coordinates are determined by nonzero singular values.
Show that setting nullspace coordinates to zero gives the minimum-norm solution in the result above.
For \(\lambda>0\), ridge regression solves \[ \begin{align} \hat{\bx}_\lambda = \argmin_{\bx}\left(\|\bA\bx-\bb\|_2^2+\lambda\|\bx\|_2^2\right). \end{align} \] The normal equations become \[ \begin{align} (\bA^T\bA+\lambda\bI)\hat{\bx}_\lambda = \bA^T\bb. \end{align} \]
Given positive weights \(w_i\), weighted least squares solves \[ \begin{align} \min_{\bx}\|\mathbf{W}^{1/2}(\bA\bx-\bb)\|_2^2, \qquad \mathbf{W}=\operatorname{diag}(w_1,...,w_m). \end{align} \] Large \(w_i\) forces the fit to match observation \(i\) more strongly.
The weighted least squares problem in the result above satisfies \[ \begin{align} \bA^T\mathbf{W}\bA\hat{\bx}=\bA^T\mathbf{W}\bb. \end{align} \]
Apply the ordinary normal equations from the result above to the transformed system \[ \begin{align} \mathbf{W}^{1/2}\bA\bx\approx \mathbf{W}^{1/2}\bb. \end{align} \] This gives \((\mathbf{W}^{1/2}\bA)^T(\mathbf{W}^{1/2}\bA)\hat{\bx}=(\mathbf{W}^{1/2}\bA)^T\mathbf{W}^{1/2}\bb\), which simplifies to the stated equation.
Return to the result above. Solve the same line-fitting problem using QR. Compare coefficients, residual norms, and the conditioning of the two linear systems solved.
Regularization: If \(\bA\) is rank-deficient, \(\|\hat{\bx}\|\) can explode. Add a penalty term \(\lambda \|\bx\|_2^2\) (Ridge Regression). How do the normal equations change?
Derive the result above by differentiating the weighted objective directly.
15.5 Polynomial Regression
For nodes \(t_1, ..., t_m\), the degree-\(n\) design matrix is: \[ \begin{align} V_{ij} = t_i^{j-1}. \end{align} \] Polynomial fitting \(p(t) = \sum c_j t^j\) is solved as \(\bV\bc \approx \mathbf{f}\).
Fit a degree-2 polynomial to 6 data points. Verify the residual and \(R^2\) score.
Numerical Stability: Construct a Vandermonde matrix with \(n=10\). Compare the errors of
lstsq(A, b)vssolve(A.T @ A, A.T @ b). Relate the outcome to the result above.