How can we explain the "bad reputation" of higher-order polynomials?

Question

We all must have heard it by now - when we start learning about statistical models overfitting data, the first example we are often given is about "polynomial functions" (e.g., see the picture here):

We are warned that although higher-degree polynomials can fit training data quite well, they surely will overfit and generalize poorly to the test data.

Why does this happen? Is there a mathematical justification as to why (higher-degree) polynomial functions overfit the data? The closest explanation I could find online was something called "Runge's phenomenon", which suggests that higher-order polynomials tend to "oscillate" a lot - does this explain why polynomial functions are known to overfit data?

I understand that there is a whole field of "regularization" that tries to fix these overfitting problems (e.g., penalization can prevent a statistical model from "hugging" the data too closely) - but just using mathematical intuition, why are polynomials known to overfit the data?

In general, "functions" (e.g., the response variable you are trying to predict using machine learning algorithms) can be approximated using older methods like Fourier series, Taylor series and newer methods like neural networks. I believe that there are theorems that guarantee that Taylor series, polynomials and neural networks can "arbitrarily approximate" any function. Perhaps neural networks can promise smaller errors for simpler complexity?

But are there mathematical reasons behind higher-order polynomials (e.g., polynomial regression) being said to have a bad habit of overfitting, to the extent that they have become very unpopular? Is this solely explainable by Runge's phenomenon?

Reference:

Gelman, A. and Imbens, G. (2019) Why high order polynomials should not be used in regression discontinuity designs. Journal of Business and Economic Statistics 37(3), pp. 447-456. (An NBER working paper version is available here)

Answer 1

High degree polynomials do not overfit the data

This is a common misconception which is nonetheless found in many textbooks. In general, in order to specify a statistical model, it is necessary to specify both a hypothesis class and a fitting procedure. In order to define the Variance of the model ("variance" here in the sense of Bias-Variance Tradeoff), it is necessary to know how to fit the model, so claiming that a hypothesis class (like polynomials) can overfit the data is simply a category error.

To illustrate this point, consider that the Lagrange interpolation method is only one of several ways to fit a polynomial to a collection of points. Another method is to use Bernstein polynomials. Given $n$ observation points $x_i=i/n$ and $y_i=f(x_i)+\epsilon_i$, the Bernstein approximant is, like the Lagrange approximant, a polynomial of degree $n-1$ (i.e., as many degrees of freedom as there are data points). However, we can see that the fitted Bernstein polynomials do not exhibit the wild oscillations that the Lagrange polynomials do. Even more strikingly, the variance of the Bernstein fits actually decreases as we add more points (and thus increase the degree of the polynomials).

As you can see, the Bernstein polynomial does not pass through each point. This reflects the fact that it uses a different fitting procedure than the Lagrange polynomial. But both models have the exact same hypothesis class (the set of degree $n-1$ polynomials), which underscores that it is incoherent to talk about statistical properties of a hypothesis class without also specifying the loss function and fitting procedure.

As a second example, it is also possible to exactly interpolate the training data using a polynomial without wild oscillations, assuming we are working over complex numbers. If the observation points $x_i$ lie along the unit circle $x_i=e^{2\pi i\sqrt{-1}/N}$, then fitting a polynomial $y_i=\sum_d c_dx_i^d$ basically comes down to computing the inverse Fourier transform of $y_i$. Below I plot the complex polynomial interpolant as a function of the complex argument $i/N$. Moreover, it turns out that the variance of this interpolant is constant as a function of $n$.

Code for these simulations can be found at the end of this answer.

In summary, It is simply false that including higher degree polynomials will increase the variance of the model, or make it more prone to overfit

As a consequence, any explanation of why Lagrange interpolants do in fact overfit will need to use the actual structure of the model and fitting procedure, and cannot be based on general parameter-counting arguments.

Fortunately, it is not too hard to mathematically show why Lagrange interpolants do exhibit this overfitting. What is basically comes down to is that fitting the Lagrange interpolant requires inversion of a Vandermonde matrix, which are often ill-conditioned.

More formally, what we want to do is to show that the variance of the Lagrange interpolants grows very quickly as a function of $n$, the number of data points. It is not too hard to establish the following quantiative bound:

proposition Let $x_i=i/n$ for $i=1,\dots,n$. Let $\epsilon_i\sim N(0,1)$ be iid noise, and let $p_{\epsilon}(x)$ be the interpolating polynomial through the points $(x_i,y_i+\epsilon_i)$. Then the average variance $\int_0^1 Var_{\epsilon}p_{\epsilon}(x)dx$ grows at least as $O(e^{(1.5-\log 4)n})$.

proof of proposition

Let our interpolating polynomial be given by $p(x)=\sum_{d=0}^{n-1}w_d x^d$. Evidently p must satisfy $y_i=p(x_i)=\sum_{d=0}^{n-1}w_d x_i^d=(Vw)_i$ where the vandermonde matrix $V$ is defined by $V_{ij}=x_i^{j-1}$.

So we have the formula $w=V^{-1}y$ for the coefficients of the interpolating polynomial. As a warmup, we'll first consider the variance of the coefficient vector $w$, and then generalize the argument to analyze the variance of the polynomial $p(x)$.

coefficient variance

We want to compute the variance of $w$, that is $$Var(w):= E_{\epsilon} \|w-E_{\epsilon} w\|^2$$ If we write $w=V^{-1}y+V^{-1}\epsilon$, then the first term is constant (does not depend on epsilon), while the second term has zero mean. Because the variance is translation invariant, $$Var(w)=E_{\epsilon} \|V^{-1}\epsilon\|^2=E_{\epsilon} \epsilon^t (V^{-1})^tV^{-1}\epsilon=Tr( (V^{-1})^t V^{-1})=Tr((VV^t)^{-1})$$

Thus the variance of the coefficients is directly related to the eigenvalues of the matrix $VV^t$- the smaller these eigenvalues, the larger the variance. Now it becomes clear why the conditioning of $V$ is so important.

In the above discussion, the values of the $x_i$s was arbitrary, but let us now assume that they are evenly spaced in the unit interval: $x_i=i/n$. In this case, we can work out an explicit lower bound for the variance.

Firstly, note that $Tr(M)/n\geq det(M)^{1/n}$ for any symmetric positive definite matrix (this follows immediately from the AMGM inequality). In our case, $$Var(w)\geq n |det(V)|^{-2/n}$$

Now, the determinant of a vandermonde matrix has the well-known formula $det(V)=\prod_{i<j} x_j-x_i$. Under our definition of $x_i$, this is $$det(V)=\prod_{i<j}(j-i)/n=n^{-(n^2-n)/2}\prod_{i<j}j-i=n^{-(n^2-n)/2} \prod_{i=1}^{n-1} i!$$ (To see the last step, note that there are $n-1$ pairs $i,j$ such that $j-i=1$, there are $n-2$ pairs such that $j-i=2$ etc. )

The product of factorials is called the \textit{superfactorial}. Similarly to how the factorials are interpolated by the gamma function, the superfactorials are interpolated by an analytic function called the Barnes G function, $G$. More precisely, we have $G(n+1)= \prod_{i=1}^{n-1} i!$ when $n$ is a positive integer. Furthermore, there is an analogue of Stirling's formula: $$\log G(n+1)= .5*n^2\log n-.75*n^2+n\log \sqrt{2\pi}+O(log(n))$$ So the determinant goes like $$\log |Det(V)|\sim -(n^2-n)/2\log n+(n^2/2)\log n-3n^2/4+n\log \sqrt{2\pi}+...=-3n^2/4-(n/2)\log n+n\log \sqrt{2\pi}+...$$ So by the inequality for the variance, $$\log Var(w)\geq \log n -(2/n)\log |det(V)|=3n/2+2log n-\log 2\pi$$

Or said a bit more succinctly, $Var(w)=O(e^{3n/2})$, that is the variance grows exponentially.

prediction variance

To bound the variance of the model predictions, which is what we ultimately care about, it is necessary to modify the above argument only slightly.

As a bit of notation, define $pow(x)$ to be the vector $pow(x)_i=x^{i-1}$. Observing as before that $w$ decomposes into a constant $V^{-1}y$ and a mean-zero term $V^{-1}\epsilon$, we have that

