Real relationships are rarely straight lines. Crop yield rises with fertiliser and then falls; demand follows seasonal cycles. Does that mean we must abandon linear regression? No. The trick is to transform the inputs into new features — basis functions — and then fit a linear model on those. This simple idea is also our first clear view of overfitting.
Basis function expansion#
Replace the raw input $x$ with a vector of features $\boldsymbol{\phi}(x) = [\phi_0(x), \phi_1(x), \dots, \phi_M(x)]$ and fit
The model is still linear in $\mathbf{w}$, so all the machinery of least squares — closed form, convexity, fast solvers — still applies. Only the design matrix changes: $\boldsymbol{\Phi}_{ij} = \phi_j(x_i)$.
Common basis functions#
| Basis | Form | Character |
|---|---|---|
| Polynomial | $\phi_j(x) = x^j$ | Global; unstable at high degree and at edges |
| Radial basis (Gaussian) | $\phi_j(x) = \exp\left(-\frac{(x - c_j)^2}{2s^2}\right)$ | Local bumps centred at $c_j$ |
| Splines | Piecewise polynomials joined smoothly at knots | Flexible, stable, widely used in statistics |
| Fourier | $\sin(kx), \cos(kx)$ | Periodic patterns (seasonality) |
| Step / binning | $\mathbb{1}[x \in \text{bin}_j]$ | Piecewise constant; simple, interpretable |
With several inputs, polynomial features include interaction terms such as $x_1x_2$. The number of features grows combinatorially: degree-$p$ polynomials in $d$ variables have $\binom{d + p}{p}$ terms.
Watching overfitting happen#
Fit polynomials of increasing degree to 15 noisy samples of a sine wave:
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error
rng = np.random.default_rng(1)
f = lambda x: np.sin(2 * np.pi * x)
x_tr = np.sort(rng.random(15)); y_tr = f(x_tr) + rng.normal(0, 0.2, 15)
x_te = np.linspace(0, 1, 200); y_te = f(x_te) + rng.normal(0, 0.2, 200)
for degree in [1, 3, 5, 9, 14]:
model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
model.fit(x_tr[:, None], y_tr)
tr = mean_squared_error(y_tr, model.predict(x_tr[:, None]))
te = mean_squared_error(y_te, model.predict(x_te[:, None]))
print(f"degree {degree:>2}: train MSE {tr:.4f} test MSE {te:.4f}")A typical result:
| Degree | Training MSE | Test MSE | Diagnosis |
|---|---|---|---|
| 1 | high | high | Underfitting — a line cannot bend |
| 3 | low | low | Good fit |
| 9 | very low | higher | Starting to overfit |
| 14 | ≈ 0 | enormous | Overfitting — passes through every point, wild oscillations |
With 15 points and 15 coefficients (degree 14), the polynomial interpolates the data exactly, including its noise. Plot it and you will see huge swings between training points — especially near the edges (Runge's phenomenon). The coefficients also become enormous, a sign we will exploit with regularisation.
Choosing the degree#
from sklearn.model_selection import cross_val_score
for degree in range(1, 12):
model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
score = -cross_val_score(model, x_tr[:, None], y_tr, cv=5,
scoring="neg_mean_squared_error").mean()
print(degree, round(score, 4))Pick the degree with the lowest cross-validated error — or, better, keep a flexible basis and control complexity with regularisation (next lecture), which is usually more stable than choosing a discrete degree.
Splines: the practical choice#
High-degree global polynomials are unstable. Splines fit low-degree polynomials (typically cubic) on intervals between knots, constrained to join smoothly. They give flexible, well-behaved curves and underlie Generalised Additive Models (GAMs):
where each $f_j$ is a smooth spline. GAMs capture non-linear effects while remaining interpretable — you can plot each $f_j$ — which makes them popular in medicine and public policy.
from sklearn.preprocessing import SplineTransformer
spline_model = make_pipeline(SplineTransformer(n_knots=6, degree=3), LinearRegression())
spline_model.fit(x_tr[:, None], y_tr)
print("spline test MSE:", mean_squared_error(y_te, spline_model.predict(x_te[:, None])))From basis functions to kernels and neural networks#
Basis expansion raises a question: which basis should we choose? Two powerful answers come later in the course:
- Kernel methods implicitly use infinitely many basis functions while computing only inner products.
- Neural networks learn the basis functions from data: each hidden unit is an adaptive basis function, and the output layer is a linear model on top of them.
In this sense, a neural network is "linear regression on learned features".