Linear regression

Linear regression is where supervised learning’s whole conceptual toolkit appears for the first time: a model class, a loss, an estimator, and three different lenses that each make something visible the others hide. Most treatments pick one lens. A practitioner needs all three, because debugging a regression means knowing which lens the symptom belongs to.

The model and the loss

Data: $n$ pairs $(x_i, y_i)$ with $x_i \in \mathbb{R}^d$. Model: predictions linear in the parameters,

($X$‘s first column is all-ones for the intercept). “Linear” constrains the parameters, not the inputs — replace the columns of $X$ with polynomials, splines or interactions and everything in this lesson survives untouched; that flexibility is exactly what the widget below exploits. The estimator is ordinary least squares (OLS):

Lens 1 — optimization: the normal equations

The objective is a convex quadratic (its Hessian is $X^\top X \succeq 0$ — the convexity lesson of the Optimization module), so setting the gradient to zero is sufficient:

— the normal equations. If $X$ has full column rank, $\hat\beta = (X^\top X)^{-1} X^\top y$ exists and is unique. Two computational regimes follow:

The Optimization module is not a separate subject here: the GD route’s iteration count scales with $\kappa(X^\top X)$, and $\kappa(X^\top X) = \kappa(X)^2$ — conditioning is squared by forming the normal equations, which is why numerical libraries prefer QR and why unstandardized features (one in $[0,1]$, one in $[0, 10^4]$) can make gradient-based fitting absurdly slow. If you did exercise 5 of the gradient-descent lesson, this is the payoff.

Lens 2 — geometry: prediction is projection

$X\beta$ ranges over the column space of $X$ — a $d$-dimensional plane inside $\mathbb{R}^n$. Minimizing $\lVert y - X\beta \rVert$ means finding the point of that plane closest to $y$: orthogonal projection. The fitted values are

with $H$ the hat matrix (it puts the hat on $y$): symmetric, idempotent ($H^2 = H$), a projection. The residual $e = y - \hat y = (I - H)y$ is the perpendicular component, which forces

residuals are exactly orthogonal to every column of X — mean-zero (the intercept column), and uncorrelated with every feature, by construction. This is why “I see structure in my residuals” is such a powerful diagnostic: any pattern the residuals show against a variable is structure the model could not represent, since anything representable was projected out. And $\operatorname{tr}(H) = d$ counts the model’s degrees of freedom — a number that returns, fractionally, when ridge regression bends this projection in the next lesson.

Lens 3 — statistics: when is OLS the right estimator?

Assume the data really are $y = X\beta^\ast + \varepsilon$ with $\mathbb{E}[\varepsilon] = 0$, $\operatorname{Cov}(\varepsilon) = \sigma^2 I$. Then:

• Unbiasedness: $\mathbb{E}[\hat\beta] = \beta^\ast$ — substitute and use $(X^\top X)^{-1}X^\top X = I$.
• Covariance: $\operatorname{Cov}(\hat\beta) = \sigma^2 (X^\top X)^{-1}$. Read this matrix: nearly-collinear columns make $X^\top X$ nearly singular, so coefficient variances blow up along the associated directions — the statistical face of the same ill-conditioning that slowed gradient descent. One pathology, three lenses.
• Gauss–Markov: among all linear unbiased estimators, OLS has minimal variance — it is BLUE. Note what is not assumed: normality. Gaussian errors add only the exact finite-sample distributions ($t$-tests on coefficients) and make OLS equal to maximum likelihood, since minimizing squared error is maximizing a Gaussian log-likelihood.

The fine print is where real datasets live: heteroscedastic errors break the covariance formula (weighted least squares or robust/sandwich errors), correlated errors break it worse (time series — its own track), and squared loss’s sensitivity to outliers is a choice, inherited from the Gaussian assumption, swappable for Huber or absolute loss at the price of closed forms. And Gauss–Markov’s “best” is conditional on unbiased — the next lesson exists precisely because a little bias, bought wisely, can be a bargain.

import numpy as np

# The three computational routes to the same estimator
beta_ne, *_ = np.linalg.lstsq(X, y, rcond=None)        # QR/SVD route — use this
beta_ch = np.linalg.solve(X.T @ X, X.T @ y)            # normal equations (κ²!)

beta_gd = np.zeros(X.shape[1])                          # the optimization route
for _ in range(2000):
    beta_gd -= 1e-2 * X.T @ (X @ beta_gd - y) / len(y)

Exercises