Least Squares: normal equations

The normal equations are often the first computational solution we meet for least squares. They feel almost too neat: write the model as a matrix product, do a little algebra, set a derivative to zero, and the optimal parameters appear in one compact expression.

This direct solution is possible when the model is linear in its parameters. That means the prediction can be written as:

In that case, the least squares problem becomes:

Do you remember the assumption of linearity in the parameters from the least squares foundation? This is where it becomes computationally powerful.

Imagine a simple analytical model where a measured response depends on concentration and temperature. We run five experiments and write one row per experiment:

The first column of ones is there for the intercept. It lets the model include a constant term:

The parameter vector contains the unknowns:

Putting everything together:

The vector $\mathbf{\epsilon }$ holds the measurement errors. Because of it, there is almost never a $\mathbf{\beta }$ that makes $\mathbf{X}\mathbf{\beta }$ exactly equal to $\mathbf{y}$: five equations, three unknowns. Least squares asks for the $\mathbf{\beta }$ that makes $\mathbf{X}\mathbf{\beta }$ as close as possible to $\mathbf{y}$.

If $\mathbf{X}$ were square and invertible, we could solve:

But in chemometrics, that is rarely the situation. We usually have more samples than parameters, or sometimes many more variables than samples. A calibration with 30 standards and two predictors gives a tall matrix (more rows than columns), like the one above. A spectrum with thousands of wavelengths gives a very wide one. In both cases, the ordinary inverse is not the right tool.

So instead of solving the system exactly, we look for the best approximate solution, the one that minimizes:

This is the sum of squared residuals written as a vector length.

To find the minimum, we expand the squared norm:

Multiplying it out gives four terms:

The two middle terms are the same number. Each one is a scalar (a 1 × 1 matrix), and a scalar is equal to its own transpose: ${\mathbf{y}}^T\mathbf{X}\mathbf{\beta }=({\mathbf{y}}^T\mathbf{X}\mathbf{\beta }{)}^T={\mathbf{\beta }}^T{\mathbf{X}}^T\mathbf{y}$. So they merge:

Now we use a familiar idea from calculus:

At a minimum of a smooth function, the derivative is zero.

Here, the derivative is taken with respect to the whole parameter vector $\mathbf{\beta }$, which gives the gradient: the vector of partial derivatives, one per parameter. Two rules from vector calculus do all the work. For a constant vector $\mathbf{a}$ and a symmetric matrix $\mathbf{M}$ (and ${\mathbf{X}}^T\mathbf{X}$ is symmetric):

They are the vector versions of $d(ax)/dx=a$ and $d(m{x}^2)/dx=2mx$. Applying them term by term (the first term does not depend on $\mathbf{\beta }$):

At the minimum, the gradient is the zero vector:

Divide by 2 and rearrange:

This is the key result: the normal equations [1],[2]. Whatever the number of samples, this is a small square system: one equation per parameter.

If ${\mathbf{X}}^T\mathbf{X}$ is invertible, we can isolate the parameter vector:

This is the classic closed-form solution.

It requires ${\mathbf{X}}^T\mathbf{X}$ to be invertible, which happens only when $\mathbf{X}$ has at least as many rows as columns and no column is a combination of the others (the matrix has full column rank). For a wide spectral matrix, with more wavelengths than samples, this is impossible: ${\mathbf{X}}^T\mathbf{X}$ is a $p\times p$ matrix built from only $n<p$ samples, so it is singular, and infinitely many coefficient vectors fit the data equally well. That is exactly the problem that ridge regression, principal component regression and partial least squares were invented to solve.

Suppose you are building a UV calibration model:

where $A$ is absorbance and $c$ is concentration. The design matrix is:

and the parameter vector is:

The normal equations give the intercept and slope that minimize the squared differences between measured and predicted absorbances.

The panel below shows the direct nature of the method, one step at a time: build $\mathbf{X}$, form the small matrices ${\mathbf{X}}^T\mathbf{X}$ and ${\mathbf{X}}^T\mathbf{y}$, and solve. There is no path, no sequence of guesses, and no learning rate: the algebra jumps straight to the coefficients. The first dataset is the same one used in the code at the end of this entry.

The name has nothing to do with the normal distribution. Move ${\mathbf{X}}^T\mathbf{X}\hat{\mathbf{\beta }}$ to the right-hand side and the normal equations become:

