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.

The Plan — From Failing Line to Curved Fit

We'll start with a dataset where a straight line is structurally wrong (residuals all positive and growing), create polynomial features, fit degree 2, then compare degrees and use cross-validation to choose the right complexity.

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}")

Degree 1 — Train MSE: 7.6452

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.

✓ Step 1 complete. Degree-1 underfits: MSE = 7.64 with all-positive, growing residuals — clear structural mismatch.

Step 2 — 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])

Feature names: ['x', 'x^2']
First 3 rows of X_poly:
[[ 1.  1.]
 [ 2.  4.]
 [ 3.  9.]]

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

✓ Step 2 complete. The model is linear in parameters — $x^2$ is just a new column in the design matrix. OLS still works.

Step 3 — 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}")

w₀ = 0.5071
w₁ = 0.3929
w₂ = 0.3821

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}")

Degree 2 — Train MSE: 0.2314

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

✓ Step 3 complete. Degree-2 fit: MSE = 0.23 (33× improvement). $w_0=0.51$, $w_1=0.39$, $w_2=0.38$.

Step 4 — 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 1: Train MSE = 7.645200
Degree 2: Train MSE = 0.231400
Degree 3: Train MSE = 0.209800
Degree 7: Train MSE = 0.000000

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.

✓ Step 4 complete. Degree 2 is the sweet spot: low train + low test. Degree 7 overfits (zero train error, wild oscillations).

Step 5 — 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 1: CV MSE = 9.1 ± 4.2
Degree 2: CV MSE = 1.1 ± 0.9
Degree 3: CV MSE = 1.8 ± 1.5
Degree 4: CV MSE = 4.2 ± 6.1
Degree 5: CV MSE = 22.1 ± 31.7

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. cv=4 uses 4-fold CV (small for 8 samples, but demonstrates the method).

✓ Step 5 complete. 4-fold CV confirms degree 2: CV MSE = 1.1 ± 0.9, lowest and most stable.

Step 6 — 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])

Feature names: ['sqft', 'beds', 'sqft beds']
First row: [ 650.    2.  1300.]

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.

✓ Step 6 complete. interaction_only=True adds cross-terms without squared terms — useful when the effect of one feature depends on another.

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 why any feature transformation — square root, log, reciprocal — can be passed to linear regression. The model doesn't know or care where the features came from: $\hat{y}=w_0+w_{1}x+w_2\log (x)$ is as valid as polynomials. This is the foundation of basis function expansion.

A pragmatic frustration: choosing degree by CV on 8 samples is unreliable — 4-fold CV leaves only 6 training points per fold. The ±31.7 standard deviation for degree 5 makes the "optimal" degree ambiguous. On real datasets with hundreds of samples, CV degree selection stabilizes. For small datasets, domain physics ("the heat equation predicts $T\propto A$") should constrain the degree before CV confirms it.

Another blind spot: the feature explosion with $p$ features and degree $d$ is $O(p^d)$. For $p=10$ and $d=5$, you get 3,003 features — Ridge or Lasso becomes necessary. And extrapolation with high-degree polynomials is genuinely dangerous: a degree-7 fit that's perfect on $[1,8]$ predictably diverges at $x=9$ or $x=0$, often in the wrong direction.

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?