Polynomial regression fits a curve by adding powers (and, with several inputs, products) of the original features and then running ordinary linear regression on the expanded set. The model y = β0 + β1x + β2x2 + … + βdxd is curved in x but still linear in the parameters, so all of the least-squares machinery — closed-form fit, standard errors, regularization — carries over unchanged. It is the textbook illustration of the bias–variance trade-off: the degree d directly dials model capacity.
Part of Regression within Supervised Learning.
The three panels fit the same 14 noisy points. A straight line (d = 1) cannot bend, so it misses the trend everywhere: high bias. A cubic follows the underlying wave. A degree-12 polynomial has almost as many coefficients as there are points, so it threads through nearly every one and oscillates wildly between and beyond them: high variance. The training error is lowest for degree 12, yet its predictions on new data are the worst of the three.
The wild swings near the edges are known as Runge's phenomenon: high-degree polynomials fitted to equally spaced points oscillate strongly near the interval boundaries. It is the main reason degrees above 3–4 are rarely used in practice.
Given input x, build the design matrix with columns 1, x, x2, …, xd (a Vandermonde matrix) and solve least squares exactly as for linear regression: β̂ = argmin ‖y − Vβ‖2. With multiple inputs, PolynomialFeatures also adds interaction terms such as x1x2 and x12x2; interaction_only=True keeps only products of distinct features.
The feature count grows combinatorially: p inputs at degree d give C(p + d, d) columns including the bias. That explosion is the practical limit of the method on wide data. The same expansion trick also gives logistic regression curved decision boundaries, and the polynomial kernel in SVMs computes the same feature space implicitly.
Treat d as a hyperparameter: fit each candidate degree inside cross-validation and pick the one with the lowest validation error, or the smallest degree within one standard error of the best (the "one-SE rule", which favours simpler models). A more robust alternative is to choose a generous degree and add an L2 penalty with ridge regression, tuning the penalty strength instead of the degree — this shrinks the high-order coefficients smoothly rather than all-or-nothing.
x ranges over 0–1000, x5 reaches 1015 and columns become nearly collinear. Center and scale x before expanding, or use orthogonal polynomial bases (Legendre, Chebyshev via numpy.polynomial).x.SplineTransformer) fit low-degree pieces joined smoothly at knots, which avoids Runge oscillation and behaves well at the edges. Generalized additive models (GAMs) generalize this to many features; tree ensembles find non-linearity automatically.
import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler, SplineTransformer
from sklearn.linear_model import LinearRegression, RidgeCV
from sklearn.model_selection import cross_val_score, KFold
rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 10, 80))
y = 2 + 1.6 * np.sin(x / 1.6) + 0.15 * x + rng.normal(0, 0.35, x.size)
X = x.reshape(-1, 1)
cv = KFold(5, shuffle=True, random_state=0)
def cv_rmse(model):
return -cross_val_score(model, X, y, cv=cv, scoring="neg_root_mean_squared_error").mean()
# 1) Choose the degree by cross-validation
for d in [1, 2, 3, 4, 6, 9, 12]:
m = make_pipeline(StandardScaler(), PolynomialFeatures(d, include_bias=False),
StandardScaler(), LinearRegression())
print(f"degree {d:2d} CV RMSE {cv_rmse(m):.3f}")
# 2) High degree + ridge penalty: let regularization pick the complexity
ridge = make_pipeline(StandardScaler(), PolynomialFeatures(10, include_bias=False),
StandardScaler(), RidgeCV(alphas=np.logspace(-4, 3, 30)))
print("deg 10 + ridge CV RMSE", round(cv_rmse(ridge), 3))
# 3) Cubic splines: local, well-behaved at the edges
spline = make_pipeline(SplineTransformer(n_knots=6, degree=3), RidgeCV(alphas=np.logspace(-4, 3, 30)))
print("cubic spline CV RMSE ", round(cv_rmse(spline), 3))
Scaling before the expansion keeps powers in a sane numeric range; scaling after it gives the ridge penalty comparable columns. Expect the CV error to drop to a minimum around degree 3–4 and then rise, while the ridge and spline versions stay close to the best without hand-picking the degree.
| Strength | Limitation |
|---|---|
| Simple way to add curvature to a linear model. | Feature count explodes with many inputs. |
| Keeps OLS theory, closed-form fit and fast inference. | Unstable at high degree; Runge oscillation at edges. |
| Low-degree terms are interpretable (e.g. diminishing returns). | Dangerous extrapolation outside the data range. |
Good uses: a known smooth physical or dose–response relationship, a quadratic term to capture diminishing returns, or adding a few interaction terms to an otherwise linear model. For flexible curve fitting, prefer splines; for many interacting features, prefer gradient-boosted trees.
PolynomialFeatures outside the pipeline and leaking test data into scaling.Take four points (0, 1), (1, 2), (2, 5), (3, 10). Ordinary least squares on a straight line gives slope Sxy / Sxx = 15 / 5 = 3 and intercept 4.5 − 3·1.5 = 0, so ŷ = 3x. The residuals are +1, −1, −1, +1 and the sum of squared errors is 4. The sign pattern matters more than the size: positive at both ends and negative in the middle is the classic U shape that says the model is missing curvature. Adding an x2 column and refitting gives ŷ = 1 + x2, which passes through every point.
In real data the residuals will never be exactly zero, but the diagnostic is the same: plot residuals against each input and against the fitted values. Systematic curvature suggests a power term or a spline; a funnel shape suggests heteroscedasticity (try a log transform of y); isolated large residuals deserve inspection before any polynomial is allowed to bend toward them.
Individual coefficients of correlated power terms are hard to read, but low-degree fits have useful closed-form summaries. For a quadratic y = β0 + β1x + β2x2:
β1 + 2β2x, so the effect of x depends on where you are. Report it at a few meaningful values (median, quartiles) rather than quoting β1 alone.x* = −β1 / (2β2). Check that x* lies inside the data range before claiming an optimum exists.x1x2 term means the slope of one input changes with the level of the other. Centering inputs first makes the main-effect coefficients interpretable as slopes at the mean.| Need | Polynomial | Spline / GAM | Tree ensemble |
|---|---|---|---|
| Smooth curve, 1–3 inputs | Good at low degree | Best default | Step-shaped |
| Many inputs, unknown interactions | Feature explosion | Additive only, unless interactions added | Strong |
| Closed-form summary (optimum, slope) | Yes | Numerical only | No |
| Behaviour outside data range | Diverges | Linear or constant beyond boundary knots, depending on the spline | Constant |
PolynomialFeatures and SplineTransformer in the preprocessing module.PolynomialFeatures, SplineTransformer.