A curved relationship, fitted by a linear model
Fuel consumption against speed is not a straight line. Efficiency improves up to about 70 km/h and then worsens as drag takes over, so the curve is roughly U-shaped. A straight line through it is wrong at both ends and wrong in the middle.
Polynomial regression handles this without leaving linear regression behind. You add speed² as a second column and fit the same least-squares model on two features. The relationship with speed becomes a curve; the model is still linear in its coefficients, which is why every tool from the previous two chapters still applies.
The naming confuses people and is worth settling now. “Linear” in linear regression refers to the coefficients, not the shape of the fitted line. A model with x, x² and x³ columns is still a linear model. This is feature construction of the kind covered in creating new features, applied mechanically.
How PolynomialFeatures builds the columns
import numpy as np
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures
demo = pd.DataFrame({"a": [2.0, 3.0], "b": [5.0, 7.0], "c": [1.0, 4.0]})
pf = PolynomialFeatures(degree=2, include_bias=False)
out = pf.fit_transform(demo)
print(list(pf.get_feature_names_out()))
# ['a', 'b', 'c', 'a^2', 'a b', 'a c', 'b^2', 'b c', 'c^2']
print(out.shape) # (2, 9)
Three columns became nine: the originals, each squared, and every pairwise product. Those cross terms are interactions, and they are often the reason this transformer helps even when you did not want curves.
Two arguments change what you get. include_bias=False suppresses the constant column, which you want because the estimator fits its own intercept. interaction_only=True keeps the products but drops the powers:
print(list(PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
.fit(demo).get_feature_names_out()))
# ['a', 'b', 'c', 'a b', 'a c', 'b c']
The column count grows faster than people expect.
for d in [2, 3, 4]:
n = PolynomialFeatures(degree=d, include_bias=False).fit(
pd.DataFrame(np.zeros((2, 10)))).n_output_features_
print(f"degree {d} on 10 columns -> {n} features")
# degree 2 on 10 columns -> 65 features
# degree 3 on 10 columns -> 285 features
# degree 4 on 10 columns -> 1000 features
A thousand features from ten, at degree 4. Unless you have far more rows than that, the model has enough freedom to fit noise exactly. In practice I use degree 2 on a handful of chosen columns rather than degree 3 on everything.
Choosing the degree
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
rng = np.random.default_rng(8)
speed = np.sort(rng.uniform(20, 120, 400))
fuel = 14 - 0.22 * speed + 0.0016 * speed ** 2 + rng.normal(0, 0.45, 400)
S = speed.reshape(-1, 1)
S_tr, S_te, f_tr, f_te = train_test_split(S, fuel, test_size=0.3, random_state=0)
for d in [1, 2, 3, 9]:
pipe = make_pipeline(PolynomialFeatures(degree=d, include_bias=False),
StandardScaler(), LinearRegression()).fit(S_tr, f_tr)
print(f"degree {d}: train {mean_squared_error(f_tr, pipe.predict(S_tr)):.4f}"
f" test {mean_squared_error(f_te, pipe.predict(S_te)):.4f}")
# degree 1: train 1.5452 test 1.4424
# degree 2: train 0.1954 test 0.2349
# degree 3: train 0.1952 test 0.2381
# degree 9: train 0.1889 test 0.2495
Degree 1 is badly wrong, at roughly six times the error of the alternatives. Degree 2 is best on held-out data. Degrees 3 and 9 keep improving on training data while getting worse on test data, which is the pattern described in the main challenges of machine learning.
The underlying relationship here really is quadratic, so degree 2 winning is not a coincidence. Two practical points follow. Start at degree 2 and only go higher with evidence from held-out data. And treat the degree as something to select by validation rather than by eye, since the training error will happily recommend degree 9.
Note the StandardScaler between the two steps. It is not optional here, for a reason the next section makes concrete.
Where polynomial regression breaks down
Two failure modes, and the first is severe enough to rule the method out of some problems entirely.
for d in [2, 9]:
pipe = make_pipeline(PolynomialFeatures(degree=d, include_bias=False),
StandardScaler(), LinearRegression()).fit(S_tr, f_tr)
print(f"degree {d}: at 140 km/h {float(pipe.predict([[140.0]])[0]):8.2f}"
f" at 160 km/h {float(pipe.predict([[160.0]])[0]):10.2f}")
# degree 2: at 140 km/h 14.75 at 160 km/h 20.13
# degree 9: at 140 km/h 47.74 at 160 km/h 883.71
print(round(14 - 0.22 * 140 + 0.0016 * 140 ** 2, 2)) # 14.56 (the true value)
Training data stopped at 120 km/h. Twenty kilometres beyond it, the degree 2 model predicts 14.75 litres against a true 14.56, and the degree 9 model predicts 47.74. At 160 it predicts 884 litres per 100 kilometres.
High-degree polynomials curve violently outside the range they were fitted on, and nothing in the model warns you. If your production inputs can exceed the training range, either cap the degree at 2, clip inputs to the observed range, or use a method such as splines or a tree ensemble that cannot run away.
The second failure is numerical.
raw = PolynomialFeatures(degree=9, include_bias=False).fit_transform(S)
print([f"{raw[:, i].max():.2e}" for i in [0, 4, 8]])
# ['1.20e+02', '2.48e+10', '5.11e+18']
The first column tops out at 120 and the ninth at five quintillion. Columns on such different scales make the least-squares solve ill-conditioned and any penalty term meaningless, which is why feature scaling goes between PolynomialFeatures and the estimator rather than before it. Order matters: scale the expanded columns, not the original one.
Degree 2 with scaling, validated on held-out data, covers most real curved relationships. Beyond that the returns fall away quickly and the risks do not. When a relationship genuinely needs more flexibility than a quadratic, the better answer is usually a different model family rather than a higher degree: splines fit local curvature without the global blow-up, and gradient boosting handles arbitrary shapes without any feature construction at all. Full options are in the PolynomialFeatures documentation.
Frequently Asked Questions
Is polynomial regression linear or non-linear?
It is a linear model fitted to non-linear features. The relationship between the input and the prediction curves, but the model remains linear in its coefficients, which is why ordinary least squares solves it and every diagnostic for linear regression still applies.
What degree should you use for polynomial regression?
Start at 2 and raise it only if held-out error improves. Training error always falls as degree rises, so it cannot guide the choice. Degree 2 won on the example here, and degree 9 fitted the training data better while predicting 884 litres twenty percent outside the training range.
Why do you need to scale features for polynomial regression?
Raising a column to the ninth power produces values around five quintillion while the original column tops out near 120. Those scales make the solve numerically unstable and any regularisation penalty meaningless. Scale after the polynomial expansion, not before it.
Key Takeaways
- Add
PolynomialFeatures(degree=2, include_bias=False)before your estimator when a relationship curves, since the model stays linear in its coefficients. - Place the scaler after the expansion and before the estimator, because the high-power columns reach magnitudes the solver cannot handle.
- Select the degree on held-out error rather than training error, which keeps improving with every degree you add.
- Check your production input range before deploying anything above degree 2, as a degree 9 fit predicted 884 where the truth was near 15.
- Count the output features before raising the degree on a wide dataset, since ten columns become a thousand at degree 4.