The normal equations

Three readings arrive from a sensor. At time zero the reading is 1, at time one it is 2, and at time two it is 4. Plotted together the points rise steadily, but no straight line passes through all three: the first two suggest a slope of 1 and the last two suggest a slope of 2. You still want one line, so you need a rule for choosing the least bad one, and the rule should use every row rather than two of them.

Least squares chooses the line whose squared errors are smallest. For each row the error is the recorded value minus the predicted value, and the cost is the sum of those squared errors over all the rows [S:64:63]. Squaring is not decoration: it removes the sign that would otherwise let a positive error cancel a negative one, and it makes an error of 2 count four times as much as an error of 1 [S:21:235]. Minimising the length of the error list and minimising its square select the same line, so the square costs nothing and changes no answer. A residual is one of those errors, and the residual sum of squares is the total cost that the fit tries to reduce. !Three readings at times 0, 1 and 2. No line passes through all three, so least squares takes 0.8333+1.5x, whose residuals 1/6, −1/3 and 1/6 square to 1/6. The slope 1.5 is the average of the slopes the first two and the last two points suggest. Status: re-simulated.

Write the line as a prediction built from two columns. A column is one number per row: the first column is the constant 1 for every row, which lets the line start above zero, and the second column is the time. With coefficients $c_0$ and $c_1$ the prediction for row $i$ is $c_0+c_1{x}_i$. Stack the two columns side by side as a table $A$ and the recorded values as a list $\mathbf{b}$. The least-squares coefficients, collected as a list $\mathbf{c}$, solve one square system, the normal equations: square because it has as many equations as coefficients [S:21:239]. !The recorded list b lies off the plane spanned by the two columns of A, (1,1,1) and (0,1,2). Least squares replaces b by its projection p, and the residual b−p=(1/6,−1/3,1/6) is perpendicular to both columns. Two columns, two perpendicular conditions, two normal equations. Status: schematic.

The equation says that the residual is perpendicular to every column of the table, which is the projection condition of Section 8.2 written in coordinates. To see the claim, move everything to one side: $A^{\top }(\mathbf{b}-A\mathbf{c})=0$. The bracketed list is the residual, and multiplying it by the transpose takes its dot product with each column in turn, so the single matrix equation carries one perpendicular condition per column. Solving the system gives the coefficients, and the fitted values are the projection of the recorded values onto the span of the columns, the set of every prediction the table can produce [S:21:236]. The matrix on the left is built from the data alone, so the same matrix serves every list of readings taken at the same times.

Worked example. The table has rows (1, 0), (1, 1) and (1, 2), and the recorded values are 1, 2 and 4. Multiplying the table by its own transpose gives a 2 × 2 matrix whose entries are 3, 3, 3 and 5, and multiplying the transpose by the recorded values gives 7 and 10. The two equations are $3c_0+3c_1=7$ and $3c_0+5c_1=10$.

Subtracting the first equation from the second leaves $2c_1=3$, so the slope is 1.5. Substituting back gives $3c_0=2.5$, so the intercept is 0.8333. The predictions at the three times are 0.8333, 2.3333 and 3.8333, the residuals are 0.1667, −0.3333 and 0.1667, and the residual sum of squares is 1/6 = 0.1667. A slope of 1.5 is the average of the two slopes the data suggest, which is a pleasant accident of three equally spaced rows.

Check the answer against the geometry. The residual list is (1/6, −1/3, 1/6). Its dot product with the constant column is 1/6 − 1/3 + 1/6 = 0, and its dot product with the time column is 0 − 1/3 + 1/3 = 0. The error is perpendicular to both columns, which places it outside the column space.

Move from the solution by any amount inside that space: the new error is the old residual minus the movement, the two are perpendicular, and Pythagoras gives a squared length equal to the old squared length plus a positive square. The squared error therefore increases. That test is worth running on any answer the formula produces, because it checks the arithmetic without redoing it.

The formula has a cost that the algebra hides. Forming $A^{\top }A$ squares the condition number of the data table, $\mathrm{cond}(A^{\top }A)\approx (\mathrm{cond}A{)}^2$ [S:28:109], so when two columns are nearly copies the coefficients can be wrong in the third digit while the fitted values still look fine [S:28:110]. Conditioning decides how much of the recorded information survives the multiplication. The equations also need independent columns for a unique answer: columns are independent when none is a combination of the others, and a repeated column is the simplest failure, since it gives infinitely many coefficient lists that fit equally well. Section 8.5 keeps the geometry and asks what to do when no exact solution exists at all.

The formal statement. Let $A$ be a table with $n$ rows and $d$ columns, and let $\mathbf{b}$ be the list of $n$ readings. Suppose the columns of $A$ are linearly independent. Then the function $F(\mathbf{c})=\Vert A\mathbf{c}-\mathbf{b}{\Vert }^2$ has exactly one minimiser. That minimiser solves

$$A^{\mathsf{T}}A\mathbf{c}=A^{\mathsf{T}}\mathbf{b}.$$

