Polynomial Regression

When the relationship between input and output bends, a straight line fails — not because linear regression is wrong, but because the feature set is wrong. Polynomial regression fixes this by engineering new features from the original ones. The algorithm doesn't change. The loss function doesn't change. OLS still applies. What changes is the design matrix, and that single insight unlocks an entire class of non-linear models without leaving the linear regression framework.

Anchor dataset: a material experiment measuring temperature rise ($y$, in °C) as a function of applied area ($x$, in cm²). The physical relationship is roughly quadratic — heat scales with area, not linearly.

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

X_raw = np.array([1, 2, 3, 4, 5, 6, 7, 8]).reshape(-1, 1)
y     = np.array([1.5, 3.2, 5.9, 9.1, 12.8, 17.5, 22.4, 29.1])

Why Linear Regression Fails Here

Fit a straight line first:

model_d1 = LinearRegression()
model_d1.fit(X_raw, y)
y_pred_d1 = model_d1.predict(X_raw)
mse_d1 = mean_squared_error(y, y_pred_d1)
print(f"Degree 1 — Train MSE: {mse_d1:.4f}")

A high MSE alone doesn't prove underfit — you need to look at the residual pattern. Compute them:

All residuals are positive and growing. This is a systematic pattern — the line consistently undershoots, with the gap widening as $x$ increases. Random noise would produce positive and negative residuals with no trend. The curve is visibly bending away from the line.

Every orange dot sits above the blue line. The gap grows from left to right. This is underfit by structural mismatch — linear regression cannot model a curve.

The Key Insight: Polynomial Is Feature Engineering

The fix is not a new algorithm. The fix is a new feature.

Instead of passing $x$ to linear regression, pass $[x,x^2]$. The model becomes:

$\hat{y}=w_0+w_{1}x+w_2{x}^2$

This is still linear in the parameters $[w_0,w_1,w_2]$. OLS can still minimize SSE exactly. The non-linearity is in the input space — the curve lives in $x$ space. In the augmented feature space $[x,x^2]$, the model is a hyperplane.

Feature transformation trace — all 8 samples, degrees 2 and 3:

Each row is now a vector in a higher-dimensional space. Linear regression in that space finds the coefficients $[w_0,w_1,w_2]$ that minimize SSE — it has no idea (and doesn't need to know) that $x^2$ came from squaring $x$.

Using sklearn PolynomialFeatures

poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(X_raw)
print("Feature names:", poly.get_feature_names_out(['x']))
print("First 3 rows of X_poly:")
print(X_poly[:3])

include_bias=False drops the redundant column of ones — LinearRegression adds its own intercept. The output matches the trace table exactly.

Degree-2 Fit — Full Trace

model_d2 = make_pipeline(PolynomialFeatures(2), LinearRegression())
model_d2.fit(X_raw, y)

w = model_d2.named_steps['linearregression']
print(f"w₀ = {w.intercept_:.4f}")
print(f"w₁ = {w.coef_[0]:.4f}")
print(f"w₂ = {w.coef_[1]:.4f}")

Manual prediction verification for the first 3 samples:

The residuals are still positive but now small and no longer trend-shaped.

mse_d2 = mean_squared_error(y, model_d2.predict(X_raw))
print(f"Degree 2 — Train MSE: {mse_d2:.4f}")

From 7.64 (degree 1) to 0.23 (degree 2): a 33× improvement from adding a single engineered feature ($x^2$).

Degree Comparison

for d in [1, 2, 3, 7]:
    m = make_pipeline(PolynomialFeatures(d), LinearRegression())
    m.fit(X_raw, y)
    train_mse = mean_squared_error(y, m.predict(X_raw))
    print(f"Degree {d}: Train MSE = {train_mse:.6f}")

Degree 7 achieves zero training error because 8 parameters can exactly fit 8 points. On unseen $x$ values, it oscillates wildly between the training points.

The degree-7 curve passes through every blue dot — perfect training fit — but oscillates between them. Any new observation falling between training points would get a wildly inaccurate prediction.

Hyperparameter: Degree Selection via Cross-Validation

Training MSE always decreases as degree increases. Cross-validation separates training error from generalization error:

from sklearn.model_selection import cross_val_score

for d in [1, 2, 3, 4, 5]:
    m = make_pipeline(PolynomialFeatures(d), LinearRegression())
    cv_mse = -cross_val_score(m, X_raw, y, cv=4, scoring='neg_mean_squared_error')
    print(f"Degree {d}: CV MSE = {cv_mse.mean():.1f} ± {cv_mse.std():.1f}")

Degree 2: lowest CV MSE and lowest variance. From degree 3 onward, CV MSE rises and its standard deviation explodes — the model fits training folds well but fails on held-out folds. Degree 2 is confirmed optimal.

Interaction Terms

For $p\ge 2$ features, PolynomialFeatures generates cross-terms by default. For degree 2 on $[x_1,x_2]$:

$[x_1,x_2]\to [x_1,x_2,x_1^2,x_1{x}_2,x_2^2]$

The term $x_1{x}_2$ captures: "the effect of $x_1$ on $y$ depends on the level of $x_2$." Without it, the model assumes additive effects — each feature's contribution is independent of the others.

Applied to our house dataset ($x_1$ = sqft, $x_2$ = bedrooms):

from sklearn.datasets import fetch_california_housing
X_house = np.array([[650, 2], [850, 2], [1100, 3],
                    [1400, 3], [1600, 4], [1900, 4]])
y_house = np.array([180, 220, 280, 340, 370, 430])

# Only the cross-term, no squared terms
poly_interact = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_interact = poly_interact.fit_transform(X_house)
print("Feature names:", poly_interact.get_feature_names_out(['sqft', 'beds']))
print("First row:", X_interact[0])

The cross-term sqft * beds = 1300 for the first house. A large-bedroom house at high sqft might fetch a disproportionate premium — the interaction term lets the model capture that non-additive effect.

When Polynomial Regression Is Not the Answer

• Degree > 3: Numerical instability at high degree. Polynomials oscillate (Runge's phenomenon). Use tree models or neural networks instead.
• High dimension ($p$ features) + high degree: Feature count explodes as $O(p^d)$. For $p=10$, $d=3$: 286 features. Use regularized polynomial (Ridge) or kernel methods (RBF kernel SVM) instead.
• Extrapolation: A degree-7 polynomial fitted to $x\in [1,8]$ will diverge catastrophically for $x=9$ or $x=0$. Tree models don't extrapolate either, but neural networks can generalize better with enough data.

Polynomial vs Linear Comparison

Related Concepts and Honest Limitations

The "linear in parameters" argument explains exactly why polynomial regression doesn't require any new machinery. Any transformation of $x$ — square root, log, reciprocal — can be treated as a new feature and passed to linear regression. The model is not restricted to polynomial terms: $\hat{y}=w_0+w_{1}x+w_2\log (x)$ is equally valid.

The honest limitation: choosing degree by CV on 8 samples is unreliable — 4-fold CV leaves only 6 training points per fold, and the high variance in the CV output (degree 5 shows ±31.7) reflects this. On real datasets with hundreds of samples, CV degree selection is much more stable. For small datasets, prior domain knowledge ("the physics suggests a quadratic relationship") should constrain the degree before CV confirms it.

Test Your Understanding

1. Degree 7 achieves Train MSE = 0 on 8 samples. How many parameters does a degree-7 polynomial in one feature have, and why does this guarantee zero training error?

2. The degree-2 model has $w_2=0.38$. If you standardize $x$ before applying PolynomialFeatures, will $w_2$ change? Will the predictions on the original scale change?

3. Cross-validation shows degree-2 CV MSE = 1.1 ± 0.9 and degree-3 CV MSE = 1.8 ± 1.5. The mean for degree 3 is higher, but the confidence intervals overlap. What additional analysis would you run before committing to degree 2?

4. For 2 features with degree=2, PolynomialFeatures produces 5 features $[x_1,x_2,x_1^2,x_1{x}_2,x_2^2]$. For 3 features with degree=2, how many features does it produce? Write out the formula.

5. A degree-5 polynomial is fit to a dataset from $x\in [0,10]$. A new observation arrives at $x=15$. Why is the polynomial's prediction at $x=15$ likely unreliable, even if the model fits the training data well?