The Normal Equations and Closed-Form Solution

Suppose we have $n$ training examples, each with $d$ features, stacked into a design matrix $X \in \mathbb{R}^{n \times d}$, and a target vector $y \in \mathbb{R}^n$. We assume a linear model $\hat{y} = Xw$ for weight vector $w \in \mathbb{R}^d$, and we measure error using the residual sum of squares:

The least squares problem is to find $w^\star = \arg\min_w J(w)$. This is a convex quadratic function of $w$, so any point where its gradient vanishes is a global minimum (a fact we justify below). The key geometric insight is this: as $w$ ranges over all of $\mathbb{R}^d$, the fitted vector $Xw$ ranges over the column space of $X$, denoted $\mathrm{col}(X)$ — the set of all linear combinations of the columns of $X$. Minimizing $\|y - Xw\|$ is therefore equivalent to finding the point in $\mathrm{col}(X)$ closest to $y$ in Euclidean distance. Basic Euclidean geometry tells us that the closest point in a subspace to an external point is found by orthogonal projection: drop a perpendicular from $y$ onto the subspace. The foot of that perpendicular is $\hat{y} = Xw^\star$, and the connecting vector $y - \hat{y}$ is orthogonal to the entire subspace, not just to $\hat{y}$. "Orthogonal to the entire subspace" means orthogonal to every column of $X$. Writing $r = y - Xw^\star$ for the residual, this condition is

This system of $d$ linear equations in $d$ unknowns is called the normal equations — "normal" in the old sense of "perpendicular." When $X^\top X$ is invertible (equivalently, when $X$ has full column rank, meaning no feature is an exact linear combination of the others), we can solve directly:

The matrix $X^+ = (X^\top X)^{-1}X^\top$ is called the Moore–Penrose pseudoinverse of $X$ in the full-rank case, and $w^\star = X^+ y$ generalizes the notion of "inverting" a non-square matrix.

The algorithm is short precisely because the hard work is a linear algebra fact, not a search procedure. Given the data, we build two small matrices and solve one linear system. ``` ALGORITHM: NormalEquations(X, y) Input: X, an n×d matrix of features (one row per example), y, an n×1 vector of targets Output: w, the d×1 least-squares weight vector 1. # Optionally prepend a column of ones to X for an intercept term 2. A ← Xᵀ X # d×d matrix, the "Gram matrix" 3. b ← Xᵀ y # d×1 vector 4. if A is singular (or ill-conditioned): 5. report failure or switch to a regularized / pseudoinverse solve 6. else: 7. solve the linear system A w = b for w 8. # in practice: Cholesky factorization of A, then two triangular solves 9. return w ``` Two implementation details matter in practice. First, we never explicitly compute $(X^\top X)^{-1}$ as a matrix and then multiply it by $X^\top y$. Explicit matrix inversion is numerically wasteful and less stable. Instead we factor $A = X^\top X$ — since $A$ is symmetric positive semi-definite, a Cholesky decomposition $A = LL^\top$ exists whenever $A$ is positive definite (full column rank), and solving $Aw = b$ reduces to two cheap triangular solves: $Lz = b$ then $L^\top w = z$. Second, forming $X^\top X$ explicitly squares the condition number of the problem, which can amplify numerical error when features are nearly collinear. A more numerically robust alternative factors $X$ itself using QR decomposition, $X = QR$ with $Q$ orthogonal and $R$ upper triangular, and solves the equivalent triangular system $Rw = Q^\top y$ directly, avoiding the squaring altogether. Both routes solve the same normal equations; QR is preferred in production numerical libraries, while the $X^\top X$ route is the one most often taught and hand-computed because $d$ is small in the pedagogical case.

Correctness. We show two things: (1) the normal equations characterize every stationary point of $J(w)$, and (2) when $X$ has full column rank, that stationary point is the unique global minimizer. For (1), expand and differentiate:

Since $J$ is a quadratic form in $w$, its gradient is

Setting the gradient to zero gives exactly $X^\top X w = X^\top y$, the normal equations. This confirms the calculus route agrees with the geometric orthogonality route derived above — both say the same thing, that $X^\top(y - Xw) = 0$. For (2), we check the Hessian: $\nabla^2_w J(w) = 2X^\top X$. A matrix of the form $X^\top X$ is always positive semi-definite, because for any vector $v$,

If $X$ has full column rank (its columns are linearly independent), then $Xv = 0$ implies $v = 0$, so the inequality above is strict for $v \neq 0$: $X^\top X$ is positive definite. A positive definite Hessian means $J$ is strictly convex, and a strictly convex function has at most one stationary point, which must be its unique global minimum. Combined with (1), this proves $w^\star = (X^\top X)^{-1}X^\top y$ is the unique minimizer of the least-squares objective. The projection argument gives the same conclusion more directly and more memorably: among all vectors in $\mathrm{col}(X)$, the one at minimum Euclidean distance from $y$ is unique (Euclidean projection onto a subspace is always unique), and it is characterized exactly by the orthogonality of the residual — which is the normal equations.

