The Normal Equations

The square symmetric system A^T A x = A^T b that characterizes the least-squares solution.

Derivation

To minimize the squared residual |Ax - b|^2, expand it and take the gradient with respect to x. Setting the gradient to zero yields A^T A x = A^T b, the normal equations. The name reflects geometry: the residual b - Ax must be normal (orthogonal) to the column space of A, so A^T (b - Ax) = 0, which rearranges to the same system.

Structure of the system

The matrix A^T A is square, symmetric, and positive semidefinite; it is positive definite exactly when A has linearly independent columns. In that case the normal equations have a unique solution, and because the matrix is symmetric positive definite, Cholesky factorization solves them efficiently.

The conditioning caveat

The Gram matrix A^T A has condition number equal to the square of the condition number of A. Forming it explicitly can therefore lose about twice as many digits of accuracy as working with A directly. For well-conditioned problems this is harmless, but for ill-conditioned ones the normal equations should be avoided in favor of QR or SVD-based methods.

Regularized form

Adding a penalty on the size of x, as in ridge regression, changes the system to (A^T A + alpha I) x = A^T b. The added term shifts every eigenvalue of A^T A upward by alpha, guaranteeing a positive-definite, well-conditioned matrix even when A^T A alone is nearly singular, at the price of a small bias in the solution.

import numpy as np
A = np.array([[1.0, 0.0], [1.0, 1.0], [1.0, 2.0]])
b = np.array([1.0, 2.0, 2.0])
x = np.linalg.solve(A.T @ A, A.T @ b)   # normal equations
print(x)

Parameter estimation from redundant physics measurements is usually posed as normal equations, often with regularization to stabilize the fit when observations barely constrain some parameters.