$$Var(x)=Var_{\epsilon} \sum_d (V^{-1}\epsilon)_d x^d=Var_{\epsilon}pow(x)^t (V^{-1}\epsilon)=\|V^{-1}pow(x)\|^2$$ Integrating over $x$, \begin{eqnarray*} Var&=&\int Var(x)dx=\int Tr(V^{-1}pow(x)pow(x)^t(V^{-1})^t))dx\\ &=& Tr(V^{-1}(\int pow(x)pow(x)^tdx)(V^{-1})^t))\\ & = & Tr(V^{-1}H_n(V^{-1})^t) \end{eqnarray*} where $H_n$ is the Hilbert matrix given by $(H_n)_{ij}=1/(i+j-1)=\int_0^1 x^{i+j-2}dx$.

We want to bound this from below, arguing as before we can reduce this to a determinant calculation: $$Var=Tr(V^{-1}H_n(V^{-1})^t)\geq n det (V^{-1}H_n(V^{-1})^t)^{1/n}=n |det V|^{-2/n}|det H|^{1/n}$$

Now, the Hilbert determinant has a well-known asumptotic expansion: $$\log det (H)\sim -n^2 \log(4)+n\log 2\pi +...$$ Thus the dominant term in $|det(H)|^{1/n}$ will be $e^{-n \log 4}$. On the other hand, we saw above that $n |det(V)|^{-n/2}$ has dominant term $e^{3n/2}$. Thus we can conclude that $Var=O(e^{(3/2-\log 4)n})$. Note crucially that $3/2>\log 4$, so that the variance indeed increases exponentially. This concludes the proof.

Edit: analysis of Bernstein variance As an interesting comparison, it is not too hard to work out the variance for the Bernstein approximator. This polynomial is defined by $\hat{B}(x)=\sum_i b_{i,n}(x)y_i$, where $b_{i,n}$ are the Bernstein basis polynomials, as defined in the Wikipedia link. Arguing similarly as before, the variance with respect to different samplings of $y$ is $\sum_i b_{i,n}(x)^2$. As in the simulation, I'll average over $x\in [0,1]$. Evidently $\int_0^1 b_{i,n}(x)^2dx={n\choose i}^2B(2i+1,2n-2i+1)$, with $B$ the beta function. We have to sum this over $i$. I claim the sum comes to: $$\int_0^1 Var(\hat{B}(x))dx={\frac {4^n}{(2n+1){2n\choose n}}}$$

(proof is below)

This formula does indeed closely match the empirical results from the simulation (note that for direct comparison, you will have to multiply by $.1^2$, which was the noise variance used in the simulation). By Stirling's formula, the variance of the Bernstein approximator goes as $O(n^{-.5})$.

Proof of variance formula The first step is to use the identity $B(2i+1,2n-2i+1)=(2i)!(2n-2i)!/(2n)!$. So we get

\begin{eqnarray*} {n\choose i}^2 B(2i+1,2n-2i+1) & = & {\frac {(n!)^2}{(i!)^2(n-i)!^2}}{\frac {(2i)!(2n-2i)!}{(2n+1)!}}\\ & = & {\frac 1 {2n+1}}{2n\choose n}^{-1}{2i\choose i}{2n-2i\choose n-i} \end{eqnarray*}

Using the generating function $\sum_{i=0}^{\infty} {2i\choose i}x^i=(1-4x)^{-1/2}$, we see that $\sum_{i=0}^n {2i\choose i}{2n-2i\choose n-i}$ is the coefficient of $x^n$ in $((1-4x)^{-1/2})^2=(1-4x)^{-1}$, that is $4^n$. Summing over $i$ we get the claimed formula for the variance.

Edit: effect of distribution of x points on Lagrange variance As shown in the very interesting answer of Robert Mastragostino, the bad behavior of the Lagrange polynomial can be avoided through judicious choice of $x_i$. This raises the possibility that the lagrange polynomial does uncharacteristically badly for uniformly sampled points. In principle, given iid sample points $x_i\sim \mu$ it is possible that the lagrange polynomials behave sensibly for most choices of $\mu$, with the uniform distribution being just an ``unlucky choice". However, this is not the case, and it turns out that the Lagrange variance still grows exponentially for any choice of $\mu$, with the single exception of the arcsine distribution $d\mu(x)=\textbf{1}_{(-1,1)}\pi^{-1}(1-x^2)^{-1/2}dx$. Note that in Robert Mastragostino's answer, the $x_i$ points are chosen to be Chebyshev nodes. As $n\to\infty$, the density of these points converges to none other than the arcsine distribution. So in a sense, the situation in his answer is the only exception to the rule of exponential growth of the variance.

More precisely,take any continuous distribution $\mu(dx)=p(x)dx$ supported on $(-1,1)$ and consider an infinite sequence of points $\{(x_i,f(x_i)+\epsilon_i)\}_{i=1,\infty}$ where $x_i\sim \mu$ are iid, $f$ is a continuous function and $\epsilon_i$ are iid normal samples. Let $\hat{f_n}(x)$ denote the Lagrange interpolant fitted on the first $n$ points. And let $V_n=\int Var(\hat{f_n})(x)p(x)dx$ denote the corresponding average variance.

Claim: Assume that (1): $p(x)>0$ for all $x\in (-1,1)$ and (2): $\int p(x)^2\sqrt{1-x^2}dx<\infty$. Then $V_n$ grows at least like $O(e^{\epsilon n})$ for some $\epsilon>0$, unless $p$ is the arcsine distribution.

To prove this, note that, as before, $Var(\hat{f_n}(x))=\sum_{i=1}^nl_{i,n}^2$, where now $l_{i,n}$ denotes the ith Lagrange basis polynomial wrt the first $n$ points. Now,

\begin{eqnarray*} \log \int l_{i,n}^2(x)d\mu(x)&\geq& 2\int \log |l_{i,n}(x)|p(x)dx\\ & = & 2\int \sum_{j\neq i, j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx \end{eqnarray*} where we used Jensen's inequality in the first line, and the definition of $l_{i,n}$ in the second.

Let us introduce the following notation: \begin{eqnarray*} f_p(x)&=&\int \log|x-y|p(y)dy\\ c_p& = & \int \log |x-y|p(x)p(y)dxdy \end{eqnarray*}

Then we get

\begin{eqnarray*}\lim\inf n^{-1}\log \int l_{i,n}^2(x)p(x)dx&\geq& 2\lim\inf_n n^{-1}\int \sum_{j\neq i, j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx\\ & \geq & 2 \int \lim\inf_n n^{-1}\sum_{j\neq i, 1,j\leq n} \left(\log |x-x_j|-\log |x_i-x_j|\right)p(x)dx\\ & = & 2\int \left(f_p(x)-f_p(x_i)\right)p(x)dx\\ & = & 2(c_p-f_p(x_i)) \end{eqnarray*} using Fatou's lemma, and then the strong law of large numbers.

Note that $f_p'(x)=H_p(x)$, where $H_p(x)=\int {\frac {p(y)}{x-y}}dy$ is the Hilbert transform of $p$. In particular, $f_p$ is continuous. Now, suppose it is the case that $f_p(x)$ is non-constant. In this case, there exists $\epsilon>0$ as well as an interval $[a,b]$ such that $\inf_{x'\in [a,b]}c_p-f_p(x')\geq\epsilon/2$. By assumption (1), $\int_a^b p(x)dx>0$, so accordingly, $P(\exists i: x_i\in [a,b])=1$. We can reorder finitely many samples without changing the distribution, so WLOG we may assume that $x_1\in [a,b]$. As a consequence $2(c_p-f_p(x_1))\geq\epsilon$.

Combining with the above, we get
$$n^{-1}\log \int l_{1,n}^2(x)p(x)dx\geq \epsilon$$ for all $n$ sufficiently large. By simple rearrangement, this is equivalent to $\int l_{1,n}^2(x)p(x)dx\geq e^{n\epsilon}$ and since $V_n=\int \sum_{i=1}^n l_{i,n}^2(x)p(x)dx \geq \int l_{1,n}^2(x)p(x)dx$ we get the claimed exponential growth, assuming $f_p$ is non-constant.

So the last step is to show that $f_p$ is non-constant if $p$ is not the arcsine distribution, or equivalently, that the hilbert transform $H_p$ is not identically zero. Under assumption (2), this follows from a theorem of Tricomi. Thus we have proved the Claim.

Code for Bernstein polynomials

from scipy.special import binom
import numpy as np
import matplotlib.pyplot as plt
variances=[]
n_range=[5,10,20,50,100]
for n in n_range:
    
    preds=[]
    for _ in range(1000):
        xs=np.linspace(0,1,100)

        X=np.linspace(0,1,n+1)
        Y=np.sin(8*X)+np.random.randn(n+1)*.1
        nu=np.arange(n+1)
        bern=binom(n,nu)[:,None]*(xs[None,:])**(nu[:,None])*(1-xs[None,:])**(n-nu[:,None])*Y[:,None]

        pred=bern.sum(0)
        preds.append(pred)
    preds=np.array(preds)
    variances.append(preds.var(0).mean())
variances=np.array(variances)
plt.scatter(n_range,variances)
plt.xlabel('n')
plt.ylabel('variance')
plt.ylim(0,.005)
plt.title('bernstein polynomial variance')

Code for complex polynomials

n=64
xarg=np.linspace(0,2*np.pi,n+1,endpoint=True)[:n]
y=np.sin(2.5*xarg)+.25*np.random.randn(n)
xcomp=np.exp(1j*xarg)

xs=np.linspace(0,2*np.pi,200)
ys=np.sin(2.5*xs)

X=np.array([xcomp**i for i in range(len(xarg))])
w=np.linalg.solve(X,y)


#only need to go up to middle frequency, bc later ones are complex conjugates
interp=2*np.dot(w[1:n//2+1],np.array([np.exp(1j*xs*i) for i in range(1,n//2+1)]))+w2[0]

plt.plot(xs,interp.real,label='complex polynomial interp')
plt.plot(xs,ys,c='gray',linestyle='--',label='true function')
plt.scatter(xarg,y)
plt.xlabel('complex argument')
plt.legend()

Answer 2

Before I attempt to answer this, I'd like to just point out that what you are observing here (overfitting in a linear regression) is just a specific example of a more general phenomenon, bias-variance trade-off, which is also observed in more "modern" machine learning contexts as well as the "classical" setting of regression fitting. Your question could probably be expanded and rephrased more generally as something like, "why do high complexity models (e.g., those with larger numbers of adjustable parameters, such as higher degree polynomials) usually exhibit higher variance?" I've never seen a really good general mathematical explanation or "proof" of this phenomenon; most discussion that I've encountered seems to treat it as more of an empirical observation rather than theoretical result. This may be related to the fact that the field has a huge number of alternative model selection criteria currently in widespread use, indicating that there isn't really a good theoretical consensus yet on how to properly quantify a model's complexity (nor the closely related bias-variance trade-off) in the first place.

With those preliminary remarks out of the way, what do we mean by the term "overfitting" in the context of a statistical regression against a high degree polynomial? I submit that it may be broken down into two separate phenomena:

• As the degree of the polynomial in the regression increases, the resulting curve fits each data point ever more closely
• Simultaneously, the oscillations between the data points become more numerous and larger in amplitude--meaning that if we later acquire new data points by making additional observations, those new data points will generally tend to fit the previously fitted curves more and more poorly, as the degree of the fit model increases

How might we explain the above two points, mathematically?

For a polynomial curve fit model specifically, one way of understanding bullet one mathematically is to observe that a lower degree polynomial model is just a higher degree polynomial in which most of the fitted coefficients have been artificially pre-constrained to be equal to zero. For example, if we model a data set using a polynomial of degree 2: $$ f(x) = C_{0} + C_{1} x + C_{2} x^{2} $$ adjusting the parameters $C_{0}, C_{1}, C_{2}$ to result in the best possible fit, it is in some sense equivalent to saying that we have modeled it using a degree 100 polynomial: $$ f(x) = C_{0} + C_{1} x + C_{2} x^{2} + C_{3} x^{3} + ... + C_{100} x^{100} $$ in which the upper 98 coefficents have all been constrained such that: $$ C_{3} = C_{4} = ... C_{99} = C_{100} = 0 $$ Now, looking at it from that point of view, consider what happens when we add a degree to our polynomial model function, increasing from degree 2 to degree 3: it is in some sense like removing a constraint: where $C_{3}$ was previously forced to be $0$, now the fit algorithm is free to adjust it. Essentially, adding another adjustable parameter gives the fitter an additional "knob" that can be adjusted to allow the resulting fitted curve to more closely track the underlying data. Another crucial point to bear in mind: by removing the constraint that $C_{3} = 0$, you also allow the fitter a slightly wider latitude to adjust the original three parameters, $C_{0}, C_{1}, C_{2}$, across a wider range of potential values, because there are more opportunities for additional new adjustments to compensate or partially cancel each other out. So bullet number one (the phenomenon that higher degree polynomials more closely reproduce the observation data) may be understood as a simple result of the fact that fewer constraints and more "adjustable dials" means more opportunities to adjust all the parameter values "just so" in order to obtain a nearly perfect fit to the existing observed data.

But what about all of the wildly gyrating oscillations--how may we understand those, mathematically? Well, consider that the locations of the local minima / maxima for a polynomial of degree $n$ are obtained by solving the following equation: $$ \frac{d f(x)}{dx} = 0 $$ or equivalently, $$ C_{1} + 2 C_{2} x + 3 C_{3} x^{2} + ... + n C_{n} x^{n-1} = 0 $$ In general, a polynomial of degree $n$ may have up to $n-1$ unique local minima / maxima, because the algebraic equation above may have up to $n-1$ unique, real-valued (i.e., non-complex) roots. Because each local minimum / maximum in this case corresponds to one oscillation, it means that in the absence of a regularization constraint, increasing the degree of the fit model will naturally tend to increase the number of oscillations.

Additional bonus comment: if you program in R, there is an example script published with this question which facilitates experimenting with the practical effects of fitting real data with very high degree polynomials.

Answer 3

It isn't something special about higher-order polynomials: the same effect happens for other sets of functions with many degrees of freedom.

For example, let's call a function "special" if its graph consists of a horizontal segment, followed by a segment which slopes upwards at 45 degrees, followed by a straight line segment which is exactly twice as long as the previous segment but can be at any angle, followed by a straight line segment which can be of any angle and length, followed by a segment which slopes downwards at 60 degrees, followed by a quarter-circle arc of any size, followed by a segment which slopes upwards at 10 degrees, followed by a semicircular arc of any size, followed by another horizontal segment. When choosing a "special" function you have so many degrees of freedom that you likely can find one which closely fits your data. This is not an indication that special functions are a good modelling choice.

Similarly, with higher-order polynomials you have so many coefficients to play with that it just isn't very impressive that by picking them correctly you can achieve a close fit.

Answer 4

You have more problems than just Runge's phenomenon. Below is an example for fitting a tenth degree polynomial to 21 data points that follow the curve

$$y = sin(6\pi x^2)$$

• Runge's phenomenon: The black broken line is the least-squares fit to these 21 points if there is no noise. You see some larger error towards the edges, but when you only interpolate then it is not incredibly large.

• Variance and overfitting. On top of that phenomenon, you get that when we add noise, then the variance of the function in between data points becomes very large.

  In the image we show this by computing the standard deviation for the sample distribution of the estimate when there is white noise with variance of $\sigma = 0.2, 0.4 \text{ or } 0.6$

  In this case with fitting a 10h order polynomial you see that the error due to sample variation is relatively constant over all values for $x$. It is only for extrapolation that there is a large influence.

So this matter of overfitting does not seem to be much like Runge's phenomenon.

I have to explain the graphs a bit better. But, below is the R-code that allows to see what I did.

### create data
set.seed(1)
sigma = 0.2
### fine grained data
xs = seq(-0.1, 1.1, 0.001)
ys = sin(xs^2*6*pi)
### 21 data points
x = seq(0,1,0.05)
y = sin(x^2*6*pi)

### polynomials for modelling
M = as.data.frame(poly(x,10))
Ms = as.data.frame(predict(poly(x,10), xs))

### model of fit without noise
mod1 <- lm(y ~ ., data = M)

plot(xs,ys, type = "l", ylim = c(-2,2), lwd = 2,
     xlab = "x", ylab = "y", xlim = c(0,1))

### add confidence intervals at three levels of sigma
for (sig in 1:3) {
  pred = predict(mod1, newdata = Ms)+cor1
  ### the error of the coefficients is sigma * I because 
  ### (X^tX)^-1 is equal to the identity matrix 
  sigerr = sqrt(rowSums(Ms^2))*sigma*sig
  polygon(c(rev(xs), xs), c(rev(pred+sigerr),pred-sigerr),
          col = rgb(1,0,0,0.3), border = NA)
}

### true data
lines(xs,ys, lwd = 2)


### use this to simulate multiple fits
#for (i in 1:1000) {
#  q = y+rnorm(21,0,sigma)
#  mod <- lm(q ~ ., data = M)
#  lines(xs, predict(mod, newdata = Ms)+cor1, lty = 1, lwd = 1, col = rgb(1,0,0,0.01)) 
#}

### fit of model to 21 points
lines(xs, predict(mod1, newdata = Ms), lty = 2, lwd = 2)

## data points
points(x,y, pch = 21, col = 1, bg = 0)

title(expression(sin(6 * pi * x^2) * " with least squares fit to 21 points"))

legend(0,-1.2, c("true function with noisless data points",
                 "fit of noisless data points"), 
       lty = c(1,2), lwd = 2, pch = c(21,NA), pt.bg = c(0,0), cex = 0.8)

Answer 5

It's much worse than just overfitting.

The problems with polynomials don't become clear in examples with only 10 or 20 parameters, so I'll examine a function that we want to fit with 200,000 parameters, where 200,000 parameters really is the correct number - we aren't overfitting.

Our function is a 10 second audio clip, with pressure sampled 20,000 times per second. If we use linear interpolation, then we can interpret our function as a linear combination of 200,000 little bumps. All these bumps behave basically the same, and if we need more parameters we can just make more. We are dealing with lots of parameters, but adding 10 more parameters isn't that much harder than adding the first 10. A 5 cent computer embedded in a happy meal toy can handle this to play back a cow mooing.

Another classic approach would be to represent our function as a linear combination of sine waves. Finding 200,000 reasonable sine waves is easy enough, and using a fourier transform we can do this computation on any phone processor, without much fuss. sin(2x) isn't really that different from sin(200,000x.)

Now, if we want to model our audio clip as a polynomial, then we have to model it as a linear combination of 200,000 monomials. While finding 200,000 well behaved bumps was easy, and finding 200,000 chill sine waves was a breeze, finding 200,000 well behaved monomials is much harder! $1$ and $x$ were easy enough to work with. We start to run into trouble with $x^{20}$, but with some careful work we could wrangle it into representing audio. By the time we are dealing with $x^{1000}$ we are pulling our hair out- when we put in 0, 1, 2, it puts out 0, 1, 1.071509e+301, which barely even fits in a 64 bit floating point number. $x^{20,000}$ is not a reasonable function to work with, and yet when we resort to recruiting it, we have 180,000 functions still to go. $x^{200,000}$ is right out.

tldr: When representing a complicated thing as a linear combination of building blocks, you need a large collection of reasonable building blocks. There aren't very many reasonable monomials.

Answer 6

The ringing is an artifact of using uniformly spaced points, because the lagrange polynomials for this spacing are not tightly concentrated around the points they are trying to fit. E.g. for 11 evenly-spaced points on the interval [0,1], here is the degree 10 polynomial that is used to fit the value at x=0.4 (is zero at all other sample points):

Clearly nudging the value at x=0.4 will wildly change the fitted function at unrelated locations. As the Chebfun package (and Trefethen's other work) shows, this problem disappears when using well chosen sampling points over the interval, based on the roots of Chebyshev polynomials. The lagrange polynomials (i.e. basis functions) in this case naturally look closer to smooth kernels. With this choice of fitting points/polynomial basis the polynomial attempting to fit x=0.4 is in fact maximized at this location.

The tradeoff is that you have to fix your desired fitting interval in advance. As we can see in the example the 'concentration' property drops off very sharply outside of [0,1].

So the wild overfitting isn't really about polynomials: the issue is that on evenly spaced points the mapping between curve-distance over the interval and sample-distance on the sampled points is near-degenerate. While we intuitively think of high-degree polynomials as smooth interpolants, this is not automatically true of high-order polynomials and you have forgotten to inform your fitting procedure about what a normal measure of function-distance looks like. This is completely fixable by switching basis, and once this is done you can do very reliable numerical interpolation with 100+ degree polynomials. Of course, this solution is better suited to working numerically with in-principle continuous quantities, rather than for discrete data where we don't get to choose our sample points.

Answer 7

Is there any mathematical justification as to why (higher degree) polynomial functions overfit the data?

Sure. As others mentioned, particularly stachyra and fblundun, it's about the complexity of the hypothesis class relative to the amount of data you have. A highly complex model will always find a way to explain a small amount of data, regardless of whether that explanation generalizes correctly. A simple model won't be able to fit the training data unless it actually explains the underlying relationship.

Imagine that my data distribution looks like this: for each $x$, I pick $y$ uniformly at random from $\{0,1\}$. No model in the world will be able to predict the next data point. It's totally random.

But suppose you plan to draw $10$ training data points and fit a degree-$10$ polynomial. I can already tell you in advance what will happen: you will be able to fit the data perfectly. In other words, your approach can't tell whether you're going to generalize well or not. It always has zero training error even when the next prediction will be very bad.

Whereas with a degree-$3$ polynomial, you will immediately notice a poor fit and conclude that you will not generalize. A simple model can detect whether or not it is fitting the data. A complex one will always fit the data regardless of how it's generated.

---

Intuition like this is formalized by VC-dimension, a measure of complexity of hypothesis classes for binary classification (but there are versions for regression as well, psuedo-dimension). The theory promises that for a simple model class, if we draw relatively few data points, then the fit to the training data is representative of the fit to the actual generative model. Whereas for a more complex model class, it may fit the training data significantly better than the test data.

Answer 8

Using more terms means more degrees of freedom, hence more overfitting, but the real question is -- why are high order terms like $x^5$ so bad?

If you had enough data to estimate one parameter $a$, and no special knowledge about the domain, a practitioner would always start with $ax$ rather than $a x^5$ in their model exploration. Is it purely empirical, or is there some mathematical way to justify this choice?

When fitting discrete data you can show that low order terms are preferable to high order terms using worst case analysis. https://arxiv.org/pdf/math/0410076.pdf

For random binary vectors of the form $\langle x_0,x_1,\ldots,x_n\rangle$, you can show that the best model relying on $k$ features of the form $\{x_i\}$ (first order) will result in a better expected fit when true generating distribution is picked adversarially, after your choice of features is revealed, than the best model relying on $k$ features of the form $\{x_i x_j\}$ (second order). I'm not aware of similar results for real valued $x$. The big obstacle is that for real-valued $x$ is that there is no longer agreement as to which measure to use over this space. Different choices of measure lead to different conclusions.

Perhaps there is some universality principle which suggests that low order terms are preferable for a typical case, which more relevant for real world, rather than "adversarial case" which is what is analysed in the paper.

Answer 9

The blog post that you are discussing here and the chart have little to do with polynomial regressions per se. The author simply used polynomials to demonstrate the overfitting idea: i.e. fitting to the noise. When you leave no degrees of freedom, your fit become very rigid, i.e. very sensitive to both errors in y's and to a choice of x's in your sample. You could use any number of other basis functions to get to about the same result, it didn't have to be polynomials.

So, if you have 10 observation and use 9th order polynomial, then you leave no degrees of freedom. This polynomial will have 10 coefficients for your 10 sample points. So, your fitted curve will have to go through every point in your sample. This means that every measurement error in every Y will be in your model coefficients. You packed all the errors into your coefficients, then inevitably out of sample fit will be awful, or as ML people say "model won't generalize."

Again, this will happen with any model, not only polynomial regression. This doesn't mean that high degree polynomials don't have issues. They do, and some of them are real and some of them are due to lack of knowledge of people mis-using them, but this example is not a demonstration of these issues.