Where a Loss Function Comes From: Likelihood, Not Convention

Here are the readings and the three lines, scored three ways: squared error, which everybody reaches for; absolute error, which a statistician might; and worst error, which the machinist would, because the inspector does not care about your average — he rejects the single part that is out of tolerance.

NumPy arrives here, one chapter after the pure-Python perceptron, for one reason: by the end of this chapter we evaluate four hundred thousand candidate lines against twenty readings each, and a Python loop is the wrong tool for that. It is also the notation every source cited below is written in.

import numpy as np

# Hours since the blade was changed, and the width of the part measured then.
SHIFT = np.array([
    (0.5, 20.17), (1.0, 20.28), (1.5, 20.53), (2.0, 20.61), (2.5, 20.69),
    (3.0, 20.94), (3.5, 21.21), (4.0, 21.31), (4.5, 21.27), (5.0, 21.35),
    (5.5, 21.58), (6.0, 21.80), (6.5, 21.67), (7.0, 22.07), (7.5, 22.10),
    (8.0, 22.31), (8.5, 22.48), (9.0, 22.66), (9.5, 22.90), (10.0, 23.13),
])
h, y = SHIFT[:, 0], SHIFT[:, 1]

LINES = {"A": (20.10, 0.26), "B": (20.20, 0.28), "C": (20.30, 0.26)}

for name, (a, b) in LINES.items():
    r = y - (a + b * h)                                      
    print(f"{name}   mean square {np.mean(r**2):.5f}"
          f"   mean absolute {np.mean(np.abs(r)):.5f}"
          f"   worst {np.max(np.abs(r)):.3f}")               

The quantity in the highlighted lines is the residual: what the line said minus what the caliper said, one number per reading. Every scoring rule in this chapter, and every loss function in the twenty-eight chapters after it, is some way of squashing a list of residuals down to a single number. They differ only in how they squash.

A   mean square 0.02699   mean absolute 0.12600   worst 0.430
B   mean square 0.02524   mean absolute 0.13700   worst 0.350
C   mean square 0.03179   mean absolute 0.15000   worst 0.320

Read the columns, not the rows. Squared error says B, absolute error says A, worst error says C: three rules, three winners, on the same twenty points.

I chose these three lines so that they would disagree, and I should say so plainly. The point is how easy that was — a few minutes of searching over sensible-looking intercepts and slopes turns up hundreds of such triples. The ranking is a property of the rule you picked, not a fact about the lines, so the rule is not an implementation detail: it is the definition of the problem. Which raises the question this chapter exists to answer: on what grounds do you choose it?

First a smaller matter, because there are not three lines but infinitely many. Take squared error for now, since it is what everybody takes, and shrink the problem to a single number using the trick that saved the perceptron eleven thousand epochs in Chapter 1: subtract the mean from both columns. Once the cloud of points is centred on the origin, the best line under squared error passes exactly through the origin — so the intercept is settled and only the slope is left to choose.

u, v = h - h.mean(), y - y.mean()        # 5.25 hours, 21.553 mm

def mse(theta):
    return np.mean((v - theta * u) ** 2)

grid = np.arange(0.0, 0.6001, 0.001)
curve = np.array([mse(t) for t in grid])
print(grid.size, "candidates ->", f"theta={grid[curve.argmin()]:.3f}", f"mse={curve.min():.6f}")

601 candidates -> theta=0.293 mse=0.010115

Six hundred and one candidate slopes, one winner: 0.293 millimetres per hour against a truth of 0.300. Twenty noisy readings and a for-loop got within a hundredth of a millimetre an hour — two and a third per cent.

The interesting part is not the winner but the shape of the search. Print the whole curve, rotated so the loss runs left to right:

ts = np.arange(0.0, 0.6001, 0.04)
ls = np.array([mse(t) for t in ts])
for t, l in zip(ts, ls):
    col = round(l / ls.max() * 50)
    print(f"theta={t:.2f} |{' ' * col}*{' ' * (50 - col)}| mse={l:7.4f}")

