Solve Least-Squares Normal Equations by Elimination

The common wrong assumption is that a least-squares line fit requires some kind of curve-fitting magic that a normal equation solver cannot handle. In fact, the statistics define the linear model, and then you solve least-squares normal equations by elimination exactly as you would any 2×2 or 3×3 system. The transformation is purely algebraic: the best-fit condition A^TAx = A^Tb converts a tall, overdetermined matrix into a square, symmetric normal matrix that Gaussian elimination can handle directly. The only surprise is that the condition number of that normal matrix is the square of the condition number of the original design matrix, so what looked like a mild problem can become a numerically dangerous one.

Work a concrete case. Fit the line y = β₀ + β₁x to the points (1, 2), (2, 3), (3, 5). The design matrix A has a first column of ones and a second column of x-values: [[1, 1], [1, 2], [1, 3]]. The vector b is [2, 3, 5]^T. Compute the normal matrix A^TA: it is [[3, 6], [6, 14]]. Compute A^Tb: it is [10, 23]^T. The normal equations are [[3, 6], [6, 14]] times [β₀, β₁]^T = [10, 23]^T.

Build the augmented matrix [A^TA | A^Tb]: [[3, 6, 10], [6, 14, 23]]. The pivot in row 1 is 3. The multiplier to eliminate the entry below is 6/3 = 2. Subtract 2 times row 1 from row 2: row 2 becomes (0, 2, 3). Forward elimination is done; the system is now in row-echelon form. Back-substitute from the last row: 2β₁ = 3, so β₁ = 1.5. Then 3β₀ + 6(1.5) = 10 gives 3β₀ + 9 = 10, so β₀ = 1/3 ≈ 0.3333. The line is y = 0.3333 + 1.5x.

Check the residuals: at x = 1, y = 1.8333, residual 0.1667; at x = 2, y = 3.3333, residual -0.3333; at x = 3, y = 4.8333, residual 0.1667. Sum of squared residuals is 0.0278 + 0.1111 + 0.0278 = 0.1667. The statistics calculator that produced the normal equations from the raw data has handed you a 2×2 linear system; from that point forward the problem is pure matrix algebra.

The reason the normal equations become a linear system is that the least-squares condition sets the gradient of ||Ax - b||² to zero, which yields A^T(b - Ax) = 0, or A^TAx = A^Tb. That is a square system of n equations in n unknowns, regardless of how many data points you started with. For the line-fitting case, n = 2, so you get a 2×2 system. For a quadratic fit with three parameters you get a 3×3 system. In every case the structure is identical: symmetric, positive semi-definite, and invertible if the columns of A are linearly independent.

Gaussian elimination on these normal equations proceeds exactly as it would on any square system. Write the augmented matrix, pick the pivot, compute the multiplier, eliminate below, then back-substitute. The only extra step you need to watch for is a zero pivot. If the data points are all at the same x-coordinate, the second column of A is a scalar multiple of the first column, A^TA becomes singular, and you will hit a zero pivot during elimination. That is the failure case: the normal equations have infinitely many solutions because the model cannot distinguish slope from intercept with no x-spread. The row operation that would divide by zero simply cannot complete; a row swap will not help because the entire column is zero below the pivot.

Partial pivoting is the standard classroom strategy. For the normal equations, the normal matrix is symmetric, so partial pivoting (swapping rows only) preserves symmetry. Complete pivoting swaps columns too, which changes the variable order, and is unnecessary for the small, well-behaved systems that hand-checking typically targets. For a 3×3 system from a quadratic fit, the practical limit for hand-checking with this calculator is about a 5×5 normal matrix; beyond that the number of row operations climbs past what a person can verify reliably.

To do the calculation yourself, follow this sequence. Set up the normal equations from your data using the sums Σx, Σy, Σx², Σxy. Write the augmented matrix [A^TA | A^Tb]. Identify the pivot: the first non-zero entry in the first row. If it is zero, swap rows so the largest absolute value in that column becomes the pivot. Compute the multiplier: (entry below pivot) / (pivot). Subtract multiplier times the pivot row from each row below. Move to the next column and repeat until the matrix is in row-echelon form. Back-substitute: solve the last equation for the last variable, then substitute upward. Every row operation must be applied to the entire augmented row, including the constant column, or the solution set changes.

Failure Case: Nearly Collinear Data

For data points (0, 0), (1, 1), (2, 2), the normal equations are [[3, 3], [3, 5]] with right-hand side [3, 5]; elimination gives β₀ = 0, β₁ = 1 exactly. But change the points slightly to (0, 0.1), (1, 1), (2, 1.9). The normal matrix is still [[3, 3], [3, 5]], but the right-hand side becomes [3, 4.8]. Solving gives β₀ = 0.1, β₁ = 0.9. The condition number of the normal matrix is about 9.9, which is the square of the condition number of the original design matrix (about 3.15). For a more extreme case, say points (0, 0), (1, 0.01), (2, 2), the normal matrix is still well-conditioned, but the tiny y-value perturbation at x=1 amplifies through the elimination steps, and the resulting coefficients can shift noticeably. If you use decimal approximations instead of fractions, the round-off error from the back-substitution propagates and the calculated line can miss the true least-squares solution by several percent.

The single most practical thing to do next: before you trust any elimination result on a normal-equations problem, compute the condition number of the original design matrix A and square it. If that squared number exceeds 10⁴, switch to a QR decomposition or singular value decomposition instead of Gaussian elimination. Golub and Van Loan's Matrix Computations (Section 5.3) gives the standard warning: squaring the condition number amplifies rounding errors to the point where the elimination output may be worthless.

A second failure mode occurs when A is not full column rank. If your data has more unknowns than equations, or if two columns of the design matrix are linearly dependent, the normal matrix A^TA is singular. Elimination will produce a row of zeros on the left side. Check the augmented column: if the constant is also zero, the system has infinitely many solutions. If the constant is non-zero, the system has no solution. In either case the statistics problem is ill-posed: you cannot determine a unique best-fit line. The practical fix is to reduce the model, dropping a parameter until the columns become independent.

Newcomers forget to apply every row operation to the entire augmented row. They scale only the coefficient side of a row, or they add a multiple of one row to another but stop at the vertical bar. That changes the solution set. The augmented matrix notation [A|b] is universal precisely because it forces you to treat the constants as part of every operation. When the back-substitution gives a result that does not fit the original data, check that every row operation was applied to the full row, including the rightmost column.