Linear Regression and Least Squares

How should we predict a student's exam score from hours studied? Or a house's price from its square footage? These questions share a common structure. A response variable, we believe, depends approximately linearly on one or more explanatory variables. This relationship is contaminated by noise we cannot control. Linear regression formalizes this idea. Least squares gives the canonical method for extracting the best-fitting line — or hyperplane — from noisy data. The method is over two centuries old. It traces to Legendre and Gauss's work on fitting orbits to astronomical observations. Yet it remains the single most used tool in applied statistics, machine learning, econometrics, and the sciences. The core insight is deceptively simple. Among all lines (or hyperplanes) we could draw through a scatter of points, choose one. Pick the line that minimizes the total squared vertical distance from the data to the line. This criterion is not arbitrary. It is the maximum-likelihood choice under Gaussian noise. It has an elegant closed-form solution via calculus or projections. It also decomposes cleanly into geometric pieces — fitted values and residuals — that are orthogonal to one another. That orthogonality is the structural fact underlying nearly every downstream computation. This includes variance decompositions, hypothesis tests, and confidence intervals. This article builds the linear model from its precise probabilistic definition. It derives the least-squares estimator via a normal-equations argument and a geometric projection argument. It also proves the key algebraic identities relating fitted values, residuals, and the design matrix. We keep the treatment general, allowing an arbitrary number of predictors. Every claim is grounded in explicit computation on worked numerical examples with real, unequal numbers. No step is left to intuition alone. The article closes by isolating exactly what least squares gives you for free. This includes unbiasedness and a specific closed-form variance. It also isolates exactly what least squares does not give you. This sets up the natural next question. Under what conditions is the least-squares estimator actually optimal among unbiased estimators? That question is answered by The Gauss-Markov Theorem, which we foreshadow but do not yet prove.

Key takeaways
  • The linear model assumes Yi=xi⊤β+εiY_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i with mean-zero, uncorrelated errors of common variance σ2\sigma^2.
  • The least-squares estimator β^\hat{\boldsymbol{\beta}} minimizes the sum of squared residuals ∑i(Yi−xi⊤β)2\sum_i (Y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2.
  • When X⊤XX^\top X is invertible, the unique minimizer is β^=(X⊤X)−1X⊤Y\hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{Y}, the normal equations.
  • Fitted values Y^=Xβ^\hat{\boldsymbol{Y}} = X\hat{\boldsymbol{\beta}} are the orthogonal projection of Y\boldsymbol{Y} onto the column space of XX; residuals ε^=Y−Y^\hat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \hat{\boldsymbol{Y}} are orthogonal to that column space.
  • β^\hat{\boldsymbol{\beta}} is unbiased for β\boldsymbol{\beta} and has covariance matrix σ2(X⊤X)−1\sigma^2 (X^\top X)^{-1}, under the stated model alone.
  • Key idea: least squares turns a statistical estimation problem into a geometric projection problem, and the resulting orthogonality of fitted values and residuals drives everything downstream.

Where this lives today

Least squares regression underlies econometric estimation of demand curves, wage equations, and macroeconomic relationships, where it remains the default first tool before more elaborate causal-inference methods. In physics and engineering, it is the workhorse of curve fitting and calibration — from fitting sensor drift models to processing GPS trilateration data, which still traces its ancestry to Gauss's orbit-fitting problem. In machine learning, ordinary least squares is the exact solution to training a linear model under squared-error loss, and it is the special case (ℓ2\ell_2 penalty =0=0) of ridge regression, whose penalized normal equations (X⊤X+λI)β^=X⊤Y(X^\top X + \lambda I)\hat{\boldsymbol\beta}=X^\top\boldsymbol Y directly generalize the formula proved above. In signal processing, the same projection idea appears as the least-squares channel estimator and as the Kalman filter's update step. In statistics more broadly, the orthogonal decomposition ∥Y∥2=∥Y^∥2+∥ε^∥2\|\boldsymbol Y\|^2 = \|\hat{\boldsymbol Y}\|^2+\|\hat{\boldsymbol\varepsilon}\|^2 proved here is the seed of the ANOVA sum-of-squares decomposition and of the FF-test for nested linear models, connecting directly to Standard Tests: t, Chi-Squared, and F. In computational science, the numerically stable way to solve the normal equations — via QRQR or singular value decomposition rather than explicit matrix inversion — is standard practice in every major statistical software package.

Ask the tutor: How does the Gauss–Markov theorem upgrade "unbiased with covariance σ2(X⊤X)−1\sigma^2(X^\top X)^{-1}" into "the best linear unbiased estimator" — and what exactly does "best" mean there?

DefinitionThe linear model, design matrix, and least-squares estimator

Let n≥pn \ge p and suppose we observe a response Y1,…,Yn∈RY_1, \dots, Y_n \in \mathbb{R} together with, for each ii, a vector of predictors xi=(xi1,…,xip)⊤∈Rp\boldsymbol{x}_i = (x_{i1}, \dots, x_{ip})^\top \in \mathbb{R}^p (one coordinate is typically the constant xi1≡1x_{i1} \equiv 1, for an intercept). The linear model postulates

Yi=xi⊤β+εi,i=1,…,n,Y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i, \qquad i = 1, \dots, n,

for an unknown parameter vector β=(β1,…,βp)⊤∈Rp\boldsymbol{\beta} = (\beta_1, \dots, \beta_p)^\top \in \mathbb{R}^p and unobserved errors ε1,…,εn\varepsilon_1, \dots, \varepsilon_n satisfying

E[εi]=0,Var⁡(εi)=σ2<∞,Cov⁡(εi,εj)=0 (i≠j).\mathbb{E}[\varepsilon_i] = 0, \quad \operatorname{Var}(\varepsilon_i) = \sigma^2 < \infty, \quad \operatorname{Cov}(\varepsilon_i, \varepsilon_j) = 0 \ (i \ne j).

In matrix form, with Y=(Y1,…,Yn)⊤\boldsymbol{Y} = (Y_1,\dots,Y_n)^\top, ε=(ε1,…,εn)⊤\boldsymbol{\varepsilon} = (\varepsilon_1,\dots,\varepsilon_n)^\top, and design matrix XX the n×pn \times p matrix whose ii-th row is xi⊤\boldsymbol{x}_i^\top,

Y=Xβ+ε,E[ε]=0,Cov⁡(ε)=σ2In.\boldsymbol{Y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \mathbb{E}[\boldsymbol{\varepsilon}] = \boldsymbol{0}, \qquad \operatorname{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 I_n.

We assume XX has full column rank pp (equivalently X⊤XX^\top X is invertible), so no predictor is an exact linear combination of the others. The least-squares estimator of β\boldsymbol{\beta} is any minimizer of the residual sum of squares

S(b)=∥Y−Xb∥2=∑i=1n(Yi−xi⊤b)2,b∈Rp.S(\boldsymbol{b}) = \|\boldsymbol{Y} - X\boldsymbol{b}\|^2 = \sum_{i=1}^n \bigl(Y_i - \boldsymbol{x}_i^\top \boldsymbol{b}\bigr)^2, \qquad \boldsymbol{b} \in \mathbb{R}^p.

Given the estimator β^\hat{\boldsymbol{\beta}}, the fitted values are Y^=Xβ^\hat{\boldsymbol{Y}} = X\hat{\boldsymbol{\beta}} and the residuals are ε^=Y−Y^\hat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \hat{\boldsymbol{Y}}.

Intuition

Picture data points scattered in a plane, one predictor and one response. Least squares asks: among all straight lines, which one makes the vertical gaps between the line and the points — the residuals — collectively as small as possible, in the sense of minimizing the sum of their squares? Squaring serves two purposes: it makes positive and negative deviations count equally (no cancellation), and it penalizes large deviations disproportionately, pulling the fit toward points that are far off much more strongly than toward points that are close. There is a second, more structural way to see the same object. Stack the nn observations into the vector Y∈Rn\boldsymbol{Y} \in \mathbb{R}^n. The set of all vectors of the form XbX\boldsymbol{b}, as b\boldsymbol{b} ranges over Rp\mathbb{R}^p, is a pp-dimensional subspace C(X)⊂Rn\mathcal{C}(X) \subset \mathbb{R}^n — the column space of XX. Minimizing ∥Y−Xb∥2\|\boldsymbol{Y} - X\boldsymbol{b}\|^2 over b\boldsymbol{b} is exactly the problem of finding the point of C(X)\mathcal{C}(X) closest to Y\boldsymbol{Y} in ordinary Euclidean distance. Elementary geometry says the closest point is the orthogonal projection of Y\boldsymbol{Y} onto C(X)\mathcal{C}(X), and the connecting vector — the residual — is perpendicular to the entire subspace. This single geometric fact, proved rigorously below, is the source of every algebraic identity involving fitted values and residuals.

Figure
Least squares as orthogonal projection onto the column space
column space 𝒞(X) Y Ŷ = Xβ̂ ε̂ (residual) ε̂ ⟂ every vector in 𝒞(X)

The central theorem

TheoremLeast-squares estimation and the normal equations

Under the linear model as defined above, with XX an n×pn \times p matrix of full column rank p≤np \le n, the residual sum of squares

S(b)=∥Y−Xb∥2S(\boldsymbol{b}) = \|\boldsymbol{Y} - X\boldsymbol{b}\|^2

has a unique minimizer over b∈Rp\boldsymbol{b} \in \mathbb{R}^p, given by

β^=(X⊤X)−1X⊤Y.\hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{Y}.

Equivalently, β^\hat{\boldsymbol{\beta}} is the unique solution of the normal equations X⊤Xb=X⊤YX^\top X \boldsymbol{b} = X^\top \boldsymbol{Y}. Writing H=X(X⊤X)−1X⊤H = X(X^\top X)^{-1}X^\top (the hat matrix), the fitted values and residuals satisfy: 1. Y^=HY\hat{\boldsymbol{Y}} = H\boldsymbol{Y}, and HH is symmetric idempotent (H⊤=HH^\top = H, H2=HH^2 = H) with rank⁡(H)=p\operatorname{rank}(H) = p. 2. ε^=(I−H)Y\hat{\boldsymbol{\varepsilon}} = (I - H)\boldsymbol{Y}, and X⊤ε^=0X^\top \hat{\boldsymbol{\varepsilon}} = \boldsymbol{0}; in particular Y^⊤ε^=0\hat{\boldsymbol{Y}}^\top \hat{\boldsymbol{\varepsilon}} = 0 (fitted values and residuals are orthogonal). 3. (Pythagorean decomposition) ∥Y∥2=∥Y^∥2+∥ε^∥2\|\boldsymbol{Y}\|^2 = \|\hat{\boldsymbol{Y}}\|^2 + \|\hat{\boldsymbol{\varepsilon}}\|^2. 4. E[β^]=β\mathbb{E}[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta} (unbiasedness), and Cov⁡(β^)=σ2(X⊤X)−1\operatorname{Cov}(\hat{\boldsymbol{\beta}}) = \sigma^2 (X^\top X)^{-1}.

Proof

We give two independent derivations of the minimizer — algebraic (via calculus/completing the square) and geometric (via orthogonal projection) — then establish the stated properties.

Step 1: existence and uniqueness of the minimizer, algebraically. Let β^=(X⊤X)−1X⊤Y\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top \boldsymbol{Y}; this is well defined since XX has full column rank pp, which forces X⊤XX^\top X to be positive definite (for any v≠0\boldsymbol{v}\neq\boldsymbol{0}, v⊤X⊤Xv=∥Xv∥2>0\boldsymbol{v}^\top X^\top X \boldsymbol{v} = \|X\boldsymbol{v}\|^2 > 0 because Xv=0X\boldsymbol{v}=\boldsymbol0 only for v=0\boldsymbol v=\boldsymbol 0 by full rank), hence invertible. For an arbitrary b∈Rp\boldsymbol{b} \in \mathbb{R}^p, write b=β^+d\boldsymbol{b} = \hat{\boldsymbol{\beta}} + \boldsymbol{d} for d=b−β^\boldsymbol{d} = \boldsymbol{b} - \hat{\boldsymbol{\beta}}. Then

Y−Xb=(Y−Xβ^)−Xd=ε^−Xd,\boldsymbol{Y} - X\boldsymbol{b} = (\boldsymbol{Y} - X\hat{\boldsymbol{\beta}}) - X\boldsymbol{d} = \hat{\boldsymbol{\varepsilon}} - X\boldsymbol{d},

where ε^:=Y−Xβ^\hat{\boldsymbol{\varepsilon}} := \boldsymbol{Y} - X\hat{\boldsymbol{\beta}}. Expanding the squared norm,

S(b)=∥ε^∥2−2ε^⊤Xd+∥Xd∥2.S(\boldsymbol{b}) = \|\hat{\boldsymbol{\varepsilon}}\|^2 - 2\hat{\boldsymbol{\varepsilon}}^\top X\boldsymbol{d} + \|X\boldsymbol{d}\|^2.

We check the cross term vanishes: by definition of β^\hat{\boldsymbol{\beta}},

X⊤ε^=X⊤(Y−X(X⊤X)−1X⊤Y)=X⊤Y−X⊤Y=0,X^\top \hat{\boldsymbol{\varepsilon}} = X^\top\bigl(\boldsymbol{Y} - X(X^\top X)^{-1}X^\top \boldsymbol{Y}\bigr) = X^\top \boldsymbol{Y} - X^\top \boldsymbol{Y} = \boldsymbol{0},

so ε^⊤Xd=(X⊤ε^)⊤d=0\hat{\boldsymbol{\varepsilon}}^\top X \boldsymbol{d} = (X^\top\hat{\boldsymbol{\varepsilon}})^\top \boldsymbol{d} = 0 for every d\boldsymbol{d}. Hence

S(b)=∥ε^∥2+∥Xd∥2=S(β^)+∥X(b−β^)∥2≥S(β^),S(\boldsymbol{b}) = \|\hat{\boldsymbol{\varepsilon}}\|^2 + \|X\boldsymbol{d}\|^2 = S(\hat{\boldsymbol{\beta}}) + \|X(\boldsymbol{b}-\hat{\boldsymbol{\beta}})\|^2 \ge S(\hat{\boldsymbol{\beta}}),

with equality if and only if ∥X(b−β^)∥2=0\|X(\boldsymbol{b}-\hat{\boldsymbol{\beta}})\|^2 = 0, i.e. X(b−β^)=0X(\boldsymbol{b}-\hat{\boldsymbol{\beta}}) = \boldsymbol{0}. Since XX has full column rank, its columns are linearly independent, so Xv=0  ⟺  v=0X\boldsymbol{v}=\boldsymbol0 \iff \boldsymbol v = \boldsymbol 0; hence equality holds iff b=β^\boldsymbol{b} = \hat{\boldsymbol{\beta}}. This proves β^\hat{\boldsymbol{\beta}} is the unique global minimizer of SS.

Step 1′: the same minimizer via orthogonal projection (geometric verification). Let C(X)={Xb:b∈Rp}\mathcal{C}(X) = \{X\boldsymbol{b} : \boldsymbol{b} \in \mathbb{R}^p\}, a pp-dimensional subspace of Rn\mathbb{R}^n (dimension pp because XX's columns are independent). The projection theorem for Hilbert/Euclidean spaces states that for any subspace C\mathcal{C} and any Y\boldsymbol{Y}, there is a unique Y^∈C\hat{\boldsymbol{Y}} \in \mathcal{C} minimizing ∥Y−v∥\|\boldsymbol{Y} - \boldsymbol{v}\| over v∈C\boldsymbol{v} \in \mathcal{C}, and it is characterized by Y−Y^⊥C\boldsymbol{Y} - \hat{\boldsymbol{Y}} \perp \mathcal{C}. Writing Y^=Xβ^\hat{\boldsymbol{Y}} = X\hat{\boldsymbol{\beta}} for some β^\hat{\boldsymbol{\beta}} (possible since Y^∈C(X)\hat{\boldsymbol{Y}} \in \mathcal{C}(X)), the orthogonality condition Y−Y^⊥C(X)\boldsymbol{Y}-\hat{\boldsymbol{Y}} \perp \mathcal{C}(X) means (Y−Xβ^)⊤Xc=0(\boldsymbol{Y} - X\hat{\boldsymbol{\beta}})^\top X\boldsymbol{c} = 0 for all c∈Rp\boldsymbol{c} \in \mathbb{R}^p, i.e. X⊤(Y−Xβ^)=0X^\top(\boldsymbol{Y}-X\hat{\boldsymbol{\beta}}) = \boldsymbol{0}, i.e. X⊤Xβ^=X⊤YX^\top X \hat{\boldsymbol{\beta}} = X^\top \boldsymbol{Y} — the normal equations, solved uniquely (by invertibility of X⊤XX^\top X) by β^=(X⊤X)−1X⊤Y\hat{\boldsymbol{\beta}}=(X^\top X)^{-1}X^\top\boldsymbol{Y}, matching Step 1. This confirms the two routes agree and explains geometrically why X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon}=\boldsymbol 0 held in Step 1: it is precisely the projection theorem's orthogonality condition.

Step 2: properties of the hat matrix H=X(X⊤X)−1X⊤H = X(X^\top X)^{-1}X^\top. By definition Y^=Xβ^=X(X⊤X)−1X⊤Y=HY\hat{\boldsymbol{Y}} = X\hat{\boldsymbol{\beta}} = X(X^\top X)^{-1}X^\top \boldsymbol{Y} = H\boldsymbol{Y}, proving the first claim of part 1. Symmetry: H⊤=X((X⊤X)−1)⊤X⊤=X(X⊤X)−1X⊤=HH^\top = X\bigl((X^\top X)^{-1}\bigr)^\top X^\top = X(X^\top X)^{-1}X^\top = H, using that X⊤XX^\top X is symmetric so its inverse is too. Idempotence: H2=X(X⊤X)−1X⊤X(X⊤X)−1X⊤=X(X⊤X)−1X⊤=HH^2 = X(X^\top X)^{-1}X^\top X (X^\top X)^{-1}X^\top = X(X^\top X)^{-1}X^\top = H, since the middle X⊤X(X⊤X)−1X^\top X (X^\top X)^{-1} cancels to IpI_p. Rank: since HH is idempotent, its eigenvalues are all 00 or 11 (if Hv=λvH\boldsymbol v=\lambda\boldsymbol v with v≠0\boldsymbol v\ne 0 then λv=Hv=H2v=λHv=λ2v\lambda\boldsymbol v = H\boldsymbol v = H^2\boldsymbol v=\lambda H\boldsymbol v=\lambda^2\boldsymbol v, so λ=λ2\lambda=\lambda^2), and for a symmetric idempotent matrix, rank⁡(H)=trace⁡(H)\operatorname{rank}(H)=\operatorname{trace}(H) (sum of eigenvalues counted with multiplicity, all 00 or 11). Compute

trace⁡(H)=trace⁡(X(X⊤X)−1X⊤)=trace⁡((X⊤X)−1X⊤X)=trace⁡(Ip)=p,\operatorname{trace}(H) = \operatorname{trace}\bigl(X(X^\top X)^{-1}X^\top\bigr) = \operatorname{trace}\bigl((X^\top X)^{-1}X^\top X\bigr) = \operatorname{trace}(I_p) = p,

using the cyclic property of trace. Hence rank⁡(H)=p\operatorname{rank}(H) = p, proving part 1 completely.

Step 3: residuals and orthogonality (part 2). ε^=Y−Y^=Y−HY=(I−H)Y\hat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \hat{\boldsymbol{Y}} = \boldsymbol{Y} - H\boldsymbol{Y} = (I-H)\boldsymbol{Y}. We showed in Step 1 that X⊤ε^=0X^\top \hat{\boldsymbol\varepsilon} = \boldsymbol 0 directly from the normal equations; this also follows from I−HI - H: X⊤(I−H)=X⊤−X⊤X(X⊤X)−1X⊤=X⊤−X⊤=0X^\top(I-H) = X^\top - X^\top X (X^\top X)^{-1}X^\top = X^\top - X^\top = 0. Consequently, since Y^=Xβ^∈C(X)\hat{\boldsymbol{Y}} = X\hat{\boldsymbol{\beta}} \in \mathcal{C}(X),

Y^⊤ε^=β^⊤X⊤ε^=β^⊤0=0.\hat{\boldsymbol{Y}}^\top \hat{\boldsymbol{\varepsilon}} = \hat{\boldsymbol{\beta}}^\top X^\top \hat{\boldsymbol{\varepsilon}} = \hat{\boldsymbol{\beta}}^\top \boldsymbol{0} = 0.

This proves part 2.

Step 4: Pythagorean decomposition (part 3). Since Y=Y^+ε^\boldsymbol Y = \hat{\boldsymbol Y}+\hat{\boldsymbol\varepsilon} and Y^⊤ε^=0\hat{\boldsymbol Y}^\top\hat{\boldsymbol\varepsilon}=0,

∥Y∥2=∥Y^+ε^∥2=∥Y^∥2+2Y^⊤ε^+∥ε^∥2=∥Y^∥2+∥ε^∥2.\|\boldsymbol{Y}\|^2 = \|\hat{\boldsymbol{Y}} + \hat{\boldsymbol{\varepsilon}}\|^2 = \|\hat{\boldsymbol{Y}}\|^2 + 2\hat{\boldsymbol{Y}}^\top\hat{\boldsymbol{\varepsilon}} + \|\hat{\boldsymbol{\varepsilon}}\|^2 = \|\hat{\boldsymbol{Y}}\|^2 + \|\hat{\boldsymbol{\varepsilon}}\|^2.

Step 5: unbiasedness and covariance (part 4). Substitute the model Y=Xβ+ε\boldsymbol{Y} = X\boldsymbol{\beta}+\boldsymbol{\varepsilon} into the estimator:

β^=(X⊤X)−1X⊤(Xβ+ε)=(X⊤X)−1X⊤Xβ+(X⊤X)−1X⊤ε=β+(X⊤X)−1X⊤ε.\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top(X\boldsymbol\beta+\boldsymbol\varepsilon) = (X^\top X)^{-1}X^\top X \boldsymbol\beta + (X^\top X)^{-1}X^\top\boldsymbol\varepsilon = \boldsymbol\beta + (X^\top X)^{-1}X^\top \boldsymbol\varepsilon.

Here XX is treated as fixed (non-random), the standard setup for the linear model. Taking expectations, using E[ε]=0\mathbb{E}[\boldsymbol\varepsilon]=\boldsymbol 0,

E[β^]=β+(X⊤X)−1X⊤E[ε]=β.\mathbb{E}[\hat{\boldsymbol\beta}] = \boldsymbol\beta + (X^\top X)^{-1}X^\top \mathbb{E}[\boldsymbol\varepsilon] = \boldsymbol\beta.

For the covariance, let A=(X⊤X)−1X⊤A = (X^\top X)^{-1}X^\top, so β^−β=Aε\hat{\boldsymbol\beta}-\boldsymbol\beta = A\boldsymbol\varepsilon. Then, using Cov⁡(ε)=σ2In\operatorname{Cov}(\boldsymbol\varepsilon)=\sigma^2 I_n,

Cov⁡(β^)=Cov⁡(Aε)=ACov⁡(ε)A⊤=σ2AA⊤=σ2(X⊤X)−1X⊤X(X⊤X)−1=σ2(X⊤X)−1,\operatorname{Cov}(\hat{\boldsymbol\beta}) = \operatorname{Cov}(A\boldsymbol\varepsilon) = A\operatorname{Cov}(\boldsymbol\varepsilon)A^\top = \sigma^2 A A^\top = \sigma^2 (X^\top X)^{-1}X^\top X (X^\top X)^{-1} = \sigma^2(X^\top X)^{-1},

where we used ((X⊤X)−1)⊤=(X⊤X)−1((X^\top X)^{-1})^\top = (X^\top X)^{-1} (symmetry of the inverse of a symmetric matrix). This proves part 4 and completes the proof. ■\blacksquare

∎

Worked examples

Example 1: simple linear regression with unequal, asymmetric data. Consider n=5n=5 observations of a single predictor xx (say, hours of machine operation) and response yy (units produced, with measurement noise):

iixix_iyiy_i
112.1
223.9
346.3
4511.2
5815.5

Here p=2p=2 with design matrix rows xi⊤=(1,xi)\boldsymbol{x}_i^\top = (1, x_i) (intercept plus slope). We need

X⊤X=(n∑xi∑xi∑xi2),X⊤Y=(∑yi∑xiyi).X^\top X = \begin{pmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{pmatrix}, \qquad X^\top \boldsymbol{Y} = \begin{pmatrix} \sum y_i \\ \sum x_i y_i \end{pmatrix}.

Compute the sums: ∑xi=1+2+4+5+8=20\sum x_i = 1+2+4+5+8 = 20, ∑xi2=1+4+16+25+64=110\sum x_i^2 = 1+4+16+25+64=110, ∑yi=2.1+3.9+6.3+11.2+15.5=39.0\sum y_i = 2.1+3.9+6.3+11.2+15.5 = 39.0, and

∑xiyi=1(2.1)+2(3.9)+4(6.3)+5(11.2)+8(15.5)=2.1+7.8+25.2+56.0+124.0=215.1.\sum x_i y_i = 1(2.1)+2(3.9)+4(6.3)+5(11.2)+8(15.5) = 2.1+7.8+25.2+56.0+124.0 = 215.1.

So

X⊤X=(52020110),X⊤Y=(39.0215.1).X^\top X = \begin{pmatrix} 5 & 20 \\ 20 & 110 \end{pmatrix}, \qquad X^\top \boldsymbol{Y} = \begin{pmatrix} 39.0 \\ 215.1 \end{pmatrix}.

The determinant is 5(110)−20(20)=550−400=1505(110) - 20(20) = 550 - 400 = 150, so

(X⊤X)−1=1150(110−20−205).(X^\top X)^{-1} = \frac{1}{150}\begin{pmatrix} 110 & -20 \\ -20 & 5 \end{pmatrix}.

Then

β^=(X⊤X)−1X⊤Y=1150(110(39.0)+(−20)(215.1)(−20)(39.0)+5(215.1))=1150(4290−4302−780+1075.5)=1150(−12295.5).\hat{\boldsymbol\beta} = (X^\top X)^{-1} X^\top \boldsymbol Y = \frac{1}{150}\begin{pmatrix} 110(39.0) + (-20)(215.1) \\ (-20)(39.0) + 5(215.1) \end{pmatrix} = \frac{1}{150}\begin{pmatrix} 4290 - 4302 \\ -780 + 1075.5 \end{pmatrix} = \frac{1}{150}\begin{pmatrix} -12 \\ 295.5 \end{pmatrix}.

Hence β^1=−12/150=−0.08\hat\beta_1 = -12/150 = -0.08 (intercept) and β^2=295.5/150=1.97\hat\beta_2 = 295.5/150 = 1.97 (slope). The fitted line is

y^=−0.08+1.97x.\hat{y} = -0.08 + 1.97x.

Check the normal equation residual-orthogonality identity numerically. Fitted values:

y^1=−0.08+1.97(1)=1.89,y^2=−0.08+1.97(2)=3.86,y^3=−0.08+1.97(4)=7.80,\hat y_1 = -0.08+1.97(1) = 1.89, \quad \hat y_2 = -0.08+1.97(2)=3.86, \quad \hat y_3 = -0.08+1.97(4)=7.80,
y^4=−0.08+1.97(5)=9.77,y^5=−0.08+1.97(8)=15.68.\hat y_4 = -0.08+1.97(5)=9.77, \quad \hat y_5 = -0.08+1.97(8)=15.68.

Residuals: ε^1=2.1−1.89=0.21\hat\varepsilon_1 = 2.1-1.89=0.21, ε^2=3.9−3.86=0.04\hat\varepsilon_2=3.9-3.86=0.04, ε^3=6.3−7.80=−1.50\hat\varepsilon_3=6.3-7.80=-1.50, ε^4=11.2−9.77=1.43\hat\varepsilon_4=11.2-9.77=1.43, ε^5=15.5−15.68=−0.18\hat\varepsilon_5=15.5-15.68=-0.18. Verify ∑ε^i≈0\sum \hat\varepsilon_i \approx 0: 0.21+0.04−1.50+1.43−0.18=0.000.21+0.04-1.50+1.43-0.18 = 0.00 ✓ (this must hold exactly since the intercept column is 11's and X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon}=\boldsymbol 0 forces ∑i1⋅ε^i=0\sum_i 1\cdot\hat\varepsilon_i = 0). Verify ∑xiε^i≈0\sum x_i \hat\varepsilon_i \approx 0: 1(0.21)+2(0.04)+4(−1.50)+5(1.43)+8(−0.18)=0.21+0.08−6.00+7.15−1.44=0.001(0.21)+2(0.04)+4(-1.50)+5(1.43)+8(-0.18) = 0.21+0.08-6.00+7.15-1.44 = 0.00 ✓, confirming the second column of the orthogonality identity X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon}=\boldsymbol 0.

Example 2: multiple regression with two predictors, illustrating the projection picture directly. Suppose we model a response YY (say, plant growth in cm) using two predictors: x1x_1 = water (liters/week, no intercept term included in this design for simplicity) and x2x_2 = sunlight (hours/day), with n=4n=4 observations:

iixi1x_{i1} (water)xi2x_{i2} (sun)yiy_i (growth)
1135.0
2226.0
33514.0
4417.5

Model: Yi=β1xi1+β2xi2+εiY_i = \beta_1 x_{i1} + \beta_2 x_{i2} + \varepsilon_i (no intercept — the two-column design matrix directly).

X=(13223541),Y=(5.06.014.07.5).X = \begin{pmatrix} 1 & 3 \\ 2 & 2 \\ 3 & 5 \\ 4 & 1 \end{pmatrix}, \qquad \boldsymbol Y = \begin{pmatrix} 5.0\\6.0\\14.0\\7.5 \end{pmatrix}.

Compute X⊤XX^\top X: entry (1,1)=12+22+32+42=1+4+9+16=30(1,1) = 1^2+2^2+3^2+4^2 = 1+4+9+16=30; entry (1,2)=(2,1)=1(3)+2(2)+3(5)+4(1)=3+4+15+4=26(1,2)=(2,1) = 1(3)+2(2)+3(5)+4(1) = 3+4+15+4=26; entry (2,2)=32+22+52+12=9+4+25+1=39(2,2) = 3^2+2^2+5^2+1^2 = 9+4+25+1=39. So

X⊤X=(30262639).X^\top X = \begin{pmatrix} 30 & 26 \\ 26 & 39 \end{pmatrix}.

Compute X⊤YX^\top \boldsymbol Y: entry 1=1(5.0)+2(6.0)+3(14.0)+4(7.5)=5+12+42+30=891 = 1(5.0)+2(6.0)+3(14.0)+4(7.5) = 5+12+42+30=89; entry 2=3(5.0)+2(6.0)+5(14.0)+1(7.5)=15+12+70+7.5=104.52 = 3(5.0)+2(6.0)+5(14.0)+1(7.5) = 15+12+70+7.5=104.5. So

X⊤Y=(89104.5).X^\top \boldsymbol Y = \begin{pmatrix} 89 \\ 104.5 \end{pmatrix}.

Determinant of X⊤XX^\top X: 30(39)−26(26)=1170−676=49430(39)-26(26) = 1170 - 676 = 494. Inverse:

(X⊤X)−1=1494(39−26−2630).(X^\top X)^{-1} = \frac{1}{494}\begin{pmatrix} 39 & -26 \\ -26 & 30 \end{pmatrix}.

Then

β^=1494(39(89)+(−26)(104.5)(−26)(89)+30(104.5))=1494(3471−2717−2314+3135)=1494(754821).\hat{\boldsymbol\beta} = \frac{1}{494}\begin{pmatrix} 39(89) + (-26)(104.5) \\ (-26)(89) + 30(104.5) \end{pmatrix} = \frac{1}{494}\begin{pmatrix} 3471 - 2717 \\ -2314+3135 \end{pmatrix} = \frac{1}{494}\begin{pmatrix} 754 \\ 821 \end{pmatrix}.

So β^1=754/494≈1.5263\hat\beta_1 = 754/494 \approx 1.5263 and β^2=821/494≈1.6620\hat\beta_2 = 821/494 \approx 1.6620. Fitted values: y^1=1.5263(1)+1.6620(3)=1.5263+4.9859=6.5123\hat y_1 = 1.5263(1)+1.6620(3) = 1.5263+4.9859 = 6.5123; y^2=1.5263(2)+1.6620(2)=3.0526+3.3239=6.3765\hat y_2 = 1.5263(2)+1.6620(2)=3.0526+3.3239=6.3765; y^3=1.5263(3)+1.6620(5)=4.5789+8.3097=12.8886\hat y_3 = 1.5263(3)+1.6620(5)=4.5789+8.3097=12.8886; y^4=1.5263(4)+1.6620(1)=6.1053+1.6620=7.7673\hat y_4=1.5263(4)+1.6620(1)=6.1053+1.6620=7.7673. Residuals: ε^1=5.0−6.5123=−1.5123\hat\varepsilon_1 = 5.0-6.5123=-1.5123, ε^2=6.0−6.3765=−0.3765\hat\varepsilon_2=6.0-6.3765=-0.3765, ε^3=14.0−12.8886=1.1114\hat\varepsilon_3=14.0-12.8886=1.1114, ε^4=7.5−7.7673=−0.2673\hat\varepsilon_4=7.5-7.7673=-0.2673. Check x⋅1⊤ε^≈0\boldsymbol x_{\cdot 1}^\top \hat{\boldsymbol\varepsilon} \approx 0: 1(−1.5123)+2(−0.3765)+3(1.1114)+4(−0.2673)=−1.5123−0.7530+3.3342−1.0692=−0.0003≈01(-1.5123)+2(-0.3765)+3(1.1114)+4(-0.2673) = -1.5123-0.7530+3.3342-1.0692 = -0.0003 \approx 0 ✓ (rounding). Check x⋅2⊤ε^≈0\boldsymbol x_{\cdot 2}^\top\hat{\boldsymbol\varepsilon}\approx 0: 3(−1.5123)+2(−0.3765)+5(1.1114)+1(−0.2673)=−4.5369−0.7530+5.5570−0.2673=−0.0002≈03(-1.5123)+2(-0.3765)+5(1.1114)+1(-0.2673) = -4.5369-0.7530+5.5570-0.2673 = -0.0002 \approx 0 ✓. Both confirm X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon} = \boldsymbol 0 to rounding precision, exactly as guaranteed by the theorem — note this design has no intercept column, so unlike Example 1 the raw residuals need not sum to zero on their own; only the orthogonality to each actual column of XX is guaranteed, illustrating that "residuals sum to zero" is a special case of X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon}=0 that requires a column of ones, not a universal law.

CounterexampleDropping full column rank: non-uniqueness of the minimizer

The theorem's uniqueness claim relies on XX having full column rank pp. Suppose instead we regress on two predictors that are perfectly collinear. Take n=3n=3, p=2p=2, with

X=(122436),Y=(134),X = \begin{pmatrix} 1 & 2 \\ 2 & 4 \\ 3 & 6 \end{pmatrix}, \qquad \boldsymbol Y = \begin{pmatrix} 1 \\ 3 \\ 4 \end{pmatrix},

where the second column is exactly twice the first (xi2=2xi1x_{i2} = 2x_{i1} for every ii), so XX has rank 1<p=21 < p = 2. Then X⊤X=(14282856)X^\top X = \begin{pmatrix} 14 & 28 \\ 28 & 56\end{pmatrix} has determinant 14(56)−28(28)=784−784=014(56)-28(28) = 784-784=0, so it is singular and the formula β^=(X⊤X)−1X⊤Y\hat{\boldsymbol\beta}=(X^\top X)^{-1}X^\top \boldsymbol Y is not defined. Worse, the minimizer of S(b)S(\boldsymbol b) is no longer unique. Any b=(b1,b2)\boldsymbol b=(b_1,b_2) produces the same fitted value Xb=(b1+2b2)(123)X\boldsymbol b = (b_1+2b_2)\begin{pmatrix}1\\2\\3\end{pmatrix}, which depends on (b1,b2)(b_1,b_2) only through the scalar c:=b1+2b2c := b_1+2b_2. Minimizing over cc gives an ordinary one-variable least-squares problem with predictor column (1,2,3)⊤(1,2,3)^\top; one finds the optimal c∗c^* is unique, but every pair (b1,b2)(b_1,b_2) with b1+2b2=c∗b_1+2b_2=c^* achieves the same minimal sum of squares. For instance (b1,b2)=(c∗,0)(b_1,b_2)=(c^*,0) and (b1,b2)=(0,c∗/2)(b_1,b_2)=(0,c^*/2) both minimize SS identically. There is an entire line of minimizers in R2\mathbb{R}^2, not a single point. This shows precisely what full column rank buys us: without it, X⊤XX^\top X is not invertible, the normal equations have infinitely many solutions (a whole affine subspace), and "the" least-squares estimator is not well-defined as a single vector — one must either drop a collinear predictor or use a generalized (Moore–Penrose) inverse and select, by convention, the minimum-norm solution.

Scope, edge cases & common traps

SituationWhat happensResolution
X⊤XX^\top X singular (collinear predictors, or n<pn<p)Normal equations have infinitely many solutionsDrop redundant predictors, or use pseudoinverse for minimum-norm solution
No intercept column in XX∑iε^i\sum_i \hat\varepsilon_i need not be 00Only X⊤ε^=0X^\top\hat{\boldsymbol\varepsilon}=\boldsymbol 0 is guaranteed in general
Errors correlated or heteroscedastic (Cov⁡(ε)≠σ2I\operatorname{Cov}(\boldsymbol\varepsilon)\ne \sigma^2 I)β^\hat{\boldsymbol\beta} still unbiased, but the variance formula σ2(X⊤X)−1\sigma^2(X^\top X)^{-1} is wrongUse generalized least squares; see foreshadowed The Gauss-Markov Theorem discussion
XX random rather than fixedUnbiasedness proof needs conditioning on XXState E[β^∣X]=β\mathbb{E}[\hat{\boldsymbol\beta}\mid X]=\boldsymbol\beta
Want inference (CIs, tt-tests on β^j\hat\beta_j)Requires a distributional assumption on εi\varepsilon_iTypically assume εi∼N(0,σ2)\varepsilon_i \sim N(0,\sigma^2) i.i.d., invoking Standard Tests: t, Chi-Squared, and F

The mathematician behind it

Origins · 1 minCarl Friedrich Gauss and how this was discovered
CF
Carl Friedrich Gauss
1777–1855 · Göttingen, German astronomy and mathematics

Gauss developed the method of least squares around 1795, motivated by the problem of combining redundant astronomical observations — most famously, predicting the orbit of the newly discovered dwarf planet Ceres in 1801 after it vanished behind the sun, using only a handful of noisy angular measurements. He justified the method by showing that, under a normal (Gaussian) error law, the least-squares estimate coincides with the maximum-likelihood estimate — a foundational link between geometry and probability that anticipates modern statistical inference. Adrien-Marie Legendre independently published the method first, in 1805, leading to a long priority dispute; Gauss maintained he had used the technique privately since the 1790s, and his 1809 book Theoria Motus Corporum Coelestium gives the probabilistic justification absent from Legendre's purely algebraic treatment.

History
- 1795 — Gauss privately develops the method of least squares for combining astronomical measurements.
- 1801 — Gauss successfully predicts the reappearance of Ceres, dramatically validating the method.
- 1805 — Legendre independently publishes the first formal description of least squares.
- 1809 — Gauss publishes Theoria Motus, linking least squares to maximum likelihood under Gaussian errors.
- 1821 — Gauss publishes the theorem now called Gauss–Markov, establishing an optimality property of the estimator.
Why minimize the sum of squares rather than, say, the sum of absolute residuals?

Squaring gives a smooth, differentiable objective with a unique closed-form global minimizer whenever XX has full column rank, obtained by solving a linear system (the normal equations). Minimizing absolute deviations (least absolute deviations, or ℓ1\ell_1 regression) is a valid alternative but leads to a piecewise-linear, non-differentiable objective solved by linear programming rather than linear algebra, and it targets the conditional median rather than the conditional mean. Historically, Gauss further justified the squared criterion by showing it coincides with maximum likelihood when errors are Gaussian — see the discussion in The Sample Mean and Sample Variance regarding why squared deviations arise naturally from Gaussian likelihoods. Both criteria are legitimate; least squares is simply the one with the richest closed-form theory.

Does the least-squares theorem require the errors to be normally distributed?

No. The theorem proved here — the formula for β^\hat{\boldsymbol\beta}, its unbiasedness, and its covariance matrix σ2(X⊤X)−1\sigma^2(X^\top X)^{-1} — uses only the assumptions that E[εi]=0\mathbb{E}[\varepsilon_i]=0, Var⁡(εi)=σ2\operatorname{Var}(\varepsilon_i)=\sigma^2, and the εi\varepsilon_i are uncorrelated. Normality is a separate, additional assumption needed only when you want exact finite-sample distributions for β^\hat{\boldsymbol\beta}, for confidence intervals, or for tt- and FF-tests on the coefficients, as covered under Standard Tests: t, Chi-Squared, and F. Without normality, one typically appeals instead to large-sample (central limit theorem) approximations for inference.

What exactly is the "hat matrix" and why is it useful?

The hat matrix H=X(X⊤X)−1X⊤H = X(X^\top X)^{-1}X^\top is the matrix that maps the observed response vector to the fitted values, Y^=HY\hat{\boldsymbol Y}=H\boldsymbol Y — it "puts a hat on" Y\boldsymbol Y. It is symmetric and idempotent, which encodes precisely the fact that it represents an orthogonal projection onto the column space of XX. Its diagonal entries, called leverages, measure how much each observation's own response value influences its own fitted value, and they sum to pp (the rank of HH), which is useful for diagnosing influential data points in regression diagnostics.

Is the least-squares estimator the best possible unbiased estimator of β\boldsymbol\beta?

Not automatically, and this article deliberately does not claim it. What is proved here is only that β^\hat{\boldsymbol\beta} is unbiased with covariance σ2(X⊤X)−1\sigma^2(X^\top X)^{-1}, under the stated linear model with uncorrelated, equal-variance errors. Whether it has the smallest variance among all linear unbiased estimators is a distinct and stronger claim, established by The Gauss-Markov Theorem, which requires exactly the same error assumptions (mean zero, common variance, uncorrelated) but no normality. If those assumptions fail — for instance under heteroscedastic or correlated errors — least squares can be badly suboptimal, and generalized least squares should be used instead.

What happens numerically if X⊤XX^\top X is invertible but nearly singular?

The formula β^=(X⊤X)−1X⊤Y\hat{\boldsymbol\beta}=(X^\top X)^{-1}X^\top\boldsymbol Y remains mathematically valid, but it becomes numerically unstable: small perturbations in the data can produce large changes in β^\hat{\boldsymbol\beta}, and the entries of (X⊤X)−1(X^\top X)^{-1} can be enormous, inflating Cov⁡(β^)=σ2(X⊤X)−1\operatorname{Cov}(\hat{\boldsymbol\beta})=\sigma^2(X^\top X)^{-1}. This situation, called multicollinearity, does not violate any hypothesis of the theorem, but it is a serious practical warning sign; standard remedies include centering/scaling predictors, dropping redundant variables, or using ridge regression, which regularizes X⊤XX^\top X by adding λI\lambda I before inverting.

How does the fitted-values/residuals decomposition relate to hypothesis testing?

The Pythagorean identity ∥Y∥2=∥Y^∥2+∥ε^∥2\|\boldsymbol Y\|^2 = \|\hat{\boldsymbol Y}\|^2+\|\hat{\boldsymbol\varepsilon}\|^2 (or its mean-centered analogue, the total-sum-of-squares decomposition) is the algebraic seed of analysis of variance. Under the additional assumption εi∼N(0,σ2)\varepsilon_i \sim N(0,\sigma^2) i.i.d., ∥ε^∥2/σ2\|\hat{\boldsymbol\varepsilon}\|^2/\sigma^2 has a chi-squared distribution with n−pn-p degrees of freedom, and ratios of such sums of squares (for nested models) have FF-distributions — the basis of the FF-test for comparing regression models, and of individual coefficient tt-tests, both covered under Standard Tests: t, Chi-Squared, and F.

Can least squares be used with more predictors than observations (p>np>n)?

Not in the form given here: if p>np>n, the design matrix cannot have full column rank pp (its rank is at most n<pn<p), so X⊤XX^\top X is necessarily singular and the theorem's hypotheses fail outright. In that regime, ordinary least squares has no unique solution, and practitioners instead use regularized methods such as ridge regression or the lasso, which impose additional structure (a penalty term) to select a unique, well-behaved estimator despite p>np>n.

Before this

Where this leads