Polynomial Regression & Nonlinear Relationships

1. Introduction: The Limits of Straight Lines

In the previous lesson, we analyzed Linear Regression, establishing how Ordinary Least Squares (OLS) fits a straight line or flat hyperplane ($\hat{y} = wx + b$) through continuous observation points. The fundamental premise of linear regression is that a unit increase in the input feature ($x$) consistently produces a constant marginal shift in the target output ($y$).

However, nature, physics, and human behavior are rarely straight lines. In real-world domains:

• Diminishing Marginal Returns: Early marketing expenditures trigger rapid user growth, but customer acquisition gradually plateaus as the addressable market saturates.
• Accelerating Physical Laws: A vehicle traveling at $100\text{ km/h}$ requires vastly more than double the braking distance of a vehicle traveling at $50\text{ km/h}$, because kinetic energy scales quadratically with velocity ($E_k \propto v^2$).
• U-Shaped Seasonal Responses: Electrical grid consumption surges during freezing winter lows (heating) and sweltering summer peaks (air conditioning), forming an inverted parabolic curve against ambient temperature.

When empirical data curves, forcing an inflexible straight line through it produces systematic modeling failure—a condition formally known as Underfitting (High Bias). To model curvature without abandoning the computational efficiency and interpretability of linear algebra, we introduce Polynomial Regression.

2. Real-World Nonlinear Phenomena: A Comparative Taxonomy

The table below highlights classical domains where standard linear models fail and non-linear or polynomial representations become mandatory:

3. Mathematical Formulation of Polynomial Regression

Recall the single-feature Simple Linear Regression equation from Lesson 15.13:

Polynomial Regression generalizes this formulation by modeling the target variable $\hat{y}$ as an $n$-th degree polynomial in the input feature $x$:

Where $x$ is the original raw input feature, $b$ is the scalar intercept, $w_j$ are the learned weights (coefficients), and $d$ represents the Polynomial Degree.

4. The Linear Duality: Why Polynomial Regression is "Linear"

A frequent point of conceptual confusion among machine learning students is: "If the predicted curve bends like a parabola or wave, why is it categorized under Linear Models?"

To visualize this duality, consider transforming the single non-linear feature $x$ into new synthetic features $z_1, z_2, \dots, z_d$:

Substituting these newly synthesized feature tokens into our polynomial equation yields:

This equation is 100% equivalent to Multiple Linear Regression! The partial derivative with respect to each weight remains constant with respect to that weight ($\frac{\partial \hat{y}}{\partial w_j} = z_j$). Therefore, we can solve for the optimal weights using standard Ordinary Least Squares (OLS) linear algebra:

Where $\mathbf{Z}$ is the transformed design matrix containing the polynomial powers of the original training observations.

5. Multivariate Polynomial Regression & Interaction Features

When a dataset contains multiple input features ($x_1, x_2, \dots, x_p$), expanding to degree $d$ synthesizes not only pure polynomial powers ($x_1^2, x_2^2$), but also cross-product Interaction Features ($x_1 x_2$).

6. The Bias-Variance Dilemma: Underfitting vs. Overfitting

The polynomial degree ($d$) is a direct throttle on the structural complexity and flexibility of the model hypothesis space:

7. Hyperparameter Tuning: Selecting the Optimal Degree

Because Polynomial Degree ($d$) cannot be learned directly through analytical matrix inversion, it acts as a critical Hyperparameter. We evaluate model performance across ascending degrees using an independent Validation Set or K-Fold Cross-Validation.

Notice the classic U-shaped validation curve: as degree increases, training error drops monotonically, but validation error drops to a minimum at $d=3$ before exploding upward.

8. End-to-End Scikit-Learn Python Implementation

We build a production-grade machine learning pipeline generating a synthetic quadratic dataset ($y = 3 + 2x + 1.5x^2 + \epsilon$), fitting polynomial models, and measuring out-of-sample metrics.

# Step 1: Import Core Scientific & Machine Learning Libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

# Step 2: Synthesize Non-Linear Quadratic Ground Truth: y = 3 + 2x + 1.5x² + ε
rng = np.random.default_rng(seed=42)
X = np.linspace(-3.0, 3.0, 150).reshape(-1, 1)
noise = rng.normal(loc=0.0, scale=2.0, size=len(X))
y = 3.0 + 2.0 * X[:, 0] + 1.5 * (X[:, 0] ** 2) + noise

# Step 3: Partition Dataset into 80% Train and 20% Test Sets
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Step 4: Construct Leak-Free Polynomial Pipeline (Degree=2)
# Note: include_bias=False prevents duplicate intercept columns (fit_intercept handles b)
poly_pipeline = make_pipeline(
    PolynomialFeatures(degree=2, include_bias=False),
    LinearRegression()
)

# Step 5: Fit the Pipeline Exclusively on Training Data
poly_pipeline.fit(X_train, y_train)

# Step 6: Generate Out-of-Sample Predictions on Unseen Test Data
y_pred = poly_pipeline.predict(X_test)

# Step 7: Quantify Generalization Error Metrics
mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_test, y_pred)

print(f"=== POLYNOMIAL REGRESSION (DEGREE 2) TEST PERFORMANCE ===")
print(f"MAE  : {mae:.3f}")
print(f"MSE  : {mse:.3f}")
print(f"RMSE : {rmse:.3f}")
print(f"R²   : {r2:.4f}")

# Inspect Recovered Mathematical Parameters
lin_reg = poly_pipeline.named_steps['linearregression']
print(f"\nRecovered Intercept (b) : {lin_reg.intercept_:.3f} (True: 3.0)")
print(f"Recovered Weights (w1, w2): {lin_reg.coef_[0]:.3f}, {lin_reg.coef_[1]:.3f} (True: 2.0, 1.5)")

9. Cross-Validation Hyperparameter Tuning Across Degrees

To avoid overfitting to a single arbitrary train-test split, we execute 5-Fold Cross-Validation across degrees $d \in \{1, 2, 3, 4, 5\}$:

from sklearn.model_selection import KFold, cross_val_score

degrees = [1, 2, 3, 4, 5]
kf = KFold(n_splits=5, shuffle=True, random_state=42)

print(f"Degree | Mean CV RMSE | Std Dev  | Diagnostic Status")
print("------------------------------------------------")

for d in degrees:
    model = make_pipeline(
        PolynomialFeatures(degree=d, include_bias=False),
        StandardScaler(),
        LinearRegression()
    )
    # cross_val_score returns negative MSE per scikit-learn optimization convention
    neg_mse = cross_val_score(model, X_train, y_train, cv=kf, scoring="neg_mean_squared_error")
    rmse_scores = np.sqrt(-neg_mse)
    
    mean_rmse = np.mean(rmse_scores)
    std_rmse = np.std(rmse_scores)
    
    status = "Optimal" if d == 2 else ("Underfitting" if d == 1 else "Mild Overfitting")
    print(f"  d={d}  |   {mean_rmse:6.3f}    |  ±{std_rmse:.3f}  | {status}")

10. Why Preprocessing Pipelines Prevent Severe Leakage & Instability

Building polynomial models without Scikit-Learn's `Pipeline` exposes ML systems to two critical failure modes:

11. Five Critical Limitations of Polynomial Regression

12. Taming Polynomial Oscillations: Ridge ($L_2$) & Lasso ($L_1$)

When domain requirements demand high-degree features to capture localized nuance, standard OLS will overfit. The primary defense is Regularization—penalizing exploding parameter magnitudes:

13. Architectural Decision Framework: When to Use Polynomial Models

15. Summary & Key Takeaways