In scientific data fitting and parameter estimation, one frequently encounters linear systems \(\bA\bx \approx \bb\) with more equations than unknowns (\(m > n\)). Because the data vector \(\bb \in \fR^m\) is corrupted by measurement noise or modeling discrepancies, it rarely lies precisely within the column space \(\mathcal{R}(\bA) \subset \fR^m\). Consequently, no exact solution exists.
The method of linear least squares replaces the impossible requirement of zero residual with the optimization problem of finding the vector \(\hat{\bx} \in \fR^n\) that minimizes the Euclidean norm of the residual vector \(\br(\bx) = \bb - \bA\bx\). Geometrically, the optimal model prediction \(\bA\hat{\bx}\) is the orthogonal projection of the observation vector \(\bb\) onto the range of \(\bA\).
The Least-Squares Problem and the Normal Equations
Let \(\bA \in \fR^{m \times n}\) with \(m \geq n\) and let \(\bb \in \fR^m\). The linear least-squares problem seeks a vector \(\hat{\bx} \in \fR^n\) satisfying \[
\begin{align}
\hat{\bx} = \arg\min_{\bx \in \fR^n} \|\bA\bx - \bb\|_2^2 = \arg\min_{\bx \in \fR^n} \sum_{i=1}^m \left( \sum_{j=1} a_{ij} x_j - b_i \right)^2.
\end{align}
\]
A vector \(\hat{\bx} \in \fR^n\) minimizes \(\|\bA\bx - \bb\|_2^2\) if and only if the residual vector \(\br = \bb - \bA\hat{\bx}\) is orthogonal to the column space \(\mathcal{R}(\bA)\): \[
\begin{align}
\br \perp \mathcal{R}(\bA) \iff \bA^T \br = \bzero \iff \bA^T(\bb - \bA\hat{\bx}) = \bzero.
\end{align}
\] This orthogonality condition produces the Normal Equations: \[
\begin{align}
\bA^T\bA \hat{\bx} = \bA^T\bb.
\end{align}
\] If \(\bA\) has full column rank (\(\operatorname{rank}(\bA) = n\)), the Gram matrix \(\bA^T\bA\) is symmetric positive definite, and the unique least-squares solution is \(\hat{\bx} = (\bA^T\bA)^{-1}\bA^T\bb\).
Define the quadratic objective function \(f(\bx) = \frac{1}{2} \|\bA\bx - \bb\|_2^2 = \frac{1}{2} (\bA\bx - \bb)^T(\bA\bx - \bb) = \frac{1}{2} \bx^T\bA^T\bA\bx - \bb^T\bA\bx + \frac{1}{2} \bb^T\bb\). Computing the gradient with respect to \(\bx\) yields \[
\begin{align}
\nabla f(\bx) = \bA^T\bA\bx - \bA^T\bb = -\bA^T(\bb - \bA\bx) = -\bA^T\br(\bx).
\end{align}
\] Setting \(\nabla f(\hat{\bx}) = \bzero\) yields \(\bA^T\bA\hat{\bx} = \bA^T\bb\). Because the Hessian \(\nabla^2 f(\bx) = \bA^T\bA\) is symmetric positive semidefinite, any stationary point is a global minimizer. When \(\bA\) has full column rank, \(\bA^T\bA\) is strictly positive definite, ensuring that \(\hat{\bx}\) is unique.
(The conditioning danger of normal equations) Although forming and solving the normal equations \(\bA^T\bA\hat{\bx} = \bA^T\bb\) via Cholesky factorization requires only \(mn^2 + \frac{1}{3}n^3\) flops, it is numerically dangerous. The 2-norm condition number of the Gram matrix is the square of the condition number of \(\bA\): \[
\begin{align}
\kappa_2(\bA^T\bA) = \kappa_2(\bA)^2.
\end{align}
\] If \(\kappa_2(\bA) = 10^8\), then \(\kappa_2(\bA^T\bA) = 10^{16}\), and solving the normal equations in standard double precision arithmetic can lose all sixteen significant digits. Therefore, the normal equations should never be the default method for ill-conditioned data fitting.
Solving Least Squares via QR Factorization
The standard, numerically stable method for solving full-rank least-squares problems is the thin QR factorization, which computes the solution without ever forming the squared matrix \(\bA^T\bA\).
Let \(\bA \in \fR^{m \times n}\) with \(m \geq n\) have full column rank, and let \(\bA = \bQ_1\bR_1\) be its thin QR factorization, where \(\bQ_1 \in \fR^{m \times n}\) has orthonormal columns and \(\bR_1 \in \fR^{n \times n}\) is nonsingular and upper triangular. The unique least-squares solution is computed by solving the upper triangular system \[
\begin{align}
\bR_1 \hat{\bx} = \bQ_1^T \bb.
\end{align}
\] The condition number of this triangular system is \(\kappa_2(\bR_1) = \kappa_2(\bA)\), completely avoiding the condition number squaring of the normal equations.
Let \(\bA = \bQ\bR = \begin{pmatrix} \bQ_1 & \bQ_2 \end{pmatrix} \begin{pmatrix} \bR_1 \\ \bzero \end{pmatrix}\) be the full QR factorization. Because orthogonal transformations preserve the Euclidean norm, we multiply the residual by \(\bQ^T\): \[
\begin{align}
\|\bA\bx - \bb\|_2^2 &= \|\bQ^T(\bA\bx - \bb)\|_2^2 \\
&= \left\| \begin{pmatrix} \bQ_1^T \\ \bQ_2^T \end{pmatrix} (\bQ_1\bR_1\bx - \bb) \right\|_2^2 \\
&= \left\| \begin{pmatrix} \bR_1\bx - \bQ_1^T\bb \\ -\bQ_2^T\bb \end{pmatrix} \right\|_2^2 \\
&= \|\bR_1\bx - \bQ_1^T\bb\|_2^2 + \|\bQ_2^T\bb\|_2^2.
\end{align}
\] The second term \(\|\bQ_2^T\bb\|_2^2\) is independent of the choice of \(\bx\), representing the inherent distance from \(\bb\) to the subspace \(\mathcal{R}(\bA)\). The first term is non-negative and can be driven to zero by setting \(\bR_1\bx = \bQ_1^T\bb\). Since \(\bA\) has full column rank, \(\bR_1\) is invertible, and \(\hat{\bx} = \bR_1^{-1}\bQ_1^T\bb\) is the unique global minimizer.
(Step-by-step least-squares line fit) Suppose we wish to fit the linear model \(y = c_0 + c_1 t\) to the three data points \((0, 1)\), \((1, 2)\), and \((2, 4)\). The overdetermined system \(\bA\bc \approx \by\) is \[
\begin{align}
\begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \end{pmatrix} \approx \begin{pmatrix} 1 \\ 2 \\ 4 \end{pmatrix}.
\end{align}
\] Computing the normal equations components: \[
\begin{align}
\bA^T\bA = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 2 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{pmatrix} = \begin{pmatrix} 3 & 3 \\ 3 & 5 \end{pmatrix}, \qquad
\bA^T\by = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 2 \end{pmatrix} \begin{pmatrix} 1 \\ 2 \\ 4 \end{pmatrix} = \begin{pmatrix} 7 \\ 10 \end{pmatrix}.
\end{align}
\] Solving \(\begin{pmatrix} 3 & 3 \\ 3 & 5 \end{pmatrix} \begin{pmatrix} c_0 \\ c_1 \end{pmatrix} = \begin{pmatrix} 7 \\ 10 \end{pmatrix}\) yields \(c_1 = 3/2 = 1.5\) and \(c_0 = (7 - 3(1.5))/3 = 5/6 \approx 0.8333\). The fitted line is \(y = \frac{5}{6} + \frac{3}{2}t\).
Rank Deficiency, SVD, and Regularization
When \(\bA\) is rank deficient (\(\operatorname{rank}(\bA) < n\)), the matrix \(\bA^T\bA\) is singular, and infinitely many vectors minimize the residual norm. In this case, the Singular Value Decomposition identifies the unique solution of minimal Euclidean norm.
Let \(\bA \in \fR^{m \times n}\) with SVD \(\bA = \bU\bsigma\bV^T = \sum_{i=1}^r \sigma_i \bu_i \bv_i^T\), where \(r = \operatorname{rank}(\bA)\). Among all vectors that minimize \(\|\bA\bx - \bb\|_2\), the unique vector with the smallest Euclidean norm \(\|\bx\|_2\) is given by the Moore-Penrose pseudoinverse solution: \[
\begin{align}
\hat{\bx}_{\text{min}} = \bA^\dagger \bb = \sum_{i=1}^r \frac{\bu_i^T\bb}{\sigma_i} \bv_i.
\end{align}
\]
When \(\bA\) is ill-conditioned or nearly rank deficient, computing \(\bA^\dagger\bb\) can amplify noise in the data due to division by very small singular values \(\sigma_i \ll 1\). Tikhonov regularization stabilizes the problem by adding an energy penalty: \[
\begin{align}
\hat{\bx}_\lambda = \arg\min_{\bx \in \fR^n} \left( \|\bA\bx - \bb\|_2^2 + \lambda \|\bx\|_2^2 \right), \qquad \lambda > 0.
\end{align}
\] The regularized solution satisfies the modified normal equations \[
\begin{align}
(\bA^T\bA + \lambda \bI) \hat{\bx}_\lambda = \bA^T\bb,
\end{align}
\] which can be solved stably via the QR factorization of the augmented matrix \(\begin{pmatrix} \bA \\ \sqrt{\lambda}\bI \end{pmatrix} \bx \approx \begin{pmatrix} \bb \\ \bzero \end{pmatrix}\).
Fit a parabola \(y = c_0 + c_1 t + c_2 t^2\) to the four data points \((-1, 1)\), \((0, 0)\), \((1, 1)\), and \((2, 5)\) using thin QR factorization.
Let \(\bA \in \fR^{m \times n}\) have full column rank with singular values \(\sigma_1 \geq ... \geq \sigma_n > 0\). Prove that \(\kappa_2(\bA^T\bA) = (\sigma_1/\sigma_n)^2 = \kappa_2(\bA)^2\).
Prove that the Tikhonov regularized solution in SVD coordinates satisfies \(\hat{\bx}_\lambda = \sum_{i=1}^n \frac{\sigma_i}{\sigma_i^2 + \lambda} (\bu_i^T\bb) \bv_i\). Explain how the filter factor \(\frac{\sigma_i^2}{\sigma_i^2 + \lambda}\) prevents noise explosion for \(\sigma_i \ll \sqrt{\lambda}\).