The matrix $A^{\mathsf{T}}A$ is the Gram matrix of the columns: its entry in row $j$ and column $k$ is the dot product of column $j$ with column $k$ [S:21:239].

The derivation by calculus. Expand the square to get $F(\mathbf{c})={\mathbf{c}}^{\mathsf{T}}{A}^{\mathsf{T}}A\mathbf{c}-2{\mathbf{c}}^{\mathsf{T}}{A}^{\mathsf{T}}\mathbf{b}+{\mathbf{b}}^{\mathsf{T}}\mathbf{b}$. The gradient is $2A^{\mathsf{T}}A\mathbf{c}-2A^{\mathsf{T}}\mathbf{b}$, and a minimiser sets it to zero [S:21:238]. The second derivative is $2A^{\mathsf{T}}A$, which is positive definite when the columns are independent, so the stationary point is the unique global minimum. Setting the gradient to zero is the normal equations.

The derivation by geometry. The residual is the part of the readings that the columns cannot explain. Move $\mathbf{c}$ by any amount and the prediction moves inside the column space. The residual is shortest exactly when it is orthogonal to every direction of that space. Orthogonality to every column says $A^{\mathsf{T}}(\mathbf{b}-A\mathbf{c})=0$, which is the same equation. For any other coefficient list ${\mathbf{c}}^{\prime }$, write $\mathbf{b}-A{\mathbf{c}}^{\prime }=(\mathbf{b}-A\mathbf{c})+A(\mathbf{c}-{\mathbf{c}}^{\prime })$. The two pieces are orthogonal, so Pythagoras gives a strictly larger squared error unless ${\mathbf{c}}^{\prime }=\mathbf{c}$ [S:21:236]. Both routes describe one perpendicular foot.

The edge case. Drop independence and the theorem changes. A repeated column gives a singular Gram matrix, so the equations have infinitely many solutions and the data cannot choose between them. More columns than rows is the extreme case: when the rows are independent, $n<d$ means you can hit every reading exactly, the residual is zero, and every interpolating coefficient list ties. L1's example had three rows and two columns, so independence held and the uniqueness claim looked automatic. The same example with four rows and four coefficients is Section 8.6.

What L1 hid. L1 said the perpendicular condition "stacks row by row". The stacking runs over columns: there is one perpendicular condition per column, so one equation per coefficient. Rows enter through the residual, not through the count of equations. L1 also said $A^{\mathsf{T}}A$ "serves every list of readings taken at the same times". That is true, and it is also the warning: the matrix is built once, and its conditioning is fixed by the design, whatever the readings say [S:28:110]. The simplification was safe because the section only ever asked how the coefficients are found, not how accurately.

Where the analogy breaks. The geometric picture treats the column space as a plane you can see, and a foot of a perpendicular on a plane is unique. With dependent columns the space is still a subspace, but the coordinates of the foot are not unique. The picture stays correct and the answer stops being unique, which is the one thing a picture cannot show.

Original sources. Boyd attributes the method to Gauss and Legendre, working independently around the beginning of the nineteenth century [S:21:235]. Brunton and Kutz date Legendre's use to 1805, with Gauss developing the optimisation view later [S:52:156]. Bach traces the tool to Legendre (1805) and Gauss (1809) [S:77:61]. The normal equations are written and solved in Boyd's section 12.2 [S:21:239], and least squares opens the regression chapter of the ESL as the most popular estimation method [S:64:63].

An open question. The field has not settled which solver to use for a large, ill-conditioned least-squares problem. Forming $A^{\mathsf{T}}A$ reduces the number of rows you must store, and it squares the condition number, $\mathrm{cond}(A^{\mathsf{T}}A)\approx (\mathrm{cond}A{)}^2$ [S:28:109]; p. 110 then explains why near-coincident columns make that squared condition number large while leaving $A\mathbf{x}$ nearly unchanged [S:28:110]. A QR route avoids the squaring and costs more. An iterative route never forms the matrix and needs a good preconditioner. The right choice depends on the regime, and no rule covers every regime. Evidence that would settle a given case: backward error, forward error in the fitted values, and wall time measured on one machine across a range of condition numbers and shapes.

A frontier method. Sketch-and-solve replaces the least-squares problem with a smaller one, $\min \Vert S\mathbf{b}-SA\mathbf{c}{\Vert }^2$, where $S$ is a random projection with more rows than the table has columns and fewer than it has rows [S:77:304]. The small problem's solution approximates the large problem's solution, so the cost falls with the sketch size. It depends on the normal equations for its target and on the conditioning analysis for its error bound.

What to read next. Read Boyd's chapter 12 [S:21:235] for the derivation and the pseudoinverse, Solomon's chapter 5 [S:28:109] for why the normal equations are the wrong route in floating point, and Bach's chapter 3 [S:77:61] for the statistical view of the same estimator. Then read the sketch-and-solve section of [S:77:304] for what changes when the table is too large to factorise, and [S:52:156] for the astronomical problem that started it.