How Polynomial Regression Extends Linear Regression Using Feature Engineering

Polynomial regression extends linear regression by engineering new features that are powers and interactions of original variables, then applying the same linear algebra solver to this expanded feature space to model nonlinear relationships.

Polynomial regression extends linear regression using feature engineering by transforming raw inputs into higher-order terms while keeping the underlying learning algorithm unchanged. In the ai-engineering-from-scratch repository, this extension is implemented through modular utility functions that generate polynomial features before fitting. By expanding the feature space to include quadratic and interaction terms, the model captures complex curves using the same ordinary least squares approach applied to an augmented design matrix.

The Mathematical Mechanism

Linear regression models a target y as a weighted sum of original input features x:

ŷ = w₀ + w₁x₁ + w₂x₂ + ... + wₚxₚ

This formulation captures only linear relationships between inputs and outputs. To model nonlinear patterns, polynomial regression augments the feature space by creating new features that contain powers and pairwise products of the original variables.

When you fit a linear model to this expanded set, the overall system represents polynomial relationships while utilizing identical linear-algebra machinery. The key insight is that polynomial regression is linear regression applied to a transformed feature space.

Implementation in ai-engineering-from-scratch

The repository implements this extension through two core components that separate feature engineering from model fitting.

Feature Engineering with polynomial_features

Located in phases/02-ml-fundamentals/08-feature-engineering/code/features.py (lines 39-48), the polynomial_features function generates the expanded feature vector containing original features, their squares, and pairwise interaction terms.

The function performs three operations on a raw input vector row = (x₁, x₂, …, xₙ):

1. Copies the original components (result = list(row))
2. Appends each component squared (xᵢ²)
3. Appends each distinct pairwise product (xᵢ·xⱼ)

This yields a feature vector of length n + n + n·(n‑1)/2, covering all terms up to degree 2. The implementation is generic; increasing the degree parameter adds higher-order powers recursively.

Model Fitting with fit_polynomial and predict_polynomial

The fitting logic resides in phases/02-ml-fundamentals/10-bias-variance/code/bias_variance.py (lines 17-31). The fit_polynomial function constructs a design matrix by stacking power transformations:

X = np.column_stack([x_train ** d for d in range(degree + 1)])

This matrix contains columns for each power of x, including the intercept column (x⁰). The function then solves the normal equation w = (XᵀX)⁻¹Xᵀy or a ridge-regularized variant—exactly the same algorithm used for ordinary linear regression, but operating on the augmented matrix.

For inference, predict_polynomial rebuilds the design matrix for new inputs and computes X @ w. Because the matrix encodes nonlinear terms, the output curve bends to capture curvature and interactions impossible for a plain linear predictor.

Practical Code Example

The following example demonstrates the complete workflow using the repository's implementation:

import numpy as np
from phases.02_ml_fundamentals.08_feature_engineering.code.features import polynomial_features
from phases.02_ml_fundamentals.10_bias_variance.code.bias_variance import fit_polynomial, predict_polynomial

# 1. Create synthetic nonlinear data

rng = np.random.RandomState(0)
x = rng.uniform(-2, 2, size=30)
y = np.sin(1.5 * x) + 0.5 * x + rng.normal(0, 0.5, size=30)

# 2. Build polynomial design matrix (degree = 2)

X_poly = np.column_stack([x ** d for d in range(3)])

# 3. Fit polynomial regression

w = fit_polynomial(x, y, degree=2)

# 4. Predict on dense grid

x_test = np.linspace(-2.5, 2.5, 200)
y_pred = predict_polynomial(x_test, w)

# 5. Visualize results

import matplotlib.pyplot as plt
plt.scatter(x, y, label='data')
plt.plot(x_test, y_pred, color='red', label='poly-reg (deg=2)')
plt.legend()
plt.show()

This example generates a sinusoidal dataset that linear regression cannot capture. By expanding to quadratic terms (x²), the fitted model bends to approximate the underlying function. The same fit_polynomial routine works for any degree; increasing the degree merely adds more columns to the design matrix, enabling higher-order polynomial fits.

Architectural Benefits

Separating feature engineering from model fitting provides three distinct advantages:

• Reusability: The same fit_polynomial logic supports any degree without code changes
• Interpretability: The transformation is explicit in polynomial_features, making the feature space transparent
• Efficiency: The implementation leverages NumPy broadcasting for vectorized matrix construction, avoiding Python loops during training

According to the ai-engineering-from-scratch source code, this separation mirrors scikit-learn's transformer/estimator pattern while maintaining educational clarity through explicit implementation of the design matrix construction.

Summary

• Polynomial regression extends linear regression by engineering higher-order features from raw inputs, applying the same linear solver to a transformed space.
• The polynomial_features function in features.py (lines 39-48) generates squares and interaction terms, expanding the feature vector length from n to n + n + n·(n‑1)/2 for degree 2.
• fit_polynomial and predict_polynomial in bias_variance.py (lines 17-31) implement the linear algebra solution on the augmented design matrix using ordinary least squares or ridge regularization.
• This approach captures nonlinear relationships and feature interactions while retaining the computational simplicity and interpretability of linear models.

Frequently Asked Questions

What is the relationship between polynomial features and the design matrix?

Polynomial features become the columns of the design matrix. In fit_polynomial, the function constructs this matrix by stacking power transformations (x**d for each degree), creating a column for the intercept, linear terms, quadratic terms, and so on. The linear regression solver then treats these engineered columns exactly like native features, computing weights that minimize squared error across the expanded space.

Why does polynomial regression use the same solver as linear regression?

Polynomial regression uses the same solver because it is linear regression applied to a nonlinear transformation of the inputs. Once the features are engineered to include powers and interactions, the relationship between these new features and the target is linear. The normal equation w = (XᵀX)⁻¹Xᵀy remains valid regardless of whether the columns of X represent raw measurements or polynomial combinations thereof.

How does the repository handle higher-order polynomials beyond degree 2?

The repository handles higher orders through the degree parameter in both polynomial_features and fit_polynomial. Increasing the degree simply extends the range of the loop that generates powers: range(degree + 1) creates columns for x⁰ through x^degree. For multivariate inputs, polynomial_features would generate all cross-terms up to the specified degree, though the current implementation focuses on quadratic expansions for educational clarity.

What are the computational trade-offs of polynomial feature engineering?

Polynomial feature engineering increases the feature dimensionality from n to approximately n^d for degree d, raising both memory consumption and computational cost for the matrix inversion step. The ai-engineering-from-scratch implementation mitigates this through NumPy vectorization, but users must balance model flexibility against the risk of overfitting that accompanies high-dimensional polynomial spaces.

{
  "mcpServers": {
    "instagit": {
      "command": "npx",
      "args": ["-y", "instagit@latest"]
    }
  }
}