Maximum Likelihood = Least Squares

Look almost anywhere in science and engineering — fitting a line through data, calibrating a sensor, training a model, reconstructing an image — and you will find the same recipe: pick the answer that makes the sum of squared errors as small as possible. Least squares is everywhere. So common that it feels like a law of nature.

But why squared? Why not the sum of the plain errors, or the sum of their absolute values, or the largest single error? Is least squares just a convenient choice — squares are smooth and easy to differentiate — or is there a deeper reason it keeps winning? There is, and it is one of the most beautiful coincidences in applied mathematics:

When measurement noise is Gaussian, the least-squares fit is exactly the maximum-likelihood estimate — the single most probable explanation of the data.
Minimising \lVert Gm - d\rVert^2 and maximising the probability of the data are the same optimisation.

The humble line of best fit is secretly a statement of maximum probability. Let us see why.

From a product of bell curves to a sum of squares

Suppose each measurement d_i is the true model prediction (Gm)_i plus some independent Gaussian noise of standard deviation \sigma. The probability of seeing that one measurement is a bell curve centred on the prediction:

p(d_i \mid m) = \frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\Big(-\frac{(d_i - (Gm)_i)^2}{2\sigma^2}\Big).

Because the noise on each measurement is independent, the probability of the whole dataset is the product of these bumps:

p(d \mid m) = \prod_{i=1}^{N} \frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\Big(-\frac{(d_i - (Gm)_i)^2}{2\sigma^2}\Big).

Products are awkward to maximise, so take the logarithm — a monotone function, so whatever maximises the likelihood also maximises its log. The log of a product is a sum, and the log undoes the exponential, leaving just the exponents:

\log p(d \mid m) = \text{const} - \frac{1}{2\sigma^2}\sum_{i=1}^{N} \big(d_i - (Gm)_i\big)^2.

The constant and the \sigma do not depend on m. So maximising the likelihood is identical to minimising that sum of squared residuals — the least-squares objective \lVert d - Gm\rVert^2. The squares were never a convenience; they fell out of the Gaussian's squared exponent.

The general case: the covariance does the weighting

Written with a full noise covariance e \sim N(0, C_D), the likelihood is

p(d \mid m) \propto \exp\!\Big(-\tfrac12 (d - Gm)^{\mathsf T} C_D^{-1} (d - Gm)\Big),

so maximum likelihood is weighted least squares:

\hat m_{\text{MLE}} = \arg\min_m\, (d - Gm)^{\mathsf T} C_D^{-1} (d - Gm).

The inverse covariance C_D^{-1} is the natural weighting matrix: a measurement with small variance (trusted) gets a large weight; a noisy one gets little say. When the noise is uniform and uncorrelated, C_D = \sigma^2 I and this collapses to ordinary least squares — the special case we started from. The general estimate is

\hat m_{\text{MLE}} = (G^{\mathsf T} C_D^{-1} G)^{-1} G^{\mathsf T} C_D^{-1} d.

This is the bridge between the deterministic and statistical stories: the least-squares method we used for stability was, all along, the maximum-likelihood estimate under Gaussian noise. It still offers no help with ill-posedness, though — for that we must add a prior, which is the next step.

Gaussian noise ⇒ likelihood \propto \exp(-\tfrac12 (d-Gm)^{\mathsf T}C_D^{-1}(d-Gm)).
MLE = weighted least squares with weight C_D^{-1}; uniform noise C_D = \sigma^2 I gives ordinary least squares.
Estimate: \hat m = (G^{\mathsf T}C_D^{-1}G)^{-1}G^{\mathsf T}C_D^{-1}d.

Worked example: two noisy readings of one number

Take the simplest inverse problem imaginable. You want a single number m — say a voltage — and you measure it twice with the same instrument, getting d_1 = 4 and d_2 = 10. Each reading is m plus independent Gaussian noise. What single value is most probable?

The log-likelihood is -\tfrac{1}{2\sigma^2}\big[(4-m)^2 + (10-m)^2\big]. Differentiate and set to zero: (4-m) + (10-m) = 0 \Rightarrow m = 7 — the plain average. The two Gaussian bumps, one peaked at 4 and one at 10, multiply into a single narrower bump peaked halfway between. Slide the guess below and watch the squared-error bowl (bottom = least squares) line up perfectly with the peak of the likelihood curve.

The lesson generalises: for Gaussian noise the maximum-likelihood point sits exactly where the sum of squared errors bottoms out. Two views, one answer.

Worked example: why squared, not absolute?

Where does the squaring physically come from? Straight from the Gaussian's exponent, which is -(d-m)^2/2\sigma^2 — a squared distance. Compare three candidate error measures for the readings 4 and 10, and ask what each one calls "best":

Squared error (4-m)^2 + (10-m)^2 is minimised at the mean, m = 7. This is Gaussian MLE.
Absolute error |4-m| + |10-m| is minimised by the median (anything between 4 and 10, here also 7). This is the MLE for Laplace noise.
Maximum error \max(|4-m|,|10-m|) is minimised at the midrange — again 7 here, but it obsesses over the single worst point.

With only two symmetric points they happen to agree, but the machinery is different: the noise distribution you assume chooses the error measure you should minimise. Assume a bell curve and the arithmetic hands you squares.

Worked example: fitting a straight line

Fitting y = a + bx through noisy points (x_i, y_i) is the same story with m = (a, b). The Gaussian log-likelihood becomes

-\frac{1}{2\sigma^2}\sum_i \big(y_i - a - b x_i\big)^2,

so the most probable line is the one minimising the summed squared vertical gaps — precisely the classic least-squares regression line \hat m = (G^{\mathsf T}G)^{-1}G^{\mathsf T}d with G = [\,\mathbf 1 \;\; x\,]. Every "line of best fit" a spreadsheet has ever drawn is quietly a maximum-likelihood estimate under the assumption that the scatter about the line is Gaussian.

The equivalence "least squares = maximum likelihood" holds only under the Gaussian noise assumption. Change the noise and you change the best-fit recipe:

If the noise is heavy-tailed or has outliers — not Gaussian — least squares is no longer the maximum-likelihood choice, and it becomes badly non-robust. A single wild outlier can wreck the entire fit, because squaring an error of 100 contributes 10,000 to the objective, so the fit bends over backwards to appease that one rogue point.
For such data, use a robust loss. Minimising absolute error (\ell_1) is the MLE for Laplace-distributed noise and shrugs off outliers, because an error of 100 costs only 100, not 10,000. Huber loss and others interpolate between the two.

So before reaching for least squares, ask: is my noise really Gaussian? If it is riddled with outliers, the "optimal" least-squares fit is optimal for the wrong problem.

Essentially, yes — and the story is delicious. In the early 1800s Carl Friedrich Gauss was convinced that the arithmetic mean ought to be the "best" estimate of a repeatedly measured quantity, and that least squares was the right general method. So he ran the logic backwards: he asked which noise distribution would make least squares optimal — for which distribution is the most probable estimate exactly the least-squares one? The answer he derived was the bell curve we now call the Gaussian. He wanted least squares to be optimal, so he defined the distribution that makes it so — and then found, wonderfully, that this same curve really does describe the errors of real instruments (thanks to the central limit theorem piling up many tiny independent nudges).

That is why the bell curve and least-squares fitting are twins throughout science: they are two faces of one assumption. Swap the assumed noise distribution and the "best fit" recipe changes with it — Gaussian gives squared error (\ell_2), Laplace gives absolute error (\ell_1) — the very same distinction that later reappears as L_2 versus L_1 regularization.