Polynomial Features and Interaction Terms in Python

Contents

Figure 1: The equation doesn't change — you just hand the model more terms to fit it with

Extra columns make straight-line models bend, condition, and occasionally explode

A linear model draws straight lines, which is a problem when the real world bends. Sales don't rise forever with ad spend, dosage effects depend on body weight, and a discount means something different on a holiday than on a Tuesday.

Polynomial features let a linear model fit curves, and interaction terms let one feature's effect depend on another. The trick is that you don't change the model, you change the inputs: add columns like x² and x₁·x₂, and ordinary linear or logistic regression can capture non-linear patterns while staying fast and interpretable. This guide covers how to generate these features in scikit-learn, how to keep them from overfitting or exploding, and when splines, trees, or other tools are the better choice.

What These Features Are

Polynomial features are powers of a single feature: x, x², x³, and so on. Interaction terms are products of different features: x₁·x₂. Combining both for two inputs a and b, the full degree-2 expansion is:

[1, a, b, a², ab, b²]

The model is still "linear" because it's linear in the coefficients: prediction = w₀ + w₁a + w₂b + w₃a² + w₄ab + w₅b². The fitting machinery, regularization, and interpretation tools all still apply — you've just given the model a richer set of building blocks.

When They Help

• Curved relationships: diminishing returns, U-shapes, and accelerating effects.
• Effects that depend on context: a feature's slope changes with the level of another feature or a category.
• Interpretability: unlike a black-box ensemble, you can read an interaction coefficient directly.
• Small to medium tabular data, where a regularized linear model with a few well-chosen terms can compete with heavier models.
• Inference: statistical models need explicit interaction terms, since they won't discover them for you.

Tree-based models and neural networks already discover non-linearities and interactions on their own, so adding polynomial features to them is mostly wasted effort.

Basic Usage

import numpy as np
from sklearn.preprocessing import PolynomialFeatures

X = np.array([[2, 3]])

poly = PolynomialFeatures(degree=2, include_bias=True).fit(X)
print(poly.get_feature_names_out(["a", "b"]))   # ['1' 'a' 'b' 'a^2' 'a b' 'b^2']
print(poly.transform(X))                         # [[1. 2. 3. 4. 6. 9.]]

inter = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False).fit(X)
print(inter.get_feature_names_out(["a", "b"]))   # ['a' 'b' 'a b']

The key parameters, per scikit-learn's documentation:

• degree: the maximum degree. It can also be a tuple (min_degree, max_degree) if you want only certain degrees.
• interaction_only: when True, you get only products of distinct features, so squares and higher powers of the same feature are excluded.
• include_bias: adds a constant column of ones. Set it to False when your model already fits an intercept (which most do) to avoid a redundant column.
• get_feature_names_out: gives readable names, which matters when you inspect coefficients later.

Fitting a Curve, and Watching It Overfit

Here's a noisy cubic relationship, fitted with increasing polynomial degrees:

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.default_rng(0)
x = rng.uniform(-3, 3, 120)
y = 0.4 * x**3 - x**2 + rng.normal(0, 2, x.size)
X = x.reshape(-1, 1)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.4, random_state=0)

for degree in (1, 2, 3, 5, 10, 15):
    model = make_pipeline(PolynomialFeatures(degree, include_bias=False),
                          StandardScaler(), LinearRegression())
    model.fit(Xtr, ytr)
    tr = mean_squared_error(ytr, model.predict(Xtr)) ** 0.5
    te = mean_squared_error(yte, model.predict(Xte)) ** 0.5
    print(f"degree {degree:>2}   train RMSE {tr:5.2f}   test RMSE {te:5.2f}")

Expect the textbook pattern:

Degree	Train RMSE	Test RMSE	Reading
1	high	high	underfits — a line can't follow a cubic
2–3	low	low	the sweet spot; degree 3 matches the true shape
5–10	lower	climbing	variance starts to win
15	near zero	worst	a degree-15 polynomial happily threads through the noise

Note the StandardScaler after expansion: x¹⁵ has an enormous range compared with x, which makes fitting numerically unstable without scaling.

Tame It With Regularization and Cross-Validation

