Polynomial Regression

Linear regression draws a straight line through data. When the relationship between the input and output curves — think the trajectory of a ball, the growth of a bacterial colony, or a dose-response curve — a straight line is a poor fit no matter how many features you add. Polynomial regression solves this by adding powered versions of the existing features (x², x³, …) so the model can bend to match the data's curvature while still using the ordinary least-squares machinery under the hood.

This page covers:

• How polynomial regression extends linear regression mathematically
• When to use polynomial regression and what degree to pick
• A complete scikit-learn pipeline: feature transformation, training, evaluation
• The overfitting trap and how to detect it with train/test curves
• Comparing model quality with RMSE and R²

How Polynomial Regression Works

The Equation

Linear regression fits:

y = β₀ + β₁x

Polynomial regression of degree n fits:

y = β₀ + β₁x + β₂x² + β₃x³ + … + βₙxⁿ

The key insight is that x², x³, and so on are treated as additional features. The model is still linear in its coefficients (β values); only the features are non-linear. That means ordinary least squares still works — you just pre-process the input matrix first.

Degrees and What They Mean

Choosing the right degree is the central skill. Too low and the model underfits (misses real curvature). Too high and the model overfits (memorizes noise and fails on new data).

When to Use Polynomial Regression

Use polynomial regression when:

• A scatter plot shows a clear curve in the data that a straight line cannot capture
• You want a fast, interpretable alternative to tree-based models for moderate curvature
• You already have a working linear regression baseline and it underfits

Prefer other approaches when:

• The relationship is highly complex or has many features → try decision trees or gradient boosting
• You need to predict far outside the training range (extrapolation) → high-degree polynomials diverge wildly outside the data range
• You have many input features → each feature gets n polynomial terms, so the feature matrix grows fast

Building a Polynomial Regression Model

Step 1: Import Libraries

import numpy as np
import matplotlib
matplotlib.use('Agg')  # non-interactive backend for scripts
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

scikit-learn's PolynomialFeatures transformer handles the x → [1, x, x², …, xⁿ] expansion automatically. Wrapping it with LinearRegression in a Pipeline keeps the code clean and prevents data leakage during cross-validation.

Step 2: Create Sample Data

rng = np.random.default_rng(42)
X = rng.uniform(-3, 3, 80).reshape(-1, 1)   # 80 points from -3 to 3
y = 0.5 * X.ravel()**2 - X.ravel() + 2 + rng.normal(0, 0.5, 80)

# True relationship: y ≈ 0.5x² - x + 2  (a parabola with noise)

The underlying relationship is quadratic, so a degree-2 polynomial should recover it well. This is the kind of situation where linear regression systematically underfits.

Step 3: Split into Training and Test Sets

Always split before fitting so you have held-out data to evaluate on. See Train/Test Split for a full explanation.

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

print(f"Training samples: {len(X_train)}")  # 64
print(f"Test samples:     {len(X_test)}")   # 16

Step 4: Build and Fit the Pipeline

degree = 2
model = make_pipeline(
    PolynomialFeatures(degree=degree, include_bias=False),
    LinearRegression()
)
model.fit(X_train, y_train)

make_pipeline chains the two steps: PolynomialFeatures transforms the input, then LinearRegression fits on the expanded features. You call .fit(), .predict(), and .score() on the pipeline exactly as you would on a plain estimator.

Step 5: Evaluate the Model

y_pred_train = model.predict(X_train)
y_pred_test  = model.predict(X_test)

rmse_train = np.sqrt(mean_squared_error(y_train, y_pred_train))
rmse_test  = np.sqrt(mean_squared_error(y_test,  y_pred_test))
r2_train   = r2_score(y_train, y_pred_train)
r2_test    = r2_score(y_test,  y_pred_test)

print(f"Train RMSE: {rmse_train:.4f}   Test RMSE: {rmse_test:.4f}")
print(f"Train R²:   {r2_train:.4f}   Test R²:   {r2_test:.4f}")

Expected output:

Train RMSE: 0.4605   Test RMSE: 0.4538
Train R²:   0.9480   Test R²:   0.9514

An R² near 0.95 on held-out data means the model explains about 95% of the variance — an excellent fit for data this noisy. The train and test scores are close, which indicates the model generalizes rather than overfitting.

Interpreting the metrics:

• RMSE (Root Mean Squared Error) — the average error in the same units as y. Lower is better.
• R² (coefficient of determination) — fraction of variance explained. Values closer to 1.0 are better; 0 means the model is no better than predicting the mean.

Step 6: Visualize the Fit

X_line = np.linspace(-3.5, 3.5, 200).reshape(-1, 1)
y_line = model.predict(X_line)

plt.figure(figsize=(7, 5))
plt.scatter(X_train, y_train, alpha=0.6, label='Train', color='steelblue', s=25)
plt.scatter(X_test,  y_test,  alpha=0.8, label='Test',  color='orange',    s=40, zorder=5)
plt.plot(X_line, y_line, color='red', linewidth=2, label=f'Degree-{degree} fit')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Polynomial Regression (degree 2)')
plt.legend()
plt.tight_layout()
plt.savefig('poly_reg_fit.png', dpi=120)
print("Plot saved.")

The red curve should thread smoothly through the scatter — close enough to capture the parabolic shape without snaking through individual noise points.

The Overfitting Trap

Polynomial regression's biggest risk is choosing too high a degree. A degree-15 polynomial can memorize all 80 training points perfectly (RMSE ≈ 0) but will oscillate wildly between them and fail on new data.

The standard diagnostic is a validation curve: plot train and test error as a function of degree.

degrees = range(1, 12)
train_rmses, test_rmses = [], []

for d in degrees:
    m = make_pipeline(PolynomialFeatures(d, include_bias=False), LinearRegression())
    m.fit(X_train, y_train)
    train_rmses.append(np.sqrt(mean_squared_error(y_train, m.predict(X_train))))
    test_rmses.append(np.sqrt(mean_squared_error(y_test,  m.predict(X_test))))

plt.figure(figsize=(7, 4))
plt.plot(degrees, train_rmses, 'o-', label='Train RMSE', color='steelblue')
plt.plot(degrees, test_rmses,  's-', label='Test RMSE',  color='orange')
plt.xlabel('Polynomial Degree')
plt.ylabel('RMSE')
plt.title('Validation Curve: Choosing the Best Degree')
plt.legend()
plt.tight_layout()
plt.savefig('poly_reg_validation_curve.png', dpi=120)
print("Plot saved.")

What you will see:

• Degree 1 (linear): both train and test RMSE are high — underfitting.
• Degree 2: both drop sharply — the model captures the true shape.
• Degree 5+: train RMSE keeps falling, but test RMSE rises — the gap between train and test signals overfitting.

The best degree is where test RMSE is lowest before it starts climbing again. In this example that is degree 2, which matches the true data-generating process.

For a more robust degree selection, use cross-validation instead of a single train/test split.

Using numpy.polyfit (Quick Alternative)

For simple univariate problems, NumPy's polyfit function offers a one-line fit without scikit-learn. It is useful for exploratory analysis but does not integrate with pipelines or cross-validation.

import numpy as np

x = np.array([1, 2, 3, 4, 5], dtype=float)
y = np.array([2.1, 4.0, 9.2, 16.1, 25.0])

# Fit a degree-2 polynomial
# Returns coefficients from highest to lowest degree: [a2, a1, a0]
coefficients = np.polyfit(x, y, 2)
poly = np.poly1d(coefficients)

print("Coefficients (highest to lowest degree):", np.round(coefficients, 3))
print("Prediction at x=6:", round(poly(6), 2))

Expected output:

Coefficients (highest to lowest degree): [ 1.121 -0.939  1.76 ]
Prediction at x=6: 36.5

The fitted curve is approximately y ≈ 1.12x² - 0.94x + 1.76, which is close to the true y = x² relationship in this data. The coefficients are not exactly [1, 0, 0] because there are only five noisy data points to fit.

np.poly1d wraps the coefficients so you can call the polynomial like a function: poly(6) evaluates 1.121(36) - 0.939(6) + 1.76 ≈ 36.5.

Complete Pipeline (All Steps Together)

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

# 1. Generate data (true relationship: y = 0.5x² - x + 2 + noise)
rng = np.random.default_rng(42)
X = rng.uniform(-3, 3, 80).reshape(-1, 1)
y = 0.5 * X.ravel()**2 - X.ravel() + 2 + rng.normal(0, 0.5, 80)

# 2. Split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 3. Build and fit
model = make_pipeline(PolynomialFeatures(degree=2, include_bias=False), LinearRegression())
model.fit(X_train, y_train)

# 4. Evaluate
rmse = np.sqrt(mean_squared_error(y_test, model.predict(X_test)))
r2   = r2_score(y_test, model.predict(X_test))
print(f"Test RMSE: {rmse:.4f}")
print(f"Test R²:   {r2:.4f}")

# 5. Inspect learned coefficients
lr = model.named_steps['linearregression']
pf = model.named_steps['polynomialfeatures']
feature_names = pf.get_feature_names_out(['x'])
for name, coef in zip(feature_names, lr.coef_):
    print(f"  {name}: {coef:.4f}")
print(f"  intercept: {lr.intercept_:.4f}")

Expected output:

Test RMSE: 0.4538
Test R²:   0.9514
  x: -0.9625
  x^2: 0.4909
  intercept: 1.9654

The recovered coefficients (x^2 ≈ 0.49, x ≈ -0.96, intercept ≈ 1.97) closely match the true values (0.5, -1, 2), which confirms the model has learned the correct shape.

Feature Scaling and Polynomial Regression

When working with higher-degree polynomials, feature values grow rapidly — x = 100 produces x² = 10,000 and x³ = 1,000,000. This can make the least-squares problem numerically unstable.

The standard fix is to scale features with StandardScaler inside the pipeline, before the polynomial expansion:

from sklearn.preprocessing import StandardScaler

model = make_pipeline(
    StandardScaler(),
    PolynomialFeatures(degree=3, include_bias=False),
    LinearRegression()
)

Placing StandardScaler first means it is fitted only on training data and applied consistently to test data — no leakage. See Feature Scaling for details on why scaling matters.

Common Pitfalls

Skipping the train/test split. If you fit and evaluate on the same data, high-degree polynomials will appear to perform perfectly while completely failing on new inputs. Always evaluate on held-out data.

Extrapolating outside the training range. Polynomial curves oscillate and diverge outside the range of training data. A model trained on x ∈ [0, 5] can give absurd predictions at x = 10. Linear regression extrapolates more conservatively.

Not scaling before high-degree terms. Large feature values combined with high-degree expansion can produce numerical overflow or ill-conditioned matrices. Use StandardScaler in the pipeline.

Choosing degree by training error alone. Training RMSE always decreases as degree increases. Use test RMSE or cross-validation to find the degree that generalizes.

Forgetting include_bias=False. PolynomialFeatures adds a constant column (intercept term) by default. LinearRegression also adds its own intercept. Passing include_bias=False to PolynomialFeatures avoids the redundant column.

Next Steps

• Linear Regression — the straight-line foundation that polynomial regression builds on
• Multiple Regression — combine many features (polynomial regression applied to multiple inputs produces interaction terms too)
• Train/Test Split — the correct way to evaluate any regression model
• Cross-Validation — more robust than a single split for tuning the polynomial degree
• Feature Scaling — why and how to standardize inputs before fitting