theta=0.00 |                                              *    | mse= 0.7244
theta=0.04 |                                  *                | mse= 0.5428
theta=0.08 |                        *                          | mse= 0.3878
theta=0.12 |                *                                  | mse= 0.2593
theta=0.16 |          *                                        | mse= 0.1575
theta=0.20 |     *                                             | mse= 0.0822
theta=0.24 |  *                                                | mse= 0.0336
theta=0.28 | *                                                 | mse= 0.0116
theta=0.32 | *                                                 | mse= 0.0161
theta=0.36 |   *                                               | mse= 0.0473
theta=0.40 |       *                                           | mse= 0.1050
theta=0.44 |            *                                      | mse= 0.1894
theta=0.48 |                   *                               | mse= 0.3004
theta=0.52 |                            *                      | mse= 0.4379
theta=0.56 |                                      *            | mse= 0.6021
theta=0.60 |                                                  *| mse= 0.7928

That is a valley, seen from the side. It has one bottom, the walls rise smoothly on both sides, and — this is the part the staircase of Chapter 1 could not offer — at every single point on it there is a well-defined direction of "downhill". Remember that shape. Chapter 3 is entirely about walking down it without visiting all six hundred and one points, and about what changes when a valley has more than one bottom.

We have a valley because we squared. Absolute error would have given it a kink at the bottom; worst error would have given it flat stretches where moving the line changes nothing at all. Squaring is undeniably convenient — and convenience is roughly the reason most courses give, dressed up four ways: it makes errors positive (so does absolute value); it punishes big errors more (why should it?); it is differentiable (so is the fourth power); it is what everyone uses (it is, and that is not an argument).

Here is the honest position. Squared error selected line B and absolute error selected line A. One of those is right for this factory and the other is wrong, and nothing said so far can tell you which. To choose the rule you need to know something about how the readings came to differ from the line, and that is a question about the world, not about mathematics. Answering it needs one small piece of machinery.

Here is the claim that turns "which line is better" into a question with an answer.

Assume the width of a part is the line plus a random error, and assume that error is drawn from a Gaussian — the bell curve — with mean zero and standard deviation $\sigma$:

The Gaussian's density is

Now do something the perceptron could not. For a given candidate slope $\theta$, every reading has a residual, and the formula above turns that residual into a number: how plausible is an error of exactly that size, if this slope is the truth? A reading on the line gets a big number, a reading half a millimetre off a small one.

The readings are independent — the caliper does not remember the last part — so the product rule says the plausibility of the whole notebook is the product of the individual densities. That product is the likelihood of $\theta$.¹ Note the direction, because it is the direction Bayes' rule is about: the data is fixed and known, and it is the parameter that varies. This is not "the probability of the slope". It is the probability the model assigns to the data you actually got, read as a function of the slope.

SIGMA = 0.12

def gaussian(r, sigma):
    return np.exp(-r ** 2 / (2 * sigma ** 2)) / (sigma * np.sqrt(2 * np.pi))

def likelihood(theta):
    return np.prod(gaussian(v - theta * u, SIGMA))          

for t in (0.25, 0.293, 0.35):
    print(f"theta={t}   likelihood = {likelihood(t):.6g}")

theta=0.25   likelihood = 521.952
theta=0.293   likelihood = 2.42028e+07
theta=0.35   likelihood = 0.190312

A slope of 0.293 makes this notebook forty-six thousand times more plausible than 0.25, and a hundred and twenty-seven million times more plausible than 0.35. Maximum likelihood is the principle that you pick the parameter making what you actually observed as unsurprising as possible. It is not a theorem but a proposal about what "best" ought to mean — a proposal with content, because it forces you to state your assumption about the noise before you are allowed to score anything.

Run the same three lines of code on a month of shifts instead of one, and the method falls over.

rng = np.random.default_rng(7)
u_big = rng.uniform(-5.25, 5.25, 2000)                     # 2000 readings, not 20
v_big = 0.30 * u_big + 0.12 * rng.standard_normal(2000)

print("2000 readings, sigma = 0.12 mm :", np.prod(gaussian(v_big - 0.30 * u_big, 0.12)))
noisy = 0.30 * u_big + 2.0 * rng.standard_normal(2000)
print("2000 readings, sigma = 2.00 mm :", np.prod(gaussian(noisy - 0.30 * u_big, 2.0)))
print("largest float64 :", np.finfo(np.float64).max)

RuntimeWarning: overflow encountered in reduce
2000 readings, sigma = 0.12 mm : inf
2000 readings, sigma = 2.00 mm : 0.0
largest float64 : 1.7976931348623157e+308

Two thousand multiplications and the answer is inf. Change one constant — a sloppier caliper, so the densities come out smaller than 1 instead of larger — and the same code returns 0.0. Both answers are wrong, in opposite directions, neither raises an exception you can catch, and the second does not even print a warning.

Nothing is wrong with the mathematics. The likelihood at those settings is a perfectly well-defined finite number: its natural logarithm is 1400.91, so the number itself is about $10^{608}$. The problem is that your computer does not have that number, and it is worth understanding exactly which numbers it does have, because this is not the last time it will decide the outcome.

The fix for the exploding product is the usual one: take logarithms. The logarithm turns products into sums, it is strictly increasing so it cannot move the location of the maximum, and a sum of two thousand moderate numbers is something float64 handles without complaint. By convention we take the negative log-likelihood, so that better means smaller. Now substitute the Gaussian density and watch what happens.

That is the result the chapter exists for, and it deserves to be stated without hedging: squared error is not a convention. It is the negative log-likelihood of a Gaussian, with the constants removed. Minimising squared error is precisely the same act as asserting that your errors are Gaussian and asking which parameter makes your data least surprising. You were making that assertion all along; you were just not being told.

The equivalence is checkable, so check it: scan the same six hundred and one slopes with the full negative log-likelihood, constants and all, and with plain squared error.

N = v.size

def nll(theta):
    r = v - theta * u
    return N * np.log(SIGMA * np.sqrt(2 * np.pi)) + np.sum(r ** 2) / (2 * SIGMA ** 2)

nlls = np.array([nll(t) for t in grid])
mses = np.array([mse(t) for t in grid])
print(f"argmin of the negative log-likelihood : theta={grid[nlls.argmin()]:.3f}  nll={nlls.min():.6f}")
print(f"argmin of the mean squared error      : theta={grid[mses.argmin()]:.3f}  mse={mses.min():.6f}")
print("same index:", nlls.argmin() == mses.argmin())

argmin of the negative log-likelihood : theta=0.293  nll=-17.001977
argmin of the mean squared error      : theta=0.293  mse=0.010115
same index: True

Different numbers on the vertical axis, and one of them is negative, which a sum of squares never is: a negative log-likelihood may go below zero, because a density can exceed 1. Same bottom of the same valley, to the last grid point.

We fixed $\sigma$ at 0.12 by fiat, and nobody at the plant knows the spread of their caliper's error. Treat it as a second unknown and let maximum likelihood decide it too. Here the constant term we just discarded comes back, because it is the only thing standing between the model and a claim of perfect precision.

r = v - 0.293 * u
sigmas = np.arange(0.01, 1.0001, 0.0001)
nll_sigma = N * np.log(sigmas * np.sqrt(2 * np.pi)) + np.sum(r ** 2) / (2 * sigmas ** 2)

print("best sigma on the grid       :", round(float(sigmas[nll_sigma.argmin()]), 4))
print("sqrt(mean squared residual)  :", round(float(np.sqrt(np.mean(r ** 2))), 4))

best sigma on the grid       : 0.1006
sqrt(mean squared residual)  : 0.1006

The two agree to four decimals, and not by accident: differentiating that expression and setting it to zero gives $\hat{\sigma}^2 = \frac{1}{N}\sum r_i^2$ exactly. So the mean squared error is not merely like a variance. Under this model it is the maximum-likelihood estimate of the variance of the noise — the number you have been minimising all along was an estimate of how noisy your sensor is.

One wrinkle, cheap to state and expensive to rediscover later: that estimate is biased low, because the residuals were measured against a fit that was itself chosen to make them small. Simulate it — two hundred thousand notebooks of twenty readings each, drawn from a distribution whose true variance is exactly 1, with the one parameter of the fit estimated from the readings themselves. Dividing the sum of squares by $N$ gives an average of 0.9501; dividing by $N-1$ gives 1.0001; and $(N-1)/N$ is 0.95 on the nose. Every parameter you fit costs one degree of freedom, and this is the smallest visible instance of a much larger problem: a model always looks better on the data it was fitted to. Chapter 4 turns that into the discipline of holding data back, and Chapter 6 gives the effect its name.

If squared error asserts that the noise is Gaussian, the next question is what happens when the assertion is false. Not slightly false — false in the way real measurements are false.

On the shop floor, most caliper readings are good to a tenth of a millimetre, and once or twice a shift a chip of swarf gets under the jaw and the reading is off by several millimetres. Errors like that are heavy-tailed: small most of the time, occasionally enormous, and enormous far more often than a bell curve allows. The Cauchy distribution is the standard clean model of that behaviour, and its density is as simple as the Gaussian's:

The difference is the tail: the Gaussian falls off like $e^{-\varepsilon^2}$, brutally fast, and the Cauchy like $1/\varepsilon^2$, barely at all. The consequence is easier to see than to say:

rng = np.random.default_rng(3)
g = 0.12 * rng.standard_normal(10 ** 6)          # Gaussian noise
c = 0.12 * rng.standard_cauchy(10 ** 6)          # Cauchy noise, same scale
for k in (10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
    print(f"{k:>9,} samples   gaussian var {g[:k].var():.4f}   cauchy var {c[:k].var():10.2f}")

      100 samples   gaussian var 0.0164   cauchy var       0.26
    1,000 samples   gaussian var 0.0146   cauchy var      59.88
   10,000 samples   gaussian var 0.0145   cauchy var     358.17
  100,000 samples   gaussian var 0.0144   cauchy var    3097.98
1,000,000 samples   gaussian var 0.0144   cauchy var   32886.10

The Gaussian's sample variance settles on 0.0144, which is $0.12^2$, and stays there. The Cauchy's climbs, and keeps climbing for as long as you sample, because there is nothing for it to converge to: the Cauchy distribution has no variance, and no mean either. Squared error, whose whole business is minimising an average of squares, is being asked for a quantity that does not exist.

So here is one shift where the caliper was fooled. Same twenty hours, same blade, same drift of 0.30 millimetres per hour — only the noise is now Cauchy. Fit it twice: once by minimising squared residuals, once by minimising the negative log-likelihood of the noise that actually generated the data. The centring trick is no help here — it pins the intercept only for squared error — so both fits go by brute force over a grid of intercepts and slopes, since we still have no way to find the bottom of a valley except by visiting it.

SWARF = np.array([
    (0.5, 20.08), (1.0, 21.95), (1.5, 20.86), (2.0, 27.51), (2.5, 20.64),
    (3.0, 20.75), (3.5, 21.01), (4.0, 21.03), (4.5, 21.37), (5.0, 20.60),
    (5.5, 22.03), (6.0, 21.95), (6.5, 21.98), (7.0, 22.01), (7.5, 21.73),
    (8.0, 22.97), (8.5, 22.60), (9.0, 22.66), (9.5, 22.44), (10.0, 22.78),
])
hs, ys = SWARF[:, 0], SWARF[:, 1]

A = np.arange(18.0, 22.001, 0.005)      # 801 intercepts
B = np.arange(-0.20, 0.8001, 0.002)     # 501 slopes
R = ys - (A[:, None, None] + B[None, :, None] * hs)      # every line against every point

SCALE = 0.12
square = np.sum(R ** 2, axis=2)                          # least squares          
cauchy = np.sum(np.log(1 + (R / SCALE) ** 2), axis=2)    # Cauchy likelihood      

for name, surface in (("least squares", square), ("Cauchy likelihood", cauchy)):
    i, j = np.unravel_index(surface.argmin(), surface.shape)
    print(f"{name:>18}:  width = {A[i]:.3f} + {B[j]:.4f} * hours"
          f"   -> 23.5 mm at hour {(23.5 - A[i]) / B[j]:.2f}")
print(f"{'the truth':>18}:  width = 20.000 + 0.3000 * hours"
      f"   -> 23.5 mm at hour {(23.5 - 20.0) / 0.30:.2f}")
print(f"{A.size * B.size:,} candidate lines evaluated")

The two highlighted lines are the entire difference between the fits. Take the log of the Cauchy density, drop the constants exactly as before, and $\sum \log\left(1 + (r/s)^2\right)$ is what survives. Same recipe, different claim about the noise.

     least squares:  width = 21.380 + 0.1080 * hours   -> 23.5 mm at hour 19.63
 Cauchy likelihood:  width = 19.935 + 0.3020 * hours   -> 23.5 mm at hour 11.80
         the truth:  width = 20.000 + 0.3000 * hours   -> 23.5 mm at hour 11.67
401,301 candidate lines evaluated

Least squares reports a drift of 0.108 millimetres an hour, roughly a third of the real rate, and concludes that the blade is good until hour 19.6. The true answer is hour 11.7. Acting on that fit, the plant runs the press for eight extra hours making parts that are out of tolerance, on the authority of the most standard loss function in the field. The Cauchy fit, using the same twenty readings, the same grid, and one line's difference in the code, lands on hour 11.8.

Two objections deserve answers, because both are the first thing a good engineer says.

The outlier is obvious — just delete it. You can, and it helps, and it is not enough. Deleting the single worst reading moves the least-squares slope from 0.108 to 0.239, which still puts the blade change at hour 13.1, an hour and a half late; deleting the worst, refitting, and deleting whatever is worst now gets you to 0.286 — and note that this is already a procedure, not an observation: delete the two largest residuals of the original fit instead and you land on 0.223. But you have now made judgement calls you cannot write down or defend, and automating the rule does not rescue it: drop-the-largest-residual-then-refit, run over a thousand simulated shifts, has a median slope error of 0.0177 against the likelihood fit's 0.0100, and is off by more than 0.05 in 14.7% of shifts against 1.3%. Deletion is a patch on top of a wrong assumption. The likelihood needs no patch, because it never assumed the outlier was impossible.

You picked a lucky dataset. That objection is exactly right, which is why the last experiment simulates a thousand independent shifts and refits both ways on each.

A = np.arange(18.0, 22.001, 0.02)        # a coarser grid: a thousand fits to do
B = np.arange(-0.20, 0.8001, 0.005)
lines = A[:, None, None] + B[None, :, None] * hs
rng = np.random.default_rng(2026)
err_sq, err_ca = [], []

for _ in range(1000):                                        # 1000 independent shifts
    ys = 20.00 + 0.30 * hs + SCALE * rng.standard_cauchy(hs.size)
    R = ys - lines
    _, j = np.unravel_index(np.sum(R ** 2, axis=2).argmin(), (A.size, B.size))
    _, q = np.unravel_index(np.sum(np.log1p((R / SCALE) ** 2), axis=2).argmin(), (A.size, B.size))
    err_sq.append(abs(B[j] - 0.30))
    err_ca.append(abs(B[q] - 0.30))

err_sq, err_ca = np.array(err_sq), np.array(err_ca)
for name, e in (("least squares", err_sq), ("Cauchy likelihood", err_ca)):
    print(f"{name:>18}: median slope error {np.median(e):.4f} mm/h"
          f"   off by more than 0.05 in {100 * np.mean(e > 0.05):4.1f}% of shifts"
          f"   worst {e.max():.3f}")
print(f"the likelihood fit is the closer of the two in {100 * np.mean(err_ca < err_sq):.1f}% of shifts")

     least squares: median slope error 0.0350 mm/h   off by more than 0.05 in 40.4% of shifts   worst 0.500
 Cauchy likelihood: median slope error 0.0100 mm/h   off by more than 0.05 in  1.3% of shifts   worst 0.090
the likelihood fit is the closer of the two in 75.6% of shifts

Median, not mean, for the same reason as everything else in this section: the least-squares errors are driven by a Cauchy, so their average is not a stable thing to report. Least squares is badly wrong two shifts in five; the likelihood fit is badly wrong one shift in seventy-seven, and its worst failure across a thousand shifts is less than a fifth of least squares' worst.

None of this makes squared error bad. It makes it specific, and the arithmetic says exactly why. Take a residual of 0.1 mm and one of 7 mm. Squared, the bad reading contributes 4,900 times as much to the total as the good one, so the line is dragged bodily toward it; under the Cauchy log-likelihood the same two residuals contribute 0.527 and 8.133, a ratio of 15.4. The bad reading still counts, it just does not get to decide. This is the beginning of robust statistics, where Huber's 1964 loss splits the difference by behaving quadratically for small residuals and linearly for large ones,⁷ and where Tukey had already shown how little contamination it takes to make the sample variance a worse tool than the mean absolute deviation.⁸

One historical note, too good to leave out. Least squares was published first, by Legendre in 1805, as a convenient algebraic device with no justification beyond that it worked.⁹ Four years later Gauss ran the argument backwards: he took it as given that the arithmetic mean is the right way to combine repeated measurements, asked which error distribution makes the mean the most probable value, and showed that essentially only one does — the one now named after him.¹⁰ The derivation in this chapter is his, it is more than two centuries old, and it is still the part most courses leave out.

Look back at the valley, though. Standing at $\theta = 0.20$ with a loss of 0.0822, the direction of "downhill" is no mystery — you can see it on the page, the curve slopes down to the right. If you could ask the loss function which way it slopes at the point you are standing on, without evaluating it anywhere else, you could take a step that way, ask again, and repeat until the ground is flat.

That question has a name. The slope of a function at a point is its derivative, and for a function of many parameters the collection of slopes in every direction at once is the gradient. Chapter 1 could not use one, because the perceptron's error was a staircase with no slope to ask about. This chapter has built something better: a loss that is smooth everywhere and that came from a stated assumption rather than a preference.

So the question for Chapter 3 is no longer whether a slope exists. It is how to compute it, why moving against it goes downhill rather than uphill — a sign almost every course asks you to take on faith — and how far to step before asking again, which turns out to be the one number that decides whether a training run converges, oscillates around the answer forever, or runs off to infinity.

Also worth reading alongside this chapter: Prince, Understanding Deep Learning §5.1–5.2 and Appendix C, which builds every loss in the book from maximum likelihood in the order used here; Goodfellow, Bengio and Courville, Deep Learning §3.1–3.11 and §5.5, whose maximum-likelihood section also derives the KL divergence that Chapter 4 needs; Murphy, Probabilistic Machine Learning: An Introduction chapter 2 and §4.2, on what maximum likelihood does and does not guarantee; Deisenroth, Faisal and Ong, Mathematics for Machine Learning §6.1–6.4 for the sum rule, product rule and Bayes' rule done properly; Tom Mitchell's short CMU note Estimating Probabilities: MLE and MAP (2016); and §22.7 of Dive into Deep Learning, which reaches the same result in runnable code.