Rather than guessing the degree, tune it together with a ridge penalty:

from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV

pipe = make_pipeline(PolynomialFeatures(include_bias=False), StandardScaler(), Ridge())
grid = GridSearchCV(
    pipe,
    {"polynomialfeatures__degree": [1, 2, 3, 4, 5, 8, 10],
     "ridge__alpha": np.logspace(-3, 3, 7)},
    cv=5, scoring="neg_root_mean_squared_error",
)
grid.fit(Xtr, ytr)

test_rmse = mean_squared_error(yte, grid.predict(Xte)) ** 0.5
print(grid.best_params_, "| CV RMSE:", round(-grid.best_score_, 2),
      "| test RMSE:", round(test_rmse, 2))

The ridge penalty shrinks the wild high-degree coefficients, so a moderately high degree with strong regularization often performs about as well as the right degree without it — and standardizing first matters, because penalties treat all coefficients equally, exactly as the ridge vs. lasso guide explains.

Interaction Terms in Practice

Imagine a marketing dataset where advertising works better during holidays. The effect of ad spend isn't constant — it depends on the holiday flag, and only an interaction term can capture that:

import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(1)
n = 500
df = pd.DataFrame({"ad_spend": rng.uniform(0, 100, n),
                   "holiday": rng.integers(0, 2, n)})
df["sales"] = (50 + 0.8 * df.ad_spend + 20 * df.holiday
               + 0.6 * df.ad_spend * df.holiday + rng.normal(0, 10, n))

X, y = df[["ad_spend", "holiday"]], df["sales"]

models = {
    "main effects only": LinearRegression(),
    "with interaction": make_pipeline(
        PolynomialFeatures(degree=2, interaction_only=True, include_bias=False),
        LinearRegression()),
}
for name, m in models.items():
    r2 = cross_val_score(m, X, y, cv=5, scoring="r2").mean()
    print(f"{name:<18} CV R² = {r2:.3f}")

fit = models["with interaction"].fit(X, y)
names = fit[0].get_feature_names_out(["ad_spend", "holiday"])
print(dict(zip(names, fit[-1].coef_.round(2))))

The main-effects model assumes ads work equally well on holidays and normal days, so it leaves variance on the table. The interaction model recovers the structure we planted: the ad_spend holiday coefficient is the extra sales per ad dollar on holidays. That's the interpretive payoff — you can say "each ad dollar yields about 0.8 units on normal days and about 1.4 on holidays."

Interpreting Interactions Correctly

With an interaction in the model, the coefficient on each main effect no longer means "the effect of this feature." It means the effect when the other feature is zero — often meaningless, such as zero temperature or zero age. Centering numeric features (subtracting the mean) before creating the interaction fixes this: the main-effect coefficients then describe the effect at the average level of the other feature, and the correlation between main effects and the interaction term drops. Keep the centering inside your pipeline so the mean is learned from training data only.

Interactions With Categories

Interactions between a numeric feature and a categorical one let each category have its own slope. One-hot encode the category and multiply by the numeric feature — in a pipeline, OneHotEncoder followed by PolynomialFeatures(interaction_only=True) creates all pairwise products among encoded columns. It's useful but expands quickly, so apply it to targeted columns with a ColumnTransformer rather than exploding the whole feature space.

Statsmodels for Inference

If you need standard errors, confidence intervals, and p-values rather than just predictions, statsmodels' formula interface writes interactions concisely:

import statsmodels.formula.api as smf

res = smf.ols("sales ~ ad_spend * holiday", data=df).fit()
print(res.summary().tables[1])
# "a * b" expands to a + b + a:b ; use "a:b" for the interaction alone,
# and I(x**2) for a squared term, e.g. "sales ~ ad_spend + I(ad_spend**2)"

Remember that standard inference assumes you chose the terms before looking at the data — hunting through many interactions invalidates the p-values.

The Feature Explosion Problem

The number of generated features grows quickly. For n features and degree d, the count (excluding the bias column) is C(n + d, d) − 1:

from math import comb

for n_features in (5, 20, 100):
    for d in (2, 3):
        print(f"{n_features:>3} features, degree {d}: {comb(n_features + d, d) - 1:>8,} columns")