where $\mathbf{r}=\mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }}$ is the vector of residuals. Each row of this equation is the dot product of one column of $\mathbf{X}$ with the residuals, and it is zero: the residual vector is perpendicular, or normal, to every column of $\mathbf{X}$.

Geometrically, all the vectors $\mathbf{X}\mathbf{\beta }$ that the model can produce form a flat subspace (a plane, in the two-column case). Least squares picks the point of that subspace closest to the measurements $\mathbf{y}$, and the shortest way from a point to a plane is the perpendicular. The fitted values are the orthogonal projection of $\mathbf{y}$ onto the columns of $\mathbf{X}$:

$\mathbf{H}$ is called the hat matrix, because it puts the hat on $\mathbf{y}$. Its diagonal elements $h_{ii}$ measure how strongly each sample pulls the fit towards itself, its leverage. You will meet them again in multiple linear regression and in outlier detection.

When $\mathbf{X}$ has full column rank, the expression:

is the Moore-Penrose pseudoinverse. It lets us write:

That notation is useful because it reminds us that we are not simply inverting $\mathbf{X}$. We are building the best least squares solution for a system that is usually not exactly solvable.

When $\mathbf{X}$ is wide or rank deficient, the formula above fails, but the pseudoinverse still exists. It is built from the singular value decomposition (SVD), $\mathbf{X}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$:

where ${\mathbf{\Sigma }}^{+}$ inverts the nonzero singular values and leaves the zero ones at zero [3],[4]. Among the infinitely many least squares solutions, ${\mathbf{X}}^{+}\mathbf{y}$ is the one with the smallest length $\Vert \mathbf{\beta }\Vert$, the minimum-norm solution. If you also set the small singular values to zero, not only the exact zeros, you get a truncated pseudoinverse, which is principal component regression in disguise.

Using a pseudoinverse is standard practice in chemometrics. The real question is how sensitive the estimated coefficients are to noise.

A classic example is the Beer-Lambert mixture model used in classical least squares (CLS) calibration [5]:

Here $\mathbf{D}$ holds the measured spectra (samples × wavelengths), $\mathbf{C}$ the concentrations (samples × components) and $\mathbf{K}$ the pure-component spectra (components × wavelengths). If the pure spectra are known, we estimate the concentrations with the pseudoinverse of $\mathbf{K}$. But $\mathbf{K}$ is wide (few components, many wavelengths), so its pseudoinverse is written the other way round:

These are the normal equations again, written for the transposed problem: each measured spectrum is regressed on the pure spectra, and $\mathbf{K}{\mathbf{K}}^T$ plays the role of ${\mathbf{X}}^T\mathbf{X}$.

The trouble starts when two columns of the design matrix carry almost the same information, for example when two component spectra in $\mathbf{K}$ are very similar. The fitted values can still be good, but several different coefficient combinations produce almost the same fitted response, so a tiny change in the data can produce a large change in $\mathbf{\beta }$. The estimate is still a valid least squares answer; it is just very uncertain, because the data do not strongly distinguish one component from the other.

The size of this sensitivity is measured by the condition number $\kappa$, the ratio between the largest and the smallest singular value of a matrix. A large $\kappa (\mathbf{X})$ means that some direction in coefficient space is barely determined by the data. Forming ${\mathbf{X}}^T\mathbf{X}$ makes things worse, because it squares the condition number:

If $\kappa (\mathbf{X})=10^4$, then $\kappa ({\mathbf{X}}^T\mathbf{X})=10^8$, and about eight of the sixteen significant digits of double precision are lost before the solve even starts. Methods that work on $\mathbf{X}$ directly, such as QR or the SVD, only suffer $\kappa (\mathbf{X})$.

The demo below fits the same model to two designs and then repeats the fit after adding a little noise to the responses.

With independent predictors, the coefficients barely move: with this noise level the standard deviation of ${\beta }_1$ is about 0.03. With two nearly duplicated wavelengths (correlation 0.99998), it jumps to about 3, a hundred times larger, and ${\beta }_2$ moves just as much in the opposite direction. Their sum, which is what the fitted response depends on, stays as stable as before (about 0.03). That is the warning sign: good predictions, unreliable individual coefficients.

