Historical Context & Motivation
The story of linear regression begins long before the matrix notation we use today, rooted in astronomers' and mathematicians' desire to fit curves to observational data. For most of the nineteenth century, practitioners solved normal equations by hand, manipulating individual sums and cross-products for each new predictor added to a model. As the number of predictors grew, this approach became not merely tedious but error-prone, motivating a search for a more systematic language. The adoption of matrix algebra in the twentieth century unified these disparate calculations into a single, compact framework that could be automated on digital computers. Today, every major statistical software package—R, Python's scikit-learn, SAS, Stata—solves regression problems using matrix operations at its core.
The central question this lesson addresses is deceptively simple: given a response vector and a set of predictors, how can we express the entire regression problem—model, estimation, inference, and geometry—in a single, unified matrix equation? The payoff is enormous: derivations that span pages in scalar notation collapse to a few lines, properties like unbiasedness and efficiency emerge directly from matrix identities, and the geometric interpretation of projection provides an intuition that transcends any particular dataset.
Core Principles & Definitions
Before writing a single equation, it is essential to understand the conceptual pillars that support the matrix formulation of linear regression. Each principle below maps a familiar idea from introductory regression into the language of linear algebra, enabling us to handle models with an arbitrary number of predictors in a unified way. Mastery of these principles is prerequisite to the derivations and proofs that follow in subsequent sections.
The Design Matrix X
The Parameter Vector β
The Error Vector ε
Orthogonal Projection
The Hat Matrix H
Geometric Interpretation of OLS
The most profound insight in the matrix treatment of regression is geometric: ordinary least squares finds the point in the column space of X that is closest to y in Euclidean distance. The diagram below illustrates this projection in a simplified three-dimensional setting where the column space of X is a two-dimensional plane. The response vector y sits outside this plane, the fitted value vector ŷ is its perpendicular projection onto the plane, and the residual vector ê = y − ŷ is perpendicular to the plane. This orthogonality condition, expressed algebraically as Xᵀê = 0, is the geometric heart of the normal equations.
This geometric picture clarifies many properties that would otherwise require separate algebraic proofs. For instance, the Pythagorean theorem applied to the right triangle formed by y, ŷ, and ê yields the ANOVA decomposition ‖y‖² = ‖ŷ‖² + ‖ê‖², which in scalar notation is SST = SSR + SSE. The coefficient of determination R² is simply the squared cosine of the angle between y and ŷ. Adding a new predictor extends the column space, so ŷ can only move closer to y, which explains why R² never decreases with additional predictors.
Mathematical Framework
We now formalize the algebraic structure. All derivations below assume the classical linear model y = Xβ + ε with rank(X) = p ≤ n, E(ε) = 0, and Var(ε) = σ²Iₙ. Under these conditions, we derive the OLS estimator, its distribution, and associated variance.
Deriving the OLS Estimator
We seek β̂ that minimizes the sum of squared residuals S(β) = (y − Xβ)ᵀ(y − Xβ). Expanding, we have S(β) = yᵀy − 2βᵀXᵀy + βᵀXᵀXβ. Taking the gradient with respect to β and setting it to zero gives ∂S/∂β = −2Xᵀy + 2XᵀXβ = 0, which yields the normal equations.
Properties of the OLS Estimator
Because β̂ = (XᵀX)⁻¹Xᵀy = (XᵀX)⁻¹Xᵀ(Xβ + ε) = β + (XᵀX)⁻¹Xᵀε, we immediately see that E(β̂) = β, confirming unbiasedness. The variance-covariance matrix of β̂ follows from the linear transformation of ε.
The Hat Matrix, Fitted Values, and Residual Structure
Substituting the OLS estimator back into the model gives ŷ = Xβ̂ = X(XᵀX)⁻¹Xᵀy = Hy, where H = X(XᵀX)⁻¹Xᵀ is the hat matrix (so named because it "puts the hat on y"). This matrix has remarkable algebraic properties that control the behavior of residuals, leverage, and influence in regression diagnostics.
Leverage and the Diagonal of H
The diagonal elements of H, denoted hᵢᵢ, satisfy 0 ≤ hᵢᵢ ≤ 1 and Σhᵢᵢ = p. The value hᵢᵢ measures how much observation i's predictor values differ from the centroid of the design; large hᵢᵢ (conventionally hᵢᵢ > 2p/n) flags that observation as a high-leverage point. Since Var(êᵢ) = σ²(1 − hᵢᵢ), high-leverage observations have small residual variance, making raw residuals unreliable for outlier detection. This motivates the use of studentized residuals rᵢ = êᵢ / (s√(1 − hᵢᵢ)), which have approximately unit variance under the model assumptions.
Worked Example: Two-Predictor Regression
Consider a small dataset with n = 4 observations and p = 3 columns (intercept plus two predictors). We will compute β̂, the hat matrix H, fitted values, residuals, and the unbiased variance estimate entirely in matrix form.
Assumptions, Diagnostics, and When OLS Fails
The elegance of the matrix framework comes with conditions that must be carefully verified in practice. Violations of these assumptions do not invalidate the algebra—the OLS formula still produces numbers—but the optimality properties (unbiasedness, efficiency, valid inference) may break down. The table below summarizes each assumption, its matrix expression, what happens when it fails, and the standard diagnostic tool.
| Assumption | Matrix Expression | Consequence of Violation | Diagnostic |
|---|---|---|---|
| Linearity | E(y) = Xβ | β̂ is biased; predictions systematically wrong | Residual vs. fitted plot; partial regression plots |
| Full rank | rank(X) = p | (XᵀX) is singular; β̂ not unique | VIF > 10; condition number of XᵀX |
| Spherical errors | Var(ε) = σ²Iₙ | β̂ still unbiased but no longer BLUE; SEs wrong | Breusch–Pagan test; Durbin–Watson test |
| Normality | ε ~ N(0, σ²Iₙ) | t and F tests invalid in small samples (OK asymptotically) | Q–Q plot of residuals; Shapiro–Wilk test |
| No influential outliers | hᵢᵢ bounded; Cook's D small | Single observations dominate β̂ | Leverage hᵢᵢ; Cook's distance; DFFITS |
From OLS to GLS and Generalized Linear Models
The matrix formulation of OLS serves as the conceptual launching pad for virtually every advanced regression method. When the assumption Var(ε) = σ²Iₙ is replaced by Var(ε) = σ²Ω for a known positive-definite Ω, we arrive at Generalized Least Squares (GLS). When the response is non-Gaussian and the mean is linked to the linear predictor through a nonlinear function, we enter the domain of Generalized Linear Models (GLMs). The table below compares these frameworks.
| Feature | OLS | GLS | GLM (IRLS) |
|---|---|---|---|
| Model | y = Xβ + ε | y = Xβ + ε, Var(ε) = σ²Ω | g(μ) = Xβ, y ~ EDF |
| Estimator | (XᵀX)⁻¹Xᵀy | (XᵀΩ⁻¹X)⁻¹XᵀΩ⁻¹y | Iterative: (XᵀWX)⁻¹XᵀWz |
| Variance of β̂ | σ²(XᵀX)⁻¹ | σ²(XᵀΩ⁻¹X)⁻¹ | (XᵀWX)⁻¹ at convergence |
| Error structure | Homoscedastic, uncorrelated | Heteroscedastic and/or correlated | Variance depends on mean via V(μ) |
| Key insight | Projection onto C(X) | Projection in Ω⁻¹-inner product space | Iteratively re-weighted projection |
The progression from OLS to GLS to GLMs illustrates a recurring theme in statistics: each generalization relaxes one assumption while preserving the same matrix-algebraic skeleton. GLS replaces Iₙ with Ω⁻¹ inside the quadratic form, effectively "whitening" the errors before applying OLS. GLMs go further by iterating the GLS step with updated weights W that depend on the current estimate of the mean, converging to the maximum likelihood estimate for exponential-family responses. In each case, the hat matrix, leverage, and residual diagnostics have direct analogs expressed in the weighted inner product, so the geometric intuition of OLS carries forward intact.
Practice Problems
Summary
The classical linear model y = Xβ + ε compactly encodes a system of n equations with p unknowns. By minimizing the squared Euclidean norm ‖y − Xβ‖², we derive the OLS estimator β̂ = (XᵀX)⁻¹Xᵀy from the normal equations XᵀXβ̂ = Xᵀy. The hat matrix H = X(XᵀX)⁻¹Xᵀ projects y onto the column space of X, producing the fitted values ŷ = Hy, while the residuals ê = (I − H)y are orthogonal to every column of X.
Under the Gauss–Markov assumptions (E(ε) = 0, Var(ε) = σ²Iₙ), OLS is BLUE with Var(β̂) = σ²(XᵀX)⁻¹ and σ² is estimated unbiasedly by s² = êᵀê/(n − p). The diagonal elements hᵢᵢ of H measure leverage, and the heterogeneous residual variances Var(êᵢ) = σ²(1 − hᵢᵢ) motivate studentized residuals for diagnostic use. This matrix framework extends naturally to GLS (when errors are heteroscedastic or correlated) and to GLMs (when the response is non-Gaussian), preserving the geometric projection interpretation throughout.