Input features	Degree 2	Degree 3
5	20	55
20	230	1,770
100	5,150	176,850

A modest dataset can't support hundreds of thousands of coefficients without severe overfitting, and memory becomes an issue too. If you pass sparse data, scikit-learn's implementation accepts it and returns CSR output (CSC required for degrees 4 and above), but sparsity doesn't remove the combinatorial growth.

Controlling Complexity

• Choose terms deliberately. Domain knowledge beats brute force — add the interactions you have a reason to expect, such as price × promotion or age × dosage, instead of every pair.
• Use interaction_only=True when you want products between features without squares and cubes.
• Regularize. Ridge shrinks the many correlated terms, and lasso or elastic net can zero out useless ones — automatic selection among the expanded features, as in the ridge vs. lasso comparison.
• Respect hierarchy. The heredity principle says that if an interaction x₁x₂ is in the model, keep the main effects x₁ and x₂ — an interaction without its parents is hard to interpret and rarely justified.
• Find candidates before adding them. Partial dependence and ICE plots reveal where effects vary by another feature, and Friedman's H-statistic quantifies interaction strength. Tree-based models can also hint at which pairs matter.
• Mind multicollinearity. x and x² are strongly correlated, as are main effects and their interaction. Centering, scaling, and regularization help, and individual coefficient values can be unstable even when predictions are fine.

Splines: Often a Better Way to Curve

