Overfitting, explained by fitting a sine wave ten different ways

The first assignment of the intelligent-systems module handed us ten numbers and asked us to find the function that made them. Ten numbers is not much. The whole point of the exercise, which I did not appreciate until I plotted the ninth-order fit, is that ten numbers is exactly enough to be lied to by your own error function.

Ten points and a polynomial

The setup is one input $x$, one target $t$, and $N$ pairs of them. I want the function $f(x)$ that generated the targets, knowing that something noisy happened between $f(x)$ and the $t$ I was given.

I fit a polynomial, and the report’s justification for that choice is the one I still give: most functions you care about are well approximated by the first few terms of their Taylor expansion, so a polynomial is a reasonable general-purpose guess when you have no reason to prefer anything else.

$$y(x, \mathbf{w}) \;=\; \sum_{j=0}^{M} w_j x^j \;=\; \mathbf{w}^{\mathsf T}\boldsymbol\phi(x), \qquad \boldsymbol\phi(x) = \left(x^0, x^1, \dots, x^M\right)^{\mathsf T}$$

$M$ is the order, and it is the only knob that controls how complicated the model is allowed to be. Stack one $\boldsymbol\phi(x_i)^{\mathsf T}$ per data point and you get the design matrix $\boldsymbol\Phi$, which is $N \times (M+1)$.

In code that is a Vandermonde matrix built one row at a time, from PlotFeatures.py:5-16:

def createSQM(x,M):
    vfunc = np.vectorize(lambda x,i: x**i)
    ilist=np.arange(M+1)
    xlist=np.repeat(x, M+1)
    return vfunc(xlist,ilist)

def getXMatrix(Xlist,M):
    out = np.empty((0,M+1))
    for x in Xlist:
        out = np.append(out,createSQM(x,M).reshape((1,M+1)),axis=0)

    return out

One detail that matters for reproducing anything: the dataset file contains only the targets. The inputs are reconstructed, at PlotFeatures.py:35, as $N$ points evenly spaced on $[0, 1]$: np.arange(0, 1+1.0/len(t), 1.0/(len(t)-1.0)). The widget below does the same thing.

Least squares, in closed form

Minimise the sum of squared residuals:

$$E(\mathbf{w}) \;=\; \frac{1}{2}\sum_{i=1}^{N}\left(t_i - \mathbf{w}^{\mathsf T}\boldsymbol\phi(x_i)\right)^2$$

Differentiate with respect to $\mathbf{w}$, set the gradient to zero,

$$\frac{\partial E(\mathbf{w})}{\partial \mathbf{w}} \;=\; -\sum_{i=1}^{N}\left(t_i - \mathbf{w}^{\mathsf T}\boldsymbol\phi(x_i)\right)\boldsymbol\phi(x_i) \;=\; 0$$

and you land on the normal equations, PRML 3.15:

$$\boldsymbol\Phi^{\mathsf T}\boldsymbol\Phi\,\mathbf{w} \;=\; \boldsymbol\Phi^{\mathsf T}\mathbf{t}$$

No iteration, no learning rate, no epochs. The entire fitting procedure is three lines, at PlotFeatures.py:21-28:

def Error(W,x,t):
    return 0.5*np.sum(np.square(np.dot(x,W)-t))

def getLogLiErr(W,x,t,Beta):
    return 0.5*Beta*np.sum(np.square(np.dot(x,W)-t))+0.5*len(t)*np.log(Beta/(2.0*np.pi))

def getOptimalW(x,t):
    return np.linalg.solve(np.dot(x.T,x),np.dot(x.T,t))

getOptimalW forms the Gram matrix $\boldsymbol\Phi^{\mathsf T}\boldsymbol\Phi$ explicitly and hands it to np.linalg.solve, which is an LU factorisation with partial pivoting, not pinv, not a QR least-squares solve, and with no regularisation term anywhere. That choice comes back to bite later in this post.

Maximum likelihood gives the identical answer

Now do it properly. Assume the target at each input is the model output plus Gaussian noise of precision $\beta$:

$$p(\mathbf{t}\mid\mathbf{X},\mathbf{w},\beta) \;=\; \prod_{i=1}^{N}\mathcal{N}\!\left(t_i \mid \mathbf{w}^{\mathsf T}\boldsymbol\phi(x_i),\, \beta^{-1}\right)$$

Take the log:

$$\ln p(\mathbf{t}\mid\mathbf{X},\mathbf{w},\beta) \;=\; \frac{N}{2}\ln\beta \;-\; \frac{N}{2}\ln(2\pi) \;-\; \frac{\beta}{2}\sum_{i=1}^{N}\left(t_i - y(x_i,\mathbf{w})\right)^2$$

The only term that depends on $\mathbf{w}$ is the last one, and it is $-\beta$ times the least-squares error. Maximising the likelihood over $\mathbf{w}$ is therefore exactly minimising $E(\mathbf{w})$: same normal equations, same weights, same curve.

This is the sentence I would put on the wall of every introductory course. Least squares is not a sensible-looking heuristic that people happen to use. It is the maximum likelihood estimator for a linear model with independent Gaussian noise of constant variance. Every time you reach for it you are asserting that noise model, whether or not you say so out loud. If your noise is heavy-tailed or heteroscedastic, least squares is not “still a reasonable choice”, it is the right answer to a question you did not mean to ask.

The report plots both, side by side, for order 4. They are the same plot, and that is the entire point of the figure:

Crank it to nine

Ten data points. An order-9 polynomial has ten coefficients. From the report:

The data has 10 degrees of freedom, all of which can be accounted for by the parametric equation. Due to this, the best fit for the data is one that goes through all the points. This has a very low error function but often does not generalise well to new data.

The curve passes through every point and, in doing so, spikes to $+2.8$ just past the left edge of the data and dives to $-1.6$ near the right. Recomputing the fit, the largest coefficient at $M=9$ is about $4.5\times10^{5}$; at $M=4$ it is $29$. Huge opposing coefficients cancelling each other to five significant figures is the numerical signature of this failure, and it is why the ninth-order Gram matrix is nearly singular.

Meanwhile the training error does what training error always does:

$M$	$E(\mathbf{w})$	RMS residual
0	1.771	0.595
1	0.890	0.422
2	0.885	0.421
3	0.201	0.200
4	0.196	0.198
5	0.163	0.181
6	0.163	0.181
7	0.063	0.113
8	0.018	0.060
9	$1.2\times10^{-8}$	$5\times10^{-5}$

The plateaus are the interesting bit. Nothing happens between $M=1$ and $M=2$, and nothing happens between $M=5$ and $M=6$, because the underlying function is odd: $\sin(x) = x - x^3/3! + x^5/5! - \dots$ has no even-order terms, so the even coefficients have almost nothing to explain. Adding one buys you essentially nothing. The report makes this argument later, to justify why the Bayesian evidence peaks at $M=3$; it is already visible here, in a plot that is otherwise just a monotone slide toward zero.

And monotone is the problem. $E(\mathbf{w})$ can never increase with $M$, because the order-$M$ model space is contained in the order-$(M{+}1)$ one: the larger model can always set the new coefficient to zero and match. A quantity that cannot go up is not a model-selection criterion. It is a measurement of how much freedom you gave the model.

The sting: $\beta$ does not survive the argmax

Go back to the log-likelihood. The precision $\beta$ is right there in the expression, and it drops out of the maximisation over $\mathbf{w}$ completely, it multiplies a term whose minimiser it cannot move. So maximum likelihood over the weights hands you a curve and tells you nothing whatsoever about how much to trust it. You have to supply $\beta$ yourself, and here it was simply asserted: the report assumes the data was generated with $\beta = 11.1$ (a noise standard deviation of about $0.30$).

You can see the consequence in the repo’s second error plot. getLogLiErr is called with $\beta = 11.1$ and its output is plotted against $M$:

It is the same curve. Look at the function again: it returns $\beta \cdot E(\mathbf{w}) + \tfrac{N}{2}\ln(\beta/2\pi)$, which for fixed $N$ and $\beta$ is $E(\mathbf{w})$ multiplied by $11.1$ and shifted by $2.85$. It cannot distinguish anything the first plot could not. That was the honest lesson of the whole exercise: I had computed a likelihood, felt like I had done something more principled than curve fitting, and produced a rescaled copy of the plot I already had.

Do it yourself

Everything above is one solve of a small linear system, so it runs comfortably in the browser at pointer speed. Click on the plot to add a point, drag one to move it, shift-click to delete it. The slider is $M$.

The moment worth chasing: set $M$ to 9, then add an eleventh point anywhere. Eleven points, ten parameters: the curve stops interpolating, $E(\mathbf{w})$ jumps off zero, and the sparkline’s right-hand end lifts off the floor. Nothing about the model changed. Only the arithmetic of degrees of freedom did.

With JavaScript on, this becomes a canvas you can click to place points, with a slider for the polynomial order, live readouts of E(w) and the log-likelihood error, and a log-scale sparkline of error against M.

The second preset is the repo’s other dataset, the 80 samples from Dataset_2.txt used later in the report for Bayesian model comparison. Load it and the picture inverts: with eighty points there is nothing for a ninth-order polynomial to interpolate, the error-vs-M sparkline flattens out after $M=3$, and every order from three upward draws more or less the same sine. Overfitting is not a property of the model. It is a property of the ratio.

The widget is a direct port of PlotFeatures.py, with one deliberate departure. NumPy’s np.linalg.solve on a ninth-order Gram matrix built from arbitrary user-placed points will happily return garbage or blow up, so the widget factorises with Cholesky and adds an escalating ridge to the diagonal until the factorisation succeeds, then reports the ridge it needed. Drag two points close together at $M=9$ and you can watch it kick in. The original script never had to deal with this because its ten inputs were always evenly spaced; that is not a virtue of the code, just of the data.

Why there is no WASM here

Several posts in this series push work into Rust compiled to WebAssembly because the alternative is a page that stutters. This is not one of them, and it is worth being explicit about where the line falls. The heaviest thing this widget does is form a $10\times10$ Gram matrix and Cholesky-factorise it, a few thousand floating-point operations, ten times over for the sweep across all orders, on every pointer move. That is microseconds. Reaching for WASM here would buy nothing measurable and cost a build step, a fetch and a memory boundary. The interesting engineering question is never “is this fast in principle”, it is “is this the part that is slow”.

Next

Least squares gave a point estimate of $\mathbf{w}$ and no error bars, and $\beta$ fell out of the maths without leaving a trace. The next post puts a prior on the weights instead and gets a predictive variance out of the same ten points, including error bars that fan out where the data is missing, which is the single most persuasive figure in this repo. That is the Bayesian treatment, and it is where the overfitting problem starts to solve itself.