When it fails. If $X$ does not have full column rank — for example, if two features are duplicates or one is a linear combination of others — then $X^\top X$ is only positive semi-definite, not invertible, and infinitely many $w$ achieve the same minimal loss (they differ by vectors in the null space of $X$). In this regime, one instead uses the Moore–Penrose pseudoinverse (via the singular value decomposition) to select the minimum-norm solution, or adds an L2 penalty as in Ridge Regression and L2 Regularization, which makes $X^\top X + \lambda I$ invertible for any $\lambda > 0$ regardless of rank.

Time complexity. Forming $X^\top X$ costs $O(nd^2)$ multiply-adds (a $d \times n$ matrix times an $n \times d$ matrix). Forming $X^\top y$ costs $O(nd)$. Solving the resulting $d \times d$ system via Cholesky factorization costs $O(d^3)$. Total time is

When $n \gg d$ (many samples, modest feature count), the $nd^2$ term dominates and the cost is roughly linear in $n$ — very efficient. When $d$ is large (thousands of features), the $d^3$ factorization term dominates and becomes the bottleneck; this is precisely the regime where iterative methods like Gradient Descent for Parameter Estimation become attractive, since a single gradient step costs only $O(nd)$.

Space complexity. Storing $X$ costs $O(nd)$; storing $X^\top X$ costs $O(d^2)$, independent of $n$. This is a genuine advantage of the normal-equation approach for streaming or very large $n$: once $X^\top X$ and $X^\top y$ are accumulated (which can even be done incrementally, row by row, since both are sums over examples), the raw data need not be kept in memory.

Tradeoffs, honestly. The closed-form solution is exact and requires no learning rate, no convergence check, and no iteration count — real advantages for small-to-moderate $d$. But it does not scale gracefully to very high-dimensional feature spaces ($d^3$ becomes intolerable past a few thousand features), it requires $X^\top X$ to be well-conditioned (near-collinear features cause numerical blowup even when technically invertible), and it has no natural mechanism to prevent overfitting when $d$ approaches or exceeds $n$ — all reasons the field also relies on regularized and iterative alternatives.

Suppose we want to fit a line $\hat{y} = w_0 + w_1 x$ to four points relating hours studied ($x$) to exam score ($y$):

We prepend a column of ones to capture the intercept $w_0$:

Step 1 — form $X^\top X$.

The $(1,1)$ entry is $n=4$ (count of rows); the off-diagonal entries are $\sum x_i = 10$; the $(2,2)$ entry is $\sum x_i^2 = 1+4+9+16=30$.

Step 2 — form $X^\top y$.

Here $\sum y_i = 3+4+6+9=22$, and $\sum x_i y_i = 1(3)+2(4)+3(6)+4(9)=3+8+18+36=65$.

Step 3 — solve $\begin{bmatrix} 4 & 10 \\ 10 & 30\end{bmatrix} w = \begin{bmatrix}22\\65\end{bmatrix}$. The determinant is $4(30) - 10(10) = 120 - 100 = 20 \neq 0$, so $X^\top X$ is invertible, confirming a unique solution exists. Its inverse is

Multiply:

So $w_0 = 0.5$, $w_1 = 2.0$: the fitted line is $\hat{y} = 0.5 + 2x$.

Verification via the orthogonality condition. Compute residuals $r_i = y_i - \hat y_i$:

Check $X^\top r = 0$: the sum of residuals is $0.5 - 0.5 - 0.5 + 0.5 = 0$ (orthogonal to the ones-column, confirming residuals average to zero), and $\sum x_i r_i = 1(0.5)+2(-0.5)+3(-0.5)+4(0.5) = 0.5-1-1.5+2 = 0$ (orthogonal to the $x$-column). Both dot products vanish exactly as the theory demands — the residual really is perpendicular to the column space.

The method of least squares predates modern machine learning by two centuries. Its discovery is credited to Legendre, who published the method first, and to Gauss, who claimed prior use and gave the deeper probabilistic justification connecting it to the normal (Gaussian) distribution of errors — the historical root of the word "normal" in "normal equations."

The normal equations are the backbone of ordinary least squares regression, still the default first model fit in statistics, econometrics, and applied science whenever a linear relationship is plausible and the feature count is modest. Scientific computing libraries — the linear algebra routines behind `numpy.linalg.lstsq`, R's `lm()`, and the normal-equation solvers inside statistical packages — implement exactly this method, usually via QR or SVD factorizations for numerical stability rather than the raw $(X^\top X)^{-1}$ formula. Beyond direct regression, the same projection idea recurs throughout machine learning: it is the basis for understanding why Linear Regression and the Least Squares Objective has a unique optimum, it motivates the regularized variants in Ridge Regression and L2 Regularization used whenever features are correlated or numerous, and it provides the exact baseline against which iterative solvers like Gradient Descent for Parameter Estimation are checked for correctness on small problems. Calibration steps in sensor fusion, camera geometry (fitting lines and planes to noisy point measurements), signal processing (least-squares filter design), and the final linear layer of many machine learning pipelines all reduce, at some point, to solving a system that looks exactly like the normal equations.

Ask the tutor: How does adding an L2 penalty term change the normal equations, and why does it always restore invertibility?