High-degree polynomials have two well-known weaknesses. They oscillate wildly near the edges of the data (Runge's phenomenon), and they extrapolate badly — a degree-10 polynomial shoots toward infinity just beyond the training range. Splines fix both by fitting low-degree polynomial pieces joined smoothly at knots. SplineTransformer generates B-spline basis functions with parameters n_knots, degree, knots, extrapolation, and sparse_output. A comparison on the cubic data above makes the difference visible:

from sklearn.preprocessing import SplineTransformer

poly10 = make_pipeline(PolynomialFeatures(10, include_bias=False),
                       StandardScaler(), Ridge(alpha=1e-3)).fit(Xtr, ytr)
spline = make_pipeline(SplineTransformer(n_knots=6, degree=3),
                       Ridge(alpha=1e-3)).fit(Xtr, ytr)

for name, m in [("polynomial deg 10", poly10), ("cubic spline", spline)]:
    test_rmse = mean_squared_error(yte, m.predict(Xte)) ** 0.5
    outside = m.predict(np.array([[4.0], [5.0], [6.0]])).round(1)
    print(f"{name:<18} test RMSE {test_rmse:.2f}   predictions at x=4,5,6: {outside}")

The training range is −3 to 3, so x = 4, 5, and 6 are extrapolations. The polynomial's predictions typically balloon, while the spline with its default constant extrapolation stays flat — a more cautious failure. Neither is "right" out there, since no model can know what happens beyond its data, but one fails more gracefully. Splines also give you local control: a bump in one region doesn't distort the fit elsewhere.

Other Alternatives

• Binning (KBinsDiscretizer) turns a numeric feature into ordinal or one-hot bins, giving a piecewise-constant fit that's easy to explain — the trade-offs are laid out in the binning and discretization guide.
• Generalized additive models (GAMs) fit smooth, non-linear functions for each feature and add them up; Explainable Boosting Machines add automatic pairwise interactions with excellent interpretability.
• Tree-based ensembles discover interactions automatically, trading explicit coefficients for flexibility.
• Kernel methods such as polynomial-kernel SVMs and kernel ridge operate in a polynomial feature space implicitly without ever materializing the columns, as the SVM guide explains.
• Neural networks learn their own feature interactions given enough data.

A Practical Workflow

1. Start with a plain regularized linear model as a baseline. Plot residuals against each feature to spot curvature; use partial dependence or domain knowledge to suspect interactions.
2. Add one thing at a time — a squared term for a curved feature, or a targeted interaction — and check cross-validated performance after each change.
3. Put expansion, centering, scaling, and the model inside a single pipeline so statistics are learned on training folds only.
4. Regularize, and compare against splines and a boosted-tree model before settling.
5. Inspect coefficients with get_feature_names_out, and confirm that predictions behave sensibly at the edges of your data.

Common Pitfalls

• Skipping scaling after expansion — high-degree columns dominate numerically.
• Choosing a high degree without cross-validation or regularization.
• Expanding every pairwise interaction on a wide dataset — then drowning in features.
• Reading main-effect coefficients as "the effect" when an interaction is present — center first.
• Including an interaction without its main effects — heredity exists for a reason.
• Trusting polynomial predictions outside the training range.
• Fitting the expansion and scaler on the full dataset before splitting — leakage again.
• Adding polynomial features to tree models or neural nets and expecting a boost.
• Hunting through many interactions and reporting uncorrected p-values.
• Forgetting include_bias=False — a redundant constant column next to the model's intercept.

Unlock AI That Actually Works

Get lifetime access to GPT-6 Astra, Claude Fable 5.1, Gemini 3.5, Grok 4.5, and more — all in one platform. Build websites, apps, videos, content, and digital products from a single command. No monthly fees. No tool-hopping.

Click here to get GPTAstra Max now — one-time payment, lifetime access.

Frequently Asked Questions

When should you use polynomial features?

When a linear model needs to bend: curved relationships (diminishing returns, U-shapes), effects that change with context, or small-to-medium tabular data where a regularized linear model with a few well-chosen terms competes with heavier models. They also keep interpretability — you can read an interaction coefficient directly — and are required for statistical inference, since statsmodels won't discover interactions for you. Don't bother for tree-based models or neural networks: they find non-linearities and interactions on their own. When the curve is smooth, splines usually beat high-degree polynomials; when the truth has real thresholds, binning fits it better.

How do you stop polynomial features from overfitting?

Three moves, in order: scale after expansion (x^15 has an enormous range compared with x, which makes fitting numerically unstable), regularize with ridge so the wild high-degree coefficients get shrunk, and tune degree and alpha together with cross-validation instead of guessing — GridSearchCV over polynomialfeatures__degree and ridge__alpha inside one pipeline. A moderately high degree with strong regularization often performs about as well as the right degree without it. Keep the expansion and scaler inside the pipeline so statistics are learned on training folds only.

What is the difference between polynomial features and interaction terms?

Polynomial features are powers of a single feature (x, x², x³); interaction terms are products of different features (x₁·x₂). For two inputs a and b, the full degree-2 expansion is [1, a, b, a², ab, b²] — the model stays "linear" because it's linear in the coefficients. PolynomialFeatures(degree=2, interaction_only=True) gives you only the products (a, b, ab) with no squares, which is the choice when you believe one feature's effect depends on another but not in curvature. Respect the heredity principle: if x₁x₂ is in the model, keep x₁ and x₂.

When are splines better than high-degree polynomials?

Almost always beyond degree 3–4. High-degree polynomials oscillate wildly near the edges of the data (Runge's phenomenon) and extrapolate explosively — a degree-10 term shoots toward infinity just beyond the training range. SplineTransformer fits low-degree pieces joined at knots, so a bump in one region doesn't distort the fit elsewhere, and its default constant extrapolation fails flat instead of ballooning: a more cautious failure, since no model knows what happens beyond its data. Use polynomials for a targeted quadratic term, splines for flexible curves.

Wrapping This Up

Polynomial features and interaction terms let linear models bend and condition, using nothing more exotic than extra columns. PolynomialFeatures handles the expansion, interaction_only and degree control it, and a pipeline with scaling and regularization keeps it honest. The costs are feature explosion, collinearity, and poor extrapolation — which is why targeted, domain-informed terms beat brute-force expansion, and why splines, GAMs, or boosted trees are often the better answer for flexible curves.

How should you start? Fit the baseline, plot your residuals, add the one interaction or squared term you have a reason to believe in, and let cross-validation say whether it earned its place. Then run the spline comparison above and see how differently the two models behave once x leaves the range you trained on.

Related Articles

• Ridge vs. Lasso: Regularization Explained
• Partial Dependence and ICE Plots for Feature Effects
• Binning and Discretization: When to Turn Numbers into Categories
• Support Vector Machines Explained: Kernels and Python