Polynomial regression

Examples

A polynomial model is still linear in its parameters, so we can fit it with least squares. Only the design matrix changes.

\[ \hat y = \theta_0 + \theta_1 x + \theta_2 x^2 + \cdots + \theta_d x^d \]

Fitting

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures

rng = np.random.default_rng(4)

def f(x):
    return np.sin(2 * x) + 2 * (x - 0.5) ** 3 - 0.5 * x

n, sigma = 30, 0.15
x = rng.uniform(0, 1, n)
y = f(x) + rng.normal(0, sigma, n)
X = x.reshape(-1, 1)

model = make_pipeline(PolynomialFeatures(4), LinearRegression()).fit(X, y)
print("training MSE:", round(float(np.mean((model.predict(X) - y) ** 2)), 4))

1

PolynomialFeatures builds the design matrix and LinearRegression fits it. Putting them in a Pipeline means the transformation is refitted whenever the model is refitted, which matters as soon as we cross-validate.

training MSE: 0.0142

The degree controls the flexibility. The lower panel shows the training error and the test error for every degree, with the current one marked. Increase \(n\) and watch the best degree move to the right.

Degree 1 is too rigid. It misses the curvature. Degree 12 goes through almost every point, including the noise. Degree 4 is about right.

Fitting many training sets

The clearest way to see variance is to fit the same model to several training sets from the same generator.

Twelve fits to twelve training sets from the same generator. A rigid model gives lines that are close together and far from the truth. A flexible one gives lines that straddle the truth but scatter widely.

The lines on the left are close together and far from the truth. The lines on the right straddle the truth but are far from each other.

Training error is not enough

The training error keeps falling. The test error has a minimum, and it never goes below the irreducible error.

The curse of dimensionality

With one input and degree \(d\) we have \(d+1\) parameters. With \(p\) inputs and degree \(d\) the number of terms grows like \(p^d\). A degree 3 polynomial in 100 inputs has more than 170000 terms.

The same problem appears for any method that needs to fill the input space with data.

The same idea for classification. The true boundary is a circle, so a linear model cannot express it at any noise level. Degree 2 can.

In high dimensions there is no such thing as a nearby point. This is the main reason why methods that make stronger assumptions, like neural networks or trees, do better in high dimensions than flexible local methods.