The straight line is the easiest drawing, but it is not the real boundary of the method. The real boundary is this:

If the answer is yes, the normal equations can be used. The graph may be a line, a curve, or a higher-dimensional surface. What matters is that the unknown coefficients appear as a linear combination of columns in the design matrix.

For example, this is a straight-line model:

but this curved polynomial is also linear in the parameters:

The curve bends because one column of $\mathbf{X}$ contains $x^2$. The parameter ${\beta }_2$ still multiplies that column directly. The same idea works for multiple predictors:

where one column may represent concentration and another may represent temperature. More predictors simply mean more columns in $\mathbf{X}$.

This is why "linear regression" can sometimes look nonlinear on a plot. Polynomial regression is linear regression in an expanded design matrix. Multiple linear regression is still linear regression, just with several measured variables contributing at the same time.

The method stops being directly applicable when the parameters enter the model in a genuinely nonlinear way, for example:

Here $k$ is inside the exponential, so the model cannot be solved in one normal-equation jump. That kind of problem belongs to nonlinear least squares methods such as Gauss-Newton or Levenberg-Marquardt. Notice, though, that if you fix $k$, the model is linear in $a$ and $y_0$, with columns $e^{-kx}$ and $1$. Some nonlinear solvers exploit this, a trick known as variable projection: they search only over $k$ and solve for the linear parameters with a normal-equation step at each trial value.

There is also a numerical limit. Even when the model is linear in the parameters, solving the normal equations directly can become unstable if ${\mathbf{X}}^T\mathbf{X}$ is singular or nearly singular, and expensive when there are very many parameters (the matrix has $p^2$ entries, and forming it takes about $n{p}^2$ operations). That is where the iterative methods of the next entries, gradient descent and conjugate gradient, come in.

The examples fit the calibration line of the interactive demo in two ways: with the normal equations, as derived above, and with the least squares routine you would use in practice. Both print an intercept of about 0.009 and a slope of about 2.503.

import numpy as np

c = np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
A = np.array([0.03, 0.49, 1.01, 1.48, 2.04, 2.51])

X = np.column_stack([np.ones_like(c), c])

# Normal equations: form X^T X and X^T y, then solve
beta = np.linalg.solve(X.T @ X, X.T @ A)

# What you would use in practice: least squares straight from X (SVD)
beta_lstsq, *_ = np.linalg.lstsq(X, A, rcond=None)

print(beta)         # [0.00857143 2.50285714]
print(beta_lstsq)   # same numbers

c = [0.0; 0.2; 0.4; 0.6; 0.8; 1.0];
A = [0.03; 0.49; 1.01; 1.48; 2.04; 2.51];

X = [ones(size(c)), c];

% Normal equations: form X'X and X'y, then solve
beta = (X' * X) \ (X' * A);

% What you would use in practice: backslash on X itself (QR)
beta_qr = X \ A;

disp([beta, beta_qr])   % 0.0086 and 2.5029 in both columns

conc <- c(0.0, 0.2, 0.4, 0.6, 0.8, 1.0)
A <- c(0.03, 0.49, 1.01, 1.48, 2.04, 2.51)

X <- cbind(1, conc)

# Normal equations: form X^T X and X^T y, then solve
beta <- solve(t(X) %*% X, t(X) %*% A)

# What you would use in practice: QR, directly or through lm()
beta_qr <- qr.solve(X, A)
beta_lm <- coef(lm(A ~ conc))

print(beta)      # 0.00857 and 2.50286
print(beta_qr)   # same numbers
print(beta_lm)   # same numbers

The normal equations ${\mathbf{X}}^T\mathbf{X}\mathbf{\beta }={\mathbf{X}}^T\mathbf{y}$ come from setting the gradient of the squared error to zero. They give the answer in one step whenever the model is linear in its parameters, which includes polynomials and several predictors. Geometrically, they say that the residual is perpendicular to every column of $\mathbf{X}$, and the fit is an orthogonal projection.

They need ${\mathbf{X}}^T\mathbf{X}$ to be invertible, so they fail for wide spectral data, which is why chemometrics needs ridge, PCR and PLS. And forming ${\mathbf{X}}^T\mathbf{X}$ squares the condition number, so in practice you solve with QR or the SVD, and you watch the coefficients, not only the predictions, when the predictors are nearly collinear.