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)
 A: 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.
A: 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)

A: 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.
A: 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.
A: 
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.
A: 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.
A: 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.
A: 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.
