When fitting a non-linear model using linear regression, we typically generate new features using non-linear functions. We also know that any function, in theory, can be approximated by a sufficiently high degree polynomial. This result is known as Weierstrass approximation theorem. But many blogs, papers, and even books tell us that high polynomials should be avoided. They tend to oscilate and overfit, and regularization doesn’t help! They even scare us with images, such as the one below, when the polynomial fit using the data points (in red) is far away from the true function (in blue):

It turns out that it’s just a MYTH. There’s nothing inherently wrong with high degree polynomials, and in contrast to what is typically taught, high degree polynomials are easily controlled using standard ML tools, like regularization. The source of the myth stems mainly from two misconceptions about polynomials that we will explore here. In fact, not only they are great non-linear features, certain representations also provide us with powerful control over the shape of the function we wish to learn.

A colab notebook with the code for reproducing the above results is available here.

Vladimir Vapnik, in his famous book “The Nature of Statistical Learning Theory” which is cited more than 100,000 times as of today, coined the approximation vs. estimation balance. The approximation power of a model is its ability to represent the “reality” we would like to learn. Typically, approximation power increases with the complexity of the model – more parameters mean more power to represent any function to arbitrary precision. Polynomials are no different – higher degree polynomials can represent functions to higher accuracy. However, more parameters make it difficult to estimate these parameters from the data.

Indeed, higher degree polynomials have a higher capacity to approximate arbitrary functions. And since they have more coefficients, these coefficients are harder to estimate from data. But how does it differ from other non-linear features, such as the well-known radial basis functions? Why do polynomials have such a bad reputation? Are they truly hard to estimate from data?

It turns out that the primary source is the standard polynomial basis for n-degree polynomials (mathbb{E}_n = {1, x, x^2, …, x^n}). Indeed, any degree (n) polynomial can be written as a linear combination of these functions:

[alpha_0 cdot 1 + alpha_1 cdot x + alpha_2 cdot x^2 + cdots + alpha_n x^n]

But the standard basis (mathbb{B}_n) is awful for estimating polynomials from data. In this post we will explore other ways to represent polynomials that are appropriate for machine learning, and are readily available in standard Python packages. We note, that one advantage of polynomials over other non-linear feature bases is that the only hyperparameter is their degree. There is no “kernel width”, like in radial basis functions¹.

The second source of their bad reputation is misunderstanding of Weierstrass’ approximation theorem. It’s usually cited as “polynomials can approximate arbitrary continuous functions”. But that’s not entrely true. They can approximate arbitrary continuous functions in an interval. This means that when using polynomial features, the data must be normalized to lie in an interval. It can be done using min-max scaling, computing empirical quantiles, or passing the feature through a sigmoid. But we should avoid the use of polynomials on raw un-normalized features.

In this post we will demonstrate fitting the function

[f(x)=sin(8 pi x) / exp(x)+x]

on the interval ([0, 1]) by fitting to (m=30) samples corrupted by Gaussian noise. The following code implements the function and generates samples:

import numpy as np

def true_func(x):
  return np.sin(8 * np.pi * x) / np.exp(x) + x

m = 30
sigma = 0.1

# generate features
np.random.seed(42)
X = np.random.rand(m)
y = true_func(X) + sigma * np.random.randn(m)

For function plotting, we will use uniformly-spaced points in ([0, 1]). The following code plots the true function and the sample points:

import matplotlib.pyplot as plt

plt_xs = np.linspace(0, 1, 1000)
plt.scatter(X.ravel(), y.ravel())
plt.plot(plt_xs, true_func(plt_xs), 'blue')
plt.show()

Now let’s fit a polynomial to the sampled points using the standard basis. Namely, we’re given the set of noisy points ({ (x_i, y_i) }_{i=1}^m), and we need to find the coefficients (alpha_0, dots, alpha_n) that minimize:

[sum_{i=1}^m (alpha_0 + alpha_1 x_i + dots + alpha_n x_i^n – y_i)^2]

As expected, this is readily accomplished by transforming each sample (x_i) to a vector of features (1, x_i, dots, x_i^n), and fitting a linear regression model to the resulting features. Fortunately, NumPy has the numpy.polynomial.polynomial.polyvanderfunction. It takes a vector containing (x_1, dots, x_m) and produces the matrix

[begin{pmatrix}
1 & x_1 & x_1^2 & dots & x_1^n \
1 & x_2 & x_2^2 & dots & x_2^n \
vdots & vdots & vdots & ddots & vdots \
1 & x_m & x_m^2 & dots & x_m^n \
end{pmatrix}]

The name of the function comes from the name of the matrix – the Vandermonde matrix. Let’s use it to fit a polynomial of degree (n=50).

from sklearn.linear_model import LinearRegression
import numpy.polynomial.polynomial as poly

n = 50
model = LinearRegression(fit_intercept=False)
model.fit(poly.polyvander(X, deg=n), y)

The reason we use fit_intercept=False is because the ‘intercept’ is provided by the first column of the Vandermonde matrix. Now we can plot the function we just fit:

plt.scatter(X.ravel(), y.ravel())                                    # plot the samples
plt.plot(plt_xs, true_func(plt_xs), 'blue')                          # plot the true function
plt.plot(plt_xs, model.predict(poly.polyvander(plt_xs, deg=n)), 'r') # plot the fit model
plt.ylim([-5, 5])
plt.show()

As expected, we got the “scary” image from the beginning of this post. Indeed, the standard basis is awful for model fitting! We hope that regularization provides a remedy, but it does not. Maybe adding some L2 regularization helps? Let’s use the Ridge class from the sklearn.linear_model package to fit an L2 regularized model:

from sklearn.linear_model import Ridge

reg_coef = 1e-7
model = Ridge(fit_intercept=False, alpha=reg_coef)
model.fit(poly.polyvander(X, deg=n), y)

plt.scatter(X.ravel(), y.ravel())                                    # plot the samples
plt.plot(plt_xs, true_func(plt_xs), 'blue')                          # plot the true function
plt.plot(plt_xs,