Linear and Polynomial Regression

Overview#

Linear regression is the oldest and most fundamental algorithm in this course — nearly every other supervised method either builds on it directly (logistic regression, regularized regression) or exists partly as a reaction to its limitations (trees, SVMs, neural networks). This notebook goes deep on how linear regression actually works internally: the math behind fitting a line, the assumptions it relies on, how it's optimized (closed-form vs gradient descent), polynomial regression as a natural extension, and a thorough look at when it succeeds and fails.

1. The Algorithm, Precisely#

Linear regression models the target as a weighted sum of the input features, plus a bias term.

Mathematical Formulation
Explain FormulaCopy
y_hat = w_1*x_1 + w_2*x_2 + ... + w_n*x_n + b

Or in vector form:
y_hat = w · x + b

Python
CopyAsk AI
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=200, n_features=3, noise=10, random_state=42)

model = LinearRegression()
model.fit(X, y)

print(f"Weights (coefficients): {model.coef_}")
print(f"Bias (intercept): {model.intercept_:.3f}")

The model's entire "knowledge" is stored in just n_features + 1 numbers — the weights and the bias. This is what makes linear regression so fast to train and so easy to interpret, compared to almost every other algorithm in this course.

2. How the Weights Are Found: Ordinary Least Squares (OLS)#

The weights are chosen to minimize the sum of squared residuals (errors) between predictions and actual values.

Mathematical Formulation
Explain FormulaCopy
Loss (OLS) = sum over all samples of (y_actual - y_predicted)²

Python
CopyAsk AI
# Manually computing OLS to see what's happening under the hood
import numpy as np

X_with_bias = np.column_stack([np.ones(len(X)), X])   # add a column of 1s for the bias term

# The closed-form (normal equation) solution:
# w = (X^T X)^-1 X^T y
w = np.linalg.inv(X_with_bias.T @ X_with_bias) @ X_with_bias.T @ y

print(f"Manually computed weights: {w}")
print(f"scikit-learn's weights:    {[model.intercept_] + list(model.coef_)}")
# These should match almost exactly

Why squared error, specifically?#

text
Ask AICopy
Using SQUARED error (rather than, say, absolute error) has two big
consequences:

1. It penalizes large errors disproportionately more than small ones —
   an error of 10 contributes 100 to the loss, while two errors of 5
   contribute only 25 each (50 total) — this makes OLS sensitive to
   outliers.

2. It has a convenient mathematical property: the loss function is
   smooth and convex, meaning it has a single global minimum that can be
   found EXACTLY, in closed form (the normal equation above) — no
   iterative search required, unlike almost every other algorithm in
   this course.

The closed-form solution vs gradient descent#

Python
CopyAsk AI
# Closed-form (normal equation): exact, but requires inverting a matrix,
# which is computationally expensive for a very large number of features
# — O(n_features^3) complexity for the matrix inversion step

# Gradient descent: iterative, approximate, but scales much better to
# large numbers of features (and is the only practical option once you
# move to more complex models like neural networks)

from sklearn.linear_model import SGDRegressor

sgd_model = SGDRegressor(max_iter=1000, learning_rate="optimal", random_state=42)
sgd_model.fit(X, y)
print(f"SGD weights: {sgd_model.coef_}")   # should converge close to the OLS solution

Mathematical Formulation
Explain FormulaCopy
Gradient descent, conceptually:
1. Start with random weights
2. Compute the gradient (direction of steepest INCREASE) of the loss
   function with respect to each weight
3. Update each weight by taking a small step in the OPPOSITE direction
   (steepest decrease)
4. Repeat until the loss stops meaningfully decreasing

weight_new = weight_old - learning_rate * gradient

For small-to-moderate datasets with a manageable number of features, the closed-form solution (which LinearRegression uses internally) is exact and typically faster; for very large datasets or very high-dimensional data, gradient descent-based methods (SGDRegressor) scale better.

3. The Assumptions Behind Linear Regression#

Linear regression's interpretability and statistical guarantees rely on several assumptions — violating them doesn't necessarily break the model, but it does undermine the validity of things like confidence intervals and hypothesis tests built on top of it.

text
Ask AICopy
1. LINEARITY: the true relationship between features and target is
   (approximately) linear. Polynomial regression (section 5) relaxes
   this specific assumption.

2. INDEPENDENCE: observations are independent of each other. Violated by,
   e.g., time series data where consecutive observations are correlated.

3. HOMOSCEDASTICITY: the variance of the residuals (errors) is roughly
   constant across all levels of the predicted value. "Heteroscedasticity"
   (error variance that grows or shrinks with the prediction) violates this.

4. NORMALITY OF RESIDUALS: the residuals are approximately normally
   distributed — mainly matters for the validity of confidence intervals
   and statistical significance tests, less so for raw prediction accuracy.

5. NO MULTICOLLINEARITY: features aren't too highly correlated with each
   other — severe multicollinearity makes individual coefficient
   estimates unstable and hard to interpret, even if overall predictions
   remain reasonable.

Checking these assumptions in practice#

Python
CopyAsk AI
import matplotlib.pyplot as plt

y_pred = model.predict(X)
residuals = y - y_pred

# Residual plot: should show no obvious pattern (random scatter around 0)
plt.scatter(y_pred, residuals)
plt.axhline(y=0, color="r", linestyle="--")
plt.xlabel("Predicted values")
plt.ylabel("Residuals")
plt.title("Residual Plot — checking homoscedasticity and linearity")
plt.show()

# A funnel shape (widening spread) signals heteroscedasticity
# A curved pattern signals the true relationship isn't actually linear

Python
CopyAsk AI
# Checking multicollinearity with Variance Inflation Factor (VIF)
from statsmodels.stats.outliers_influence import variance_inflation_factor
import pandas as pd

X_df = pd.DataFrame(X, columns=[f"feature_{i}" for i in range(X.shape[1])])
vif_data = pd.DataFrame()
vif_data["feature"] = X_df.columns
vif_data["VIF"] = [variance_inflation_factor(X_df.values, i) for i in range(X_df.shape[1])]
print(vif_data)
# VIF > 5-10 is often used as a rule-of-thumb threshold for concerning
# multicollinearity, though the exact cutoff is somewhat context-dependent

4. Interpreting Coefficients#

Python
CopyAsk AI
import pandas as pd

feature_names = ["square_footage", "num_bedrooms", "age_years"]
coef_df = pd.DataFrame({"feature": feature_names, "coefficient": model.coef_})
print(coef_df)

text
Ask AICopy
Each coefficient represents: "holding all other features constant, how
much does the predicted target change for a one-unit increase in this
feature?"

Example: if the coefficient for square_footage is 150, the model predicts
an increase of $150 in price for each additional square foot — ASSUMING
num_bedrooms and age_years stay fixed.

The "holding constant" caveat, and why scaling matters for comparing coefficients#

Python
CopyAsk AI
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

model_scaled = LinearRegression().fit(X_scaled, y)
print(f"Standardized coefficients: {model_scaled.coef_}")

Raw coefficients aren't directly comparable across features measured on different scales (a coefficient of 150 for square footage vs 5000 for number of bedrooms doesn't mean bedrooms matter 33x less) — standardizing features first makes coefficients directly comparable in terms of "how many standard deviations does the target change per standard deviation change in this feature."

5. Polynomial Regression: Extending to Non-Linear Relationships#

When the true relationship is curved rather than a straight line, polynomial regression adds polynomial terms (x², x³, ...) as additional features — importantly, this is still linear regression, just fit on transformed (expanded) features.

Python
CopyAsk AI
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression

np.random.seed(42)
X_curve = np.sort(np.random.rand(50, 1) * 10, axis=0)
y_curve = 0.5 * X_curve.ravel() ** 2 - 2 * X_curve.ravel() + 5 + np.random.randn(50) * 2

poly_model = make_pipeline(PolynomialFeatures(degree=2), LinearRegression())
poly_model.fit(X_curve, y_curve)

print(f"R² score: {poly_model.score(X_curve, y_curve):.3f}")

Python
CopyAsk AI
# Inspecting what PolynomialFeatures actually does
poly = PolynomialFeatures(degree=2)
X_example = np.array([[2, 3]])
X_poly = poly.fit_transform(X_example)
print(poly.get_feature_names_out())    # ['1', 'x0', 'x1', 'x0^2', 'x0 x1', 'x1^2']
print(X_poly)                             # [[1, 2, 3, 4, 6, 9]]

PolynomialFeatures generates every polynomial combination of the original features up to the specified degree — including interaction terms (like x0 * x1) between different original features, not just powers of individual ones. The underlying LinearRegression then fits weights for each of these expanded features exactly as before — it's still solving a linear problem, just in a higher-dimensional, transformed feature space.

Why "linear" regression can fit curves: the key insight#

Mathematical Formulation
Explain FormulaCopy
"Linear" in "linear regression" refers to the model being linear IN ITS
PARAMETERS (the weights), not linear in the original input features.

y = w1*x + w2*x^2 + w3*x^3 + b

This is still a LINEAR combination of the (transformed) features
[x, x^2, x^3] — even though the resulting curve, as a function of the
ORIGINAL x, is clearly non-linear.

The overfitting risk of higher-degree polynomials#

Python
CopyAsk AI
degrees = [1, 2, 4, 10, 15]
for degree in degrees:
    poly_model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    poly_model.fit(X_curve, y_curve)
    train_r2 = poly_model.score(X_curve, y_curve)
    print(f"Degree {degree}: training R² = {train_r2:.4f}")

text
Ask AICopy
Degree 1: training R² = 0.7823
Degree 2: training R² = 0.9156
Degree 4: training R² = 0.9203
Degree 10: training R² = 0.9587
Degree 15: training R² = 0.9931   <- near-perfect training fit, but almost
                                      certainly severely overfit — this is
                                      the exact overfitting pattern
                                      covered in the Statistical Learning
                                      Foundations notebook, now visible
                                      concretely with polynomial degree as
                                      the complexity knob

Always validate polynomial degree choice with cross-validation, exactly as with any other hyperparameter — training R² alone will keep improving indefinitely as degree increases, right up until the model is essentially memorizing noise.

6. Regularized Variants: A Pointer Forward#

Plain linear/polynomial regression (OLS) has no built-in defense against overfitting, especially with many features or high polynomial degree. Ridge, Lasso, and Elastic Net (covered in full depth in the Regularization Techniques notebook of the main Machine Learning course) add a penalty term to this same OLS loss function to control that overfitting directly.

Python
CopyAsk AI
from sklearn.linear_model import Ridge

ridge_poly = make_pipeline(PolynomialFeatures(degree=10), Ridge(alpha=10))
ridge_poly.fit(X_curve, y_curve)
# The regularization penalty discourages the extreme coefficients that
# would otherwise let a degree-10 polynomial chase every noisy data point

7. Feature Engineering Beyond Polynomials#

Python
CopyAsk AI
# Log transform — useful when a feature's effect is naturally
# multiplicative rather than additive, or to reduce the influence of
# extreme values
X_log = np.log1p(X)   # log1p(x) = log(1+x), handles zero values safely

# Interaction terms — capturing that two features' COMBINED effect
# differs from the sum of their individual effects
X_interaction = X[:, 0] * X[:, 1]

# Binning a continuous feature into categories, then one-hot encoding —
# useful when the relationship is genuinely stepped rather than smooth

Thoughtful feature engineering often improves a linear model's performance more than switching to a fundamentally different algorithm — since linear regression's whole limitation is the shape of function it can represent, and clever features can reshape the problem to fit that limitation.

8. When Linear Regression Is (and Isn't) the Right Choice#

text
Ask AICopy
Linear regression tends to work well when:
- The true relationship is approximately linear (or can be made so via
  reasonable feature engineering/transforms)
- Interpretability matters — coefficients directly explain feature
  effects, unlike most other algorithms in this course
- The dataset is small-to-moderate, where more flexible models risk
  overfitting badly
- A fast, simple baseline is needed before trying more complex methods

Linear regression struggles when:
- The true relationship is highly non-linear and can't be reasonably
  approximated even with polynomial/engineered features
- There are complex interactions between many features that would
  require an impractical number of manually engineered interaction terms
- Outliers are present and haven't been addressed (squared error loss is
  highly sensitive to them)
- The number of features approaches or exceeds the number of samples
  (needs regularization, at minimum, to remain stable)

Notes#

• OLS minimizes squared error specifically because it yields a smooth, convex loss function with an exact closed-form solution — this is a deliberate mathematical choice, not an arbitrary one.
• "Linear" refers to linearity in the model's parameters, not necessarily in the original features — this is exactly why polynomial regression is still considered a form of linear regression.
• Model assumptions (linearity, homoscedasticity, independence, normality of residuals) matter most for the validity of statistical inference (confidence intervals, p-values) — violating them doesn't automatically ruin raw predictive accuracy, but it does undermine interpretation.
• Polynomial degree is a direct bias-variance dial, exactly like every other complexity hyperparameter covered elsewhere in this course — validate it with cross-validation, never training score alone.
• Coefficients are only directly comparable across features when those features have been scaled to the same range first.

Best Practices Checklist#

• Check residual plots for non-random patterns (signals of non-linearity or heteroscedasticity) rather than assuming assumptions hold
• Scale features before comparing coefficient magnitudes across different features
• Use cross-validation to select polynomial degree, never training R² alone
• Check for multicollinearity (e.g. via VIF) when coefficient interpretation matters, not just prediction accuracy
• Consider Ridge/Lasso/Elastic Net whenever using polynomial features beyond a low degree, to control overfitting
• Remember squared-error loss is sensitive to outliers — investigate or address extreme values before fitting
• Use linear regression as a fast, interpretable baseline even when you expect to eventually need a more complex model
• Engineer meaningful features (interactions, transforms) before assuming linear regression "can't handle" a non-linear problem