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

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)

• Although overfitting can be detected by poor prediction beyond the support of the dataset, overfitting is not really the same phenomenon as oscillation. After all, if the underlying phenomenon truly follows a high-degree polynomial (plus some error) and you fit exactly the correct degree, then you will not have overfit, even if that degree is an extremely high one. Thus the question seems to be a matter of how many parameters one can spend on modeling a dataset, rather than one about polynomials per se.
– whuber
Jan 13 at 19:01
• @Cagdas I'm sorry, but I can't make sense of any of that comment: I can't tell what you refer to by "this," what "information content" might mean; or what exactly is "catching up" as data points are added. It sounds like you might have some kind of flexible model in mind in which the polynomial degree is estimated as well as the polynomial, but I just can't tell.
– whuber
Jan 13 at 22:10
• BTW, any good answer will explain why polynomials differ in behavior from, say, piecewise splines that don't exhibit the same kind of radical behavior between data points. This comes down to the fact that polynomials are non-local: the behavior of a polynomial at a point $x$ is strongly controlled by the values you force it to have at points far from $x.$ In this sense polynomials cannot be "flabby" the way arbitrary continuous (or even infinitely differentiable) functions can be. The correct characterization, IMO, is in terms of the properties of analytic functions.
– whuber
Jan 13 at 23:14
• it's another poorly written blog on regression in the question. if you want to learn about regressions, stay away from blogs by computer programmers. most of them have no clue about the subject, unfortunately, but they like to blog indeed. Jan 15 at 18:20
• what is also annoying is that the blogger stole a chart in question from Bishop's book, it's p.7 Fig 1.4. if you're interested. is it that hard to attribute someone's work? so, why then read this blogger's posts when you can simply read the Bishop's book on the subject? it's very simple and explains polynomial regression properly. avoid third parties, go straight to the source Jan 15 at 18:27

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.

Assume we are given some fixed $$x_i$$ points, as well as $$y$$ values sampled from $$y_i\sim f(x_i)+\epsilon_i$$ where $$f$$ is some smooth function and $$\epsilon$$ are iid samples from a unit normal distribution. The Lagrange interpolation then takes the form $$\hat{f}_y(x)=\sum_i l_i(x)y_i$$ where $$l_i$$ are the Lagrange polynomials of the $$x_i$$. We want to calculate $$Var_y(\hat{f}(x))$$, where the variance is with respect to different sampled values of the $$y$$s. By independence of the $$\epsilon_i$$s, we have

$$Var_y(\hat{f}(x))=Var_y(\sum_i l_i(x)(f(x_i)+\epsilon_i))=\sum_i l_i(x)^2$$

Now, consider for simplicity the case where the observation points are $$x_1,\dots, x_n=1,\dots, n$$. In general, the variance $$Var(\hat{f}(x))$$ is a function of the location $$x$$, but we will consider the integrated variance $$Var=\int_0^nVar(\hat{f}(x))dx$$.

Claim: $$Var$$ increases (at least) like $$O(e^{an}n^{-9/2})$$ for some $$a>0$$.

In general, we have the formula $$l_i(x)={\frac {l(x)}{l'(x_i)(x-x_i)}}$$ where $$l(x):=\prod_i (x-x_i)$$. We can compute the denominator: $$l'(x_j)=\prod_{i\neq j}(j-i)=(j-1)(j-2)...(1)(-1)...(-(n-j))=(-1)^{n-j}(j-1)! (n-j)!$$

So the variance at $$x$$ is $$l(x)^2 \sum_j {\frac 1 {(j-1)!^2(n-j)!^2(x-x_j)^2}}$$. Now, $$|x-x_j|<=n/2$$, so $${\frac 1 {(x-x_j)^2}}\geq 4/n^2$$

$$\begin{eqnarray*} \sum_j {\frac 1 {(j-1)!^2(n-j)!^2}}{\frac 1 {(x-x_j)^2}}&\geq &4n^{-2}n!^{-2}\sum_j {n\choose j}^2\\ & = &4n^{-2} n!^{-2}{2n\choose n}\\ &\sim& n^{-7/2} 4^n(e/n)^{2n} \end{eqnarray*}$$

where we used Stirling's formula in the last line.

So in order to get a lower bound $$\int_0^n Var(\hat{f}(x))dx$$, we need to integrate $$l(x)^2$$. Noting that $$l(x)^2=n^{2n}\prod_i(x'-i/n)^2$$, where $$x'=x/n$$, we have

$$l(x)^2=n^{2n}\prod_i (x'-i/n)^2=n^{2n}e^{2\sum_i \log |x'-i/n|}\sim n^{2n}e^{2n\int_0^1 \log |x'-y|dy}=n^{2n}e^{2n(-H(x')-1)}$$ where $$H$$ is the binary entropy function.So,

$$\int_0^n l(x)^2x\sim n^{2n-1}\int_0^1 e^{2n(-H(x')-1)} dx'\geq n^{2n-1}e^{-2C(\epsilon)n}\int_{1-\epsilon}^1 dx'=n^{2n-1}e^{-2Cn}/\epsilon$$

where $$C(\epsilon)=1+H(1-\epsilon)$$ and $$0<\epsilon<1/2$$ is arbitrary.

The integrated variance $$\int_0^n Var(\hat{f}(x))dx$$ is thus at least $$O(e^{-2Cn}n^{2n-9/2}4^n(e/n)^{2n})=O(n^{-9/2}({\frac {2e}{e^C}})^{2n})$$

Now crucially, $$C(\epsilon)<1+\log(2)$$, because the binary entropy function is maximized at $$1/2$$. Therefore, $$2ee^{-C}>1$$, so we get the claimed form for the variance.

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()

• +1 Thank you for the excellent insights.
– whuber
Jan 15 at 23:31
• Thanks for the valuable answer. The Bernstein example however is misleading in my opinion. In the absense of rounding errors, the Bernstein and monomial bases would exactly lead to the same result (and thus show the oscillating behaviour seen in the OP). Thus, you seem to apply any form of regularization, which might be hidden in the algorithm. So could you please state how exactly you perform the fit? Jan 20 at 21:07
• @davidhigh You are correct that the Bernstein and monomial bases would both lead to the same result if they were both optimizing the same objective function. The Bernstein interpolant is defined as $\sum_i b_{i,n}(x)y_i$ where $b_i$ are the Bernstein basis polynomials. Equivalently, it is the solution to the optimization $argmin_{\hat{f}\in Poly_n}\sum_i E_{N\sim Bern(i/n,n)}(\hat{f}(i/n)-y_{N})^2$. So there is no explicit regularization in the Bernstein example. Instead, it is the solution to a different objective function. Jan 20 at 22:04

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.

• The first argument is good. The second, although intuitive, looks a little too hand-wavy, because all those oscillations could be pushed off into the Complex plane somewhere. Descarte's Rule of Signs might be of some help in making this argument stick. Regardless, the problem lies not in oscillation per se but in the high likelihood that they will have huge magnitudes.
– whuber
Jan 14 at 16:32
• @whuber: I agree with all of these criticisms--this is definitely not at the level of a publication quality argument yet! But, I figured a hand-wavy answer is still better than none at all; if someone else offers a better response, the rest of the community will surely vote it up to the top. Jan 14 at 18:49
• "I've never seen a really good general mathematical explanation or "proof" of this phenomenon" - Are you familiar with VC theory? en.wikipedia.org/wiki/Vapnik-Chervonenkis_theory
– usul
Jan 15 at 0:21
• @usul, could you please point to a specific result on that page which bears on the question? I am rather skeptical that there could be a proof that higher complexity model classes have higher variance, especially since it is not true (see e.g. pnas.org/content/116/32/15849, and also my answer to the OP) Jan 15 at 22:54
• @ Strachyra: Thank you so much for your answer! I have always been curious about these topics and have never been able to come across coherent answers and explanations on the history of these topics. For instance, over here I ask a question regarding how did researchers first realize that models might be overfitting and why they thought that adding a penalty term in a specific form to the optimization equation might mitigate this problem :stats.stackexchange.com/questions/560439/… . Can you please take a look at this if you have time?Thanks! Jan 16 at 7:40

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.

• this argument seems to imply that a neural network model with more parameters than data points will necessarily overfit, while there is a plethora of empiricial evidence to the contrary. Jan 15 at 23:23
• @MikeHawk I don't mean that all sets of functions with many degrees of freedom are prone to overfitting, just that some are. Jan 16 at 0:10
• fair enough, but then this does not answer OP's question, which was specifically about polynomials Jan 16 at 4:11
• @MikeHawk I think my argument applies for hypothesis classes where, for a given size of data set, it is "usually" possible to find "a close fit" and where the interpolation method prioritises doing so - which would include the most famous interpolation method for higher order polynomials. Jan 16 at 12:25

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)

• The orthogonal polynomials are merely a different basis: when you fit points with them, you will experience exactly the same problems described in the question. If you were to choose a different normalization of this basis, you could (say) limit all slopes to the range $[-1,1],$ thereby eliminating the appearance of any problem at the endpoints.
– whuber
Jan 13 at 22:12
• @whuber isn't it that making the slope at one end small will make it large on the other end? Those polynomials are functions with increasingly larger slope the further away you get from some point. Also, for an alternative basis of polynomials you might get the slope locally small, but it will be worse in other places, and still those alternative basis functions will be a composition of these 'ugly' looking ones. And as the degree increases it only get's worse. Jan 14 at 0:23
• No, one can reduce all the slopes simply by rescaling each polynomial. The problem does not lie in variation in the slopes per se; as I indicated in a comment to the question, it stems from the fact that forcing the polynomial to pass through one point $x$ can induce extraordinary changes in its value at other points $z$ far from $x.$ This is a form of "rigidity" that is best explained by the theory of entire Complex analytic functions (which are a natural generalization of polynomials and enjoy--or suffer from--many of the same properties).
– whuber
Jan 14 at 1:11
• +1 for your latest illustration and analysis: they are revealing. They suggest an examination of how the Mahalanobis distance varies as higher-degree polynomial terms are introduced into a hierarchy of models. I (strongly) suspect the large variations evident at the extremes of your explanatory variable can be attributed to a huge distance induced by the high powers. That distance, of course, is a key part of the formula for the prediction variance.
– whuber
Jan 14 at 16:42
• @ Sextus Empiricus: Thank you so much for your answer! I have always been curious about these topics and have never been able to come across coherent answers and explanations on the history of these topics. For instance, over here I ask a question regarding how did researchers first realize that models might be overfitting and why they thought that adding a penalty term in a specific form to the optimization equation might mitigate this problem :stats.stackexchange.com/questions/560439/… . Can you please take a look at this if you have time?Thanks! Jan 16 at 7:39

## 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.

• Your line of reasoning is not quite correct, I think. For example, a Fourier transformation is an expansion into polynomials, but trigonometric ones. Moreover, using 200.000 Chebyshev polynomials is likewise no big deal and will quickly converge, see for example the Chebfun-library. Same for other orthogonal polynomials. So, polynomials are good indeed, whereas monomials are unsuitable sometimed. Jan 21 at 12:57

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.

• +1. Thank you for sharing a genuinely new insight in this thread.
– whuber
Jan 18 at 15:53
• that is extremely interesting. but i would guess this is the exception, rather than the rule, no? that is, if the points are sampled from some distribution $\mu$, then (I think) you would still observe ringing, unless $\mu$ happens to be a sum of deltas at Chebyshev nodes (or maybe the limiting form of an arcsine distribution) Jan 20 at 20:08
• @MikeHawk yes, you need to tune the node locations to get this behaviour. Chebyshev are optimal in certain senses but I think any unimodal limiting density that accumulates 'quadratically' at the endpoints would show some similar benefits. Jan 29 at 19:04
• @Robert Mastragostino I actually believe that last claim is incorrect. In the limit of many points, we have $n^{-1}\log |l_j(x)|\to K(p)(x)-K(p)(x_j)$ where $l_j$ is the jth lagrange basis poly, $p$ is the limiting density, and $K(p)(x):=\int \log |x-y|p(y)dx$. Thus we should expect $l_j(x)$ to have exponentially large fluctuations, unless $K(p)$ happens to be a constant function. And there is a theorem by Tricomi that this is only true for the arcsine density. I worked out this argument in more detail in an addendum to my answer. Jan 30 at 14:46
• @MikeHawk Thanks for the update! Very illuminating. Jan 30 at 17:54

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.

• Could you indicate where in the Grunwald & Dawid paper this conclusion about "low order terms" is drawn? The paper does not discuss polynomials, nor does it use "orders" of terms explicitly, and at 70 pages it is too long to search manually.
– whuber
Jan 14 at 18:05
• Correct, this paper is not about polynomials, but rather, about a similar design choice when dealing with discrete data. Edited the post to clarify. Jan 14 at 21:25

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.

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.