Where Loss Functions Come From

(a) Density \(p(y)=\frac{1}{\hat y}e^{-y/\hat y}\) for \(y>0\). So

Note it is scale-invariant: multiplying both \(y\) and \(\hat y\) by \(c\) changes the loss by the constant \(\log c\). Appropriate for durations and waiting times, where a 10% error matters equally at 1 second and 1 hour. (This is the Itakura–Saito divergence, used in audio spectrogram modelling for exactly that reason.)

(b) \(-\log p = \frac{(y_i-\hat y_i)^2}{2\sigma_i^2} + \log\sigma_i + \text{const}\). The \(\log\sigma_i\) does not involve \(\vw\), so

Noisier observations are down-weighted, which is exactly right: they carry less information. In matrix form the solution is \(\hat{\vw}=(\mX^\top\mathbf{\Omega}\mX)^{-1}\mX^\top\mathbf{\Omega}\vy\) with \(\mathbf{\Omega}=\operatorname{diag}(1/\sigma_i^2)\).

(c) The observed label is 1 either when the true label is 1 and no flip occurs, or when the true label is 0 and a flip does:

So the loss is ordinary cross-entropy applied to \(\tilde p\):

Why this matters: \(\tilde p\) is confined to \([\epsilon, 1-\epsilon]\), so the model can never be asked to output a probability of 0 or 1. That caps the loss for any single example, which is precisely what stops noisy labels from producing unbounded gradients. This is the principled version of label smoothing — and knowing the derivation tells you the smoothing parameter should be your estimated label-noise rate, not an arbitrary 0.1.

Continuity of value. At \(r=\delta\): the quadratic branch gives \(\tfrac12\delta^2\); the linear branch gives \(\delta(\delta-\tfrac12\delta)=\tfrac12\delta^2\). Equal ✓.

Continuity of derivative. Quadratic branch: \(\ell'(r)=r\), so \(\ell'(\delta^-)=\delta\). Linear branch (for \(r>\delta\)): \(\ell(r)=\delta r-\tfrac12\delta^2\), so \(\ell'(\delta^+)=\delta\). Equal ✓. So \(\ell_\delta\in C^1\) for every \(\delta\), which is why it is safe for gradient-based optimisation where absolute error is not.

(The second derivative is discontinuous — 1 inside, 0 outside — so \(\ell_\delta\notin C^2\). That is why Newton-type methods need care with Huber, and why log-cosh, which is smooth everywhere, is sometimes preferred.)

The implied density. Exponentiating \(-\ell_\delta\) and normalising:

This is the Huber density: Gaussian in the middle, Laplace in the tails, glued together smoothly at \(\pm\delta\). Any \(\delta\) gives a valid density (the normalising constant depends on \(\delta\)), so there is no single "correct" value — \(\delta\) parameterises a family interpolating between Gaussian (\(\delta\to\infty\)) and Laplace (\(\delta\to0\)).

Compared with a Gaussian: identical near the centre, but with far heavier tails — decaying like \(e^{-\delta|r|}\) rather than \(e^{-r^2/2}\). At \(r=6\delta\) the Huber density is larger by a factor of roughly \(e^{12.5\delta^2}\). That is the entire content of "robust": the model considers a large residual surprising but possible, rather than astronomically impossible, and so declines to distort itself to eliminate one.

Choosing \(\delta\): the standard is \(\delta=1.345\hat\sigma\), which gives 95% efficiency relative to OLS when the noise really is Gaussian while retaining robustness when it is not. Note \(\hat\sigma\) must be estimated robustly (e.g. from the median absolute deviation) or the outliers you are guarding against will inflate it.

(a) Bernoulli. \(\log p(y\mid p)=y\log p+(1-y)\log(1-p)\). \(\frac{\partial^2}{\partial p^2} = -\frac{y}{p^2}-\frac{1-y}{(1-p)^2}\). Taking \(-\E\) with \(\E[Y]=p\):

Standard error \(=\sqrt{p(1-p)/n}\) — the familiar formula, now derived rather than quoted.

(b) Poisson. \(\log p = y\log\lambda-\lambda-\log y!\), second derivative \(-y/\lambda^2\), so \(\mathcal{I}(\lambda)=\E[Y]/\lambda^2 = 1/\lambda\). Standard error \(=\sqrt{\lambda/n}\) — consistent with \(\Var=\lambda\).

(c) Gaussian mean. \(\log p = -\frac{(y-\mu)^2}{2\sigma^2}+\text{const}\), second derivative \(-1/\sigma^2\), so \(\mathcal{I}(\mu)=1/\sigma^2\), constant in \(\mu\). Standard error \(=\sigma/\sqrt n\) — Phase A Lecture 13, recovered.

Most and least precise. \(\text{se}=\sqrt{p(1-p)/n}\) is maximised at \(p=\tfrac12\) (value \(0.5/\sqrt n\)) and shrinks towards zero as \(p\to0\) or \(1\). So rare events are estimated with the smallest absolute error.

But that is the wrong way round for click-through rates, and this is the practically important point. What matters is relative precision:

This blows up as \(p\to0\). At \(p=0.001\) you need \(n\approx10^6\) for a 3% relative error, and \(n\approx10^4\) gives you \(\pm32\%\). The absolute error is tiny and completely useless.

Hence the rule of thumb that you need on the order of \(1/p\) events, not observations: A/B tests on rare conversions need enormous samples, and a reported "CTR improved from 0.10% to 0.12%" on ten thousand impressions is roughly seven conversions against six. Always compute the interval.

3. The objective. With \(\mu_i=\vw^\top\vx_i\) and \(\log\sigma_i=\vv^\top\vx_i\),

Note the structure: the first term penalises claiming large noise, the second rewards it (it shrinks the residual penalty). The model must trade them off, which is what stops it from declaring infinite noise everywhere and calling the job done.

import numpy as np
from scipy.optimize import minimize
rng = np.random.default_rng(0)
n, d = 600, 3
X = np.column_stack([np.ones(n), rng.normal(size=(n, d))])
w_true = np.array([1.0, 2.0, -1.0, 0.5])
v_true = np.array([-0.5, 0.9, 0.0, 0.0])          # noise grows with feature 1
sd = np.exp(X @ v_true)
y  = X @ w_true + sd*rng.normal(size=n)

w_ols, *_ = np.linalg.lstsq(X, y, rcond=None)
res = y - X @ w_ols
print("OLS  w:", np.round(w_ols, 3), " true:", w_true)
print("funnel check  corr(|res|, x1) =",
      round(np.corrcoef(np.abs(res), X[:,1])[0,1], 4), " (should be ~0 if homoscedastic)")

def nll(theta):
    w, v = theta[:d+1], theta[d+1:]
    ls = np.clip(X @ v, -8, 8)
    return np.sum(ls + (y - X@w)**2 / (2*np.exp(2*ls)))

opt = minimize(nll, np.concatenate([w_ols, np.zeros(d+1)]), method="L-BFGS-B")
w_h, v_h = opt.x[:d+1], opt.x[d+1:]
print("\\nHet  w:", np.round(w_h, 3))
print("Het  v:", np.round(v_h, 3), " true:", v_true)
print(f"\\n||w_ols - w_true|| = {np.linalg.norm(w_ols-w_true):.4f}")
print(f"||w_het - w_true|| = {np.linalg.norm(w_h  -w_true):.4f}")

# calibration of 90% intervals
z = 1.645
cov_ols = np.mean(np.abs(y - X@w_ols) <= z*res.std())
cov_het = np.mean(np.abs(y - X@w_h)   <= z*np.exp(X@v_h))
print(f"\\n90% interval coverage:  OLS {cov_ols:.3f}   heteroscedastic {cov_het:.3f}")

2. The OLS residual plot funnels visibly, and \(\operatorname{corr}(|r|, x_1)\) comes out around \(+0.35\) rather than 0 — the quantitative version of the funnel in Lecture 02's diagnostic figure.

4. The answer to "which improves more". The uncertainty, by a large margin.

• Point predictions improve modestly. OLS remains unbiased under heteroscedasticity — the normal equations do not require constant variance for unbiasedness, only for efficiency. So \(\vw_{\text{OLS}}\) is right on average; it is merely noisier than it needs to be, because it weights the high-noise observations as heavily as the clean ones. Expect \(\|\vw-\vw^\star\|\) to fall by perhaps 20–40%.
• Intervals improve dramatically. OLS produces one constant-width interval, which is simultaneously far too wide in the quiet region and far too narrow in the noisy one. Its 90% intervals typically cover around 93–95% overall — but that average hides the failure, since coverage is roughly 99% for low-noise points and 75% for high-noise ones. The heteroscedastic model gives per-point intervals with coverage near 90% uniformly.

The general lesson. Getting the mean right and getting the uncertainty right are separate problems, and modelling the noise fixes the second far more than the first. This is exactly why modern probabilistic regression predicts both \(\mu(\vx)\) and \(\sigma(\vx)\), why quantile regression exists, and why a model that outputs a point estimate cannot tell you when to distrust it.