What you’ll learn

• why adding $x^2$ to the feature matrix lets a linear model fit a curve
• how to solve a quadratic fit by hand on five points
• how PolynomialFeaturesPolynomialFeatures explodes combinatorially, and how fast
• how to choose the degree with cross-validation and the one-standard-error rule
• how to read a learning curve to tell underfitting from overfitting
• why polynomial pipelines must scale their features

Intuition

A straight line through curved data is wrong in a specific, visible way: it misses low, then high, then low again. The residual plot shows a smile or a frown instead of a shapeless band.

The fix is disarmingly cheap. Do not change the algorithm — change the columns. Hand the model $x^2$ alongside $x$, and least squares will find the best parabola exactly the way it found the best line, because the model is still a weighted sum of whatever columns you supply. Linear regression is linear in the parameters, and never had an opinion about the features.

  flowchart LR
  X["Raw feature x"] --> PF["PolynomialFeatures(degree=2)"]
  PF --> C["Columns: x, x²"]
  C --> LR["The same LinearRegression"]
  LR --> Y["Curved prediction ŷ"]

The math

The model

A degree-$d$ polynomial in one feature:

Rename $z_1 = x,\; z_2 = x^2,\; \ldots,\; z_d = x^d$ and it reads $\hat{y} = \theta_0 + \theta_1 z_1 + \cdots + \theta_d z_d$ — an ordinary multiple linear regression. The Normal Equation applies unchanged:

That matrix has a name — the Vandermonde matrix — and a reputation: it becomes badly conditioned fast as $d$ grows, which is one reason high-degree fits misbehave numerically as well as statistically.

How many terms?

With $n$ input features and degree $d$, PolynomialFeaturesPolynomialFeatures generates every product of total degree at most $d$:

Two features at degree 3 already give you interaction terms like $x_0^2 x_1$ — which is genuinely useful, and also how a tidy 10-column dataset becomes 286 columns without anyone deciding to do that.

Worked example by hand

Five points, $x$ chosen symmetric about zero so that $\sum x_i = 0$ and $\sum x_i^3 = 0$. That kills half the terms in $\mathbf{X}^\top\mathbf{X}$ and makes the algebra tractable.

Step 1 — assemble the normal equations. With $\sum x = 0$, $\sum x^2 = 10$, $\sum x^3 = 0$ and $\sum x^4 = 34$:

Step 2 — the middle row is already solved. $10\theta_1 = -2$, so $\theta_1 = -0.2$.

Step 3 — solve the remaining 2×2 system.

Substitute into $10\theta_0 + 34\theta_2 = 78$:

The fitted curve is $\hat{y} = 1 - 0.2x + 2x^2$.

Step 4 — check.

Nothing here is new machinery. It is the Normal Equation from the previous page, applied to a design matrix that happens to contain a squared column.

See it move

Three capacities cycling over the same curved dataset — watch the degree-12 curve chase individual points instead of the trend.

figureThe same 40 points, three model capacitiesmatplotlib

From scratch

Build the Vandermonde matrix explicitly and reuse the Normal Equation:

import numpy as np
 
 
def vandermonde(x, degree):
    """Columns [1, x, x^2, ..., x^degree]."""
    x = np.asarray(x, dtype=float)
    return np.vstack([x**k for k in range(degree + 1)]).T
 
 
x = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
y = np.array([9.0, 4.0, 1.0, 2.0, 9.0])
 
X = vandermonde(x, 2)
print(X)
# [[ 1. -2.  4.]
#  [ 1. -1.  1.]
#  [ 1.  0.  0.]
#  [ 1.  1.  1.]
#  [ 1.  2.  4.]]
 
theta = np.linalg.solve(X.T @ X, X.T @ y)
print(theta.round(4))              # [ 1.  -0.2  2. ]
 
resid = y - X @ theta
print(resid.round(4))              # [-0.4  0.8  0.  -0.8  0.4]
print(round((resid ** 2).sum(), 4))  # 1.6

import numpy as np
 
 
def vandermonde(x, degree):
    """Columns [1, x, x^2, ..., x^degree]."""
    x = np.asarray(x, dtype=float)
    return np.vstack([x**k for k in range(degree + 1)]).T
 
 
x = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
y = np.array([9.0, 4.0, 1.0, 2.0, 9.0])
 
X = vandermonde(x, 2)
print(X)
# [[ 1. -2.  4.]
#  [ 1. -1.  1.]
#  [ 1.  0.  0.]
#  [ 1.  1.  1.]
#  [ 1.  2.  4.]]
 
theta = np.linalg.solve(X.T @ X, X.T @ y)
print(theta.round(4))              # [ 1.  -0.2  2. ]
 
resid = y - X @ theta
print(resid.round(4))              # [-0.4  0.8  0.  -0.8  0.4]
print(round((resid ** 2).sum(), 4))  # 1.6

Same numbers as the hand calculation, including the residuals.

With scikit-learn

Never fit PolynomialFeaturesPolynomialFeatures outside a pipeline — the expansion has to be re-applied identically at prediction time, and the scaler after it has to be fitted on training data only:

import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
 
rng = np.random.default_rng(7)
x = np.sort(rng.uniform(-3, 3, 120)).reshape(-1, 1)
y = 0.5 * x.ravel() ** 2 + x.ravel() + 2 + rng.normal(0, 1.1, 120)
 
model = make_pipeline(
    PolynomialFeatures(degree=2, include_bias=False),
    StandardScaler(),        # essential: x^5 dwarfs x, and scale wrecks conditioning
    LinearRegression(),
)
model.fit(x, y)
 
print(model.score(x, y).round(4))                   # 0.8196
print(model[0].get_feature_names_out())             # ['x0' 'x0^2']
print(model[-1].coef_.round(3))                     # [1.78  1.248]

import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
 
rng = np.random.default_rng(7)
x = np.sort(rng.uniform(-3, 3, 120)).reshape(-1, 1)
y = 0.5 * x.ravel() ** 2 + x.ravel() + 2 + rng.normal(0, 1.1, 120)
 
model = make_pipeline(
    PolynomialFeatures(degree=2, include_bias=False),
    StandardScaler(),        # essential: x^5 dwarfs x, and scale wrecks conditioning
    LinearRegression(),
)
model.fit(x, y)
 
print(model.score(x, y).round(4))                   # 0.8196
print(model[0].get_feature_names_out())             # ['x0' 'x0^2']
print(model[-1].coef_.round(3))                     # [1.78  1.248]

include_bias=Falseinclude_bias=False matters: LinearRegressionLinearRegression already fits an intercept, so leaving the constant column in gives you two intercepts and a rank-deficient design.

Choosing the degree

Training error falls forever as you add degrees. It has to — each extra term is a strictly larger family of curves. The only honest signal comes from data the model has not seen:

figureTraining error versus cross-validated error by degreematplotlib

Reading the plot

1. The cliff at degree 2 is the real signal. CV error drops from 2.79 to 1.06 and never improves meaningfully again. The data was generated by a quadratic, and the curve says so.
2. Degrees 2 through 10 are a tie. The shaded band is one standard error across folds; every point in that range sits inside it. Picking the exact minimum here would be reading noise.
3. The one-standard-error rule. Find the best CV score, add one standard error, and take the simplest model still under that threshold. Here it selects degree 2 — the right answer, and the one that generalises.
4. The gap after degree 10 is overfitting made visible. Training error keeps falling while CV error climbs.

Learning curves

The degree curve varies capacity at fixed data. A learning curve does the opposite: fixed capacity, growing data. The two failure modes look completely different.

figureLearning curves for an underfitting and an overfitting modelmatplotlib

Pitfalls

Compare

If the curvature is monotone, try transforming $y$ or $x$ with a log before reaching for a polynomial. One column, no degree to tune, and it extrapolates sensibly.

🧪 Try It Yourself

Exercise 1 – Expand features by hand

Exercise 2 – Solve the quadratic fit

Exercise 3 – Build a polynomial pipeline

Exercise 4 – Watch overfitting appear

Exercise 5 – Extrapolate and watch it explode

Recap

• Adding $x^2, x^3, \ldots$ as columns lets ordinary least squares fit curves — the model stays linear in its parameters.
• The hand-worked quadratic on five points gives $\hat{y} = 1 - 0.2x + 2x^2$ with $R^2 = 0.972$.
• PolynomialFeaturesPolynomialFeatures produces $\binom{n+d}{d}$ terms; 10 features at degree 3 is 286 columns.
• Choose the degree with cross-validation, then apply the one-standard-error rule to prefer the simplest model in the tie.
• Learning curves distinguish the two failure modes: converged-and-high is bias, wide-gap is variance.
• Scale after expanding, and never extrapolate.

Exercise 6 – Apply the one-standard-error rule

Next

Continue to Cost Functions - Mean Squared Error (MSE) — a closer look at the loss every fit on this page was silently minimising, and what changes when you pick a different one.