Least Squares & Normal Equations — Least squares estimation...

SECTION 1

Historical Context & Motivation

The problem of fitting a line—or more generally, a linear model—to observed data is among the oldest and most consequential in the mathematical sciences. Long before the formal apparatus of statistical inference existed, astronomers and surveyors grappled with the practical reality that measurements are noisy: multiple observations of the same celestial position or geographic angle yield slightly different values. The method of least squares arose as a principled answer to the question of how to reconcile conflicting observations into a single best estimate. Its development is intertwined with the birth of modern statistics, and its mathematical elegance has ensured its place at the center of regression analysis, signal processing, and machine learning to this day.

1795

Gauss's Unpublished Work

Carl Friedrich Gauss, at just 18 years old, develops the method of least squares to predict the orbit of Ceres, though he does not publish the technique until 1809 in Theoria Motus Corporum Coelestium. His derivation connects least squares to the assumption of normally distributed errors.

1805

Legendre Publishes First

Adrien-Marie Legendre formally publishes the method of least squares in Nouvelles méthodes pour la détermination des orbites des comètes, providing clear algebraic formulations. This marks the earliest publicly available description of the technique.

1809

Gauss and the Normal Distribution

Gauss publishes his derivation linking the least squares criterion to the maximum likelihood principle under Gaussian errors, establishing a deep probabilistic foundation for the method.

1822

Gauss–Markov Theorem Foundation

Gauss proves that among all linear unbiased estimators, the ordinary least squares estimator has the smallest variance—a result later formalized and extended by Andrey Markov. The Gauss–Markov theorem becomes a cornerstone of linear model theory.

1900s–present

Modern Extensions

The 20th century sees least squares generalized to weighted, generalized, and nonlinear forms. Computational advances—from Gaussian elimination to QR factorization and SVD—make large-scale least squares feasible, underpinning modern data science.

The central question that motivates this entire topic is deceptively simple: given a set of observations and a linear model that cannot perfectly reproduce them all, how should we choose the model parameters to obtain the "best" approximation? The least squares criterion offers a compelling answer by minimizing the sum of squared deviations between observed and predicted values, and the normal equations provide the algebraic machinery that translates this optimization into a system of linear equations with a unique closed-form solution.

SECTION 2

Core Principles & Definitions

Before diving into the algebra, it is essential to establish the conceptual pillars upon which least squares estimation rests. The framework begins with a linear model that relates an n × 1 response vector y to an n × p design matrix X and a p × 1 parameter vector β through the equation y = Xβ + ε, where ε captures random error. The overarching goal is to find the estimate β̂ that makes the model's predictions as close to the observed data as possible in the Euclidean sense.

The Linear Model

We posit y = Xβ + ε where y ∈ ℝⁿ, X ∈ ℝⁿˣᵖ (known), β ∈ ℝᵖ (unknown), and ε ∈ ℝⁿ (random error with E[ε] = 0). The column space of X defines the set of all achievable fitted values.

The Least Squares Criterion

Choose β̂ to minimize the residual sum of squares S(β) = (y − Xβ)ᵀ(y − Xβ) = ‖y − Xβ‖². This is a convex quadratic in β, guaranteeing a global minimum.

The Normal Equations

Setting the gradient ∂S/∂β = 0 yields XᵀXβ̂ = Xᵀy. When XᵀX is invertible (i.e., X has full column rank), the unique solution is β̂ = (XᵀX)⁻¹Xᵀy.

Geometric Interpretation

The fitted vector ŷ = Xβ̂ is the orthogonal projection of y onto the column space C(X). The residual vector e = y − ŷ is perpendicular to every column of X, which is precisely the content of the normal equations.

Gauss–Markov Optimality

Under E[ε] = 0 and Cov(ε) = σ²Iₙ, the OLS estimator β̂ is BLUE — the Best Linear Unbiased Estimator — meaning no other linear unbiased estimator has smaller variance for any estimable linear function cᵀβ.

✦ KEY TAKEAWAY

Think of the column space of X as a flat screen in high-dimensional space. The observed vector y is a point floating off the screen. The least squares solution finds the shadow of y on the screen—the point on the screen closest to y. The residual is the light beam connecting y to its shadow, and the normal equations encode the fact that this beam strikes the screen at a right angle. Just as a shadow is the best 2-D approximation of a 3-D object on a wall, ŷ is the best approximation of y within the subspace defined by the predictors.

SECTION 3

Geometric Visualization of Least Squares Projection

The deepest insight into least squares estimation comes from viewing it as an orthogonal projection. The diagram below represents this geometric picture in three dimensions. The column space C(X) is depicted as a two-dimensional plane (representing a p-dimensional subspace in general). The response vector y lives in ℝⁿ and typically does not lie on this plane. The OLS fitted value ŷ = Xβ̂ is the foot of the perpendicular from y to the plane, and the residual vector e = y − ŷ is orthogonal to every vector in C(X). This orthogonality condition is exactly what the normal equations XᵀXβ̂ = Xᵀy encode.

The plane represents the column space C(X) spanned by the predictor columns x₁ and x₂. The response vector y (pink) lies off the plane. The fitted vector ŷ (green) is the orthogonal projection of y onto the plane, and the dashed red line is the residual vector e, meeting the plane at a right angle.

This geometric perspective clarifies several facts simultaneously. First, because ŷ is a linear combination of the columns of X, it lives in C(X) by definition. Second, the orthogonality condition Xᵀ(y − Xβ̂) = 0 guarantees that no further movement within C(X) can reduce ‖y − Xβ‖². Third, the hat matrix H = X(XᵀX)⁻¹Xᵀ is the orthogonal projection matrix onto C(X), satisfying H² = H, Hᵀ = H, and ŷ = Hy. The Pythagorean theorem in ℝⁿ gives the fundamental decomposition ‖y‖² = ‖ŷ‖² + ‖e‖², which underlies the analysis of variance (ANOVA) decomposition in regression.

SECTION 4

Mathematical Framework — Derivation of the Normal Equations

We now derive the normal equations rigorously from the calculus-based minimization of the residual sum of squares. Consider the model y = Xβ + ε where X has full column rank p. The objective function is a scalar-valued quadratic form in the parameter vector β.

OBJECTIVE FUNCTION

S(β) = (y − Xβ)ᵀ(y − Xβ) = yᵀy − 2βᵀXᵀy + βᵀXᵀXβ

S(β) is the residual sum of squares (RSS). The expansion uses the symmetry of the inner product and the fact that yᵀXβ is a scalar equal to its transpose βᵀXᵀy.

Differentiating S(β) with respect to the vector β using standard matrix calculus identities, we obtain the gradient. Recall that ∂(βᵀa)/∂β = a and ∂(βᵀAβ)/∂β = 2Aβ when A is symmetric. Setting the gradient to the zero vector yields the normal equations.

GRADIENT AND FIRST-ORDER CONDITION

∂S/∂β = −2Xᵀy + 2XᵀXβ = 0 ⟹ XᵀXβ̂ = Xᵀy

The p × p matrix XᵀX is called the Gram matrix of the columns of X. The p × 1 vector Xᵀy is the cross-product of predictors with the response.

CLOSED-FORM OLS ESTIMATOR

β̂ = (XᵀX)⁻¹Xᵀy

This expression requires that XᵀX be nonsingular, equivalently that rank(X) = p (full column rank). The matrix (XᵀX)⁻¹Xᵀ is the Moore–Penrose left inverse of X when n > p.

To confirm that this critical point is indeed a minimum, we examine the Hessian matrix ∂²S/∂β∂βᵀ = 2XᵀX. Since X has full column rank, XᵀX is positive definite, which means S(β) is strictly convex, and the critical point β̂ is the unique global minimizer. This is a powerful result: there is no concern about local minima or saddle points in ordinary least squares.

COVARIANCE OF THE OLS ESTIMATOR

Cov(β̂) = σ²(XᵀX)⁻¹

Under the assumption Cov(ε) = σ²Iₙ, the covariance matrix of β̂ follows directly from Cov(Ay) = A · Cov(y) · Aᵀ with A = (XᵀX)⁻¹Xᵀ. The scalar σ² is typically estimated by s² = eᵀe/(n − p).

📐 Why "Normal" Equations?

The term "normal" does not refer to the Gaussian (normal) distribution. It comes from the geometric fact that the residual vector e = y − Xβ̂ is normal (perpendicular) to the column space of X. The condition Xᵀe = 0 states that e is orthogonal to every column of X, hence "normal" in the geometric sense.

SECTION 5

Detailed Breakdown — The Hat Matrix and ANOVA Decomposition

Once the OLS estimator β̂ = (XᵀX)⁻¹Xᵀy is obtained, several important derived quantities illuminate the structure of the fit. The hat matrix H = X(XᵀX)⁻¹Xᵀ maps the observed response y to the fitted values ŷ = Hy. It is idempotent (H² = H), symmetric (Hᵀ = H), and has trace equal to p, the number of parameters. Similarly, the residual-forming matrix I − H projects y onto the orthogonal complement of C(X), producing the residuals e = (I − H)y. The trace of I − H is n − p, which explains why the unbiased estimator of σ² divides by n − p rather than n.

The total sum of squares SST decomposes into the regression sum of squares SSR and the residual sum of squares SSE. The coefficient of determination R² measures the fraction of total variability captured by the model.

This decomposition has a beautiful Pythagorean structure. The vectors ŷ − ȳ1 and e lie in orthogonal subspaces, so their squared norms add. The coefficient of determination R² = SSR/SST is the cosine squared of the angle between y − ȳ1 and ŷ − ȳ1 in ℝⁿ. When R² = 1, y lies entirely in the column space (perfect fit); when R² = 0, the fitted values provide no improvement over the grand mean ȳ.

Key quantities derived from the OLS solution

Quantity	Formula	Properties
Hat Matrix H	`X(XᵀX)⁻¹Xᵀ`	Symmetric, idempotent, tr(H) = p, eigenvalues 0 or 1
Fitted Values ŷ	`Hy = Xβ̂`	Orthogonal projection of y onto C(X)
Residuals e	`(I − H)y`	Xᵀe = 0; E[e] = 0; Cov(e) = σ²(I − H)
σ² estimator s²	`eᵀe / (n − p)`	Unbiased: E[s²] = σ²; denominator from tr(I − H)

SECTION 6

Worked Example — Simple Linear Regression via Normal Equations

Consider a simple linear regression with n = 4 observations. We model a response y as a function of a single predictor x: yᵢ = β₀ + β₁xᵢ + εᵢ. The data are (x, y) = {(1, 2), (2, 3), (3, 5), (4, 4)}. We will construct the design matrix, form the normal equations, solve for β̂, and compute the residuals.

Solving the Normal Equations for Simple Linear Regression

Step 1 — Construct the Design Matrix X and Response Vector yThe model yᵢ = β₀ + β₁xᵢ can be written in matrix form as y = Xβ where the first column of X is a column of ones (for the intercept) and the second column contains the predictor values. Thus: X = [[1, 1], [1, 2], [1, 3], [1, 4]] and y = [2, 3, 5, 4]ᵀ

X is 4 × 2 with full column rank 2; β = [β₀, β₁]ᵀ

Step 2 — Compute XᵀX and XᵀyXᵀX = [[∑1, ∑xᵢ], [∑xᵢ, ∑xᵢ²]] = [[4, 10], [10, 30]]. The cross-product vector is Xᵀy = [∑yᵢ, ∑xᵢyᵢ]ᵀ = [14, 39]ᵀ. We verify: ∑yᵢ = 2+3+5+4 = 14 and ∑xᵢyᵢ = 1×2 + 2×3 + 3×5 + 4×4 = 2+6+15+16 = 39.

XᵀX = [[4, 10], [10, 30]], Xᵀy = [14, 39]ᵀ

Step 3 — Solve the Normal Equations XᵀXβ̂ = XᵀyWe need (XᵀX)⁻¹. For a 2×2 matrix [[a,b],[c,d]], the inverse is (1/det)[[d,−b],[−c,a]]. Here det(XᵀX) = 4×30 − 10×10 = 120 − 100 = 20. So (XᵀX)⁻¹ = (1/20)[[30, −10], [−10, 4]] = [[1.5, −0.5], [−0.5, 0.2]]. Therefore β̂ = (XᵀX)⁻¹Xᵀy = [[1.5, −0.5], [−0.5, 0.2]] × [14, 39]ᵀ = [1.5×14 + (−0.5)×39, (−0.5)×14 + 0.2×39]ᵀ = [21 − 19.5, −7 + 7.8]ᵀ = [1.5, 0.8]ᵀ.

β̂₀ = 1.5, β̂₁ = 0.8. The fitted line is ŷ = 1.5 + 0.8x.

Step 4 — Compute Fitted Values and Residualsŷ₁ = 1.5 + 0.8(1) = 2.3, ŷ₂ = 1.5 + 0.8(2) = 3.1, ŷ₃ = 1.5 + 0.8(3) = 3.9, ŷ₄ = 1.5 + 0.8(4) = 4.7. Residuals: e₁ = 2 − 2.3 = −0.3, e₂ = 3 − 3.1 = −0.1, e₃ = 5 − 3.9 = 1.1, e₄ = 4 − 4.7 = −0.7.

Check: ∑eᵢ = −0.3 − 0.1 + 1.1 − 0.7 = 0 ✓ (residuals sum to zero when an intercept is included)

Step 5 — Verify Orthogonality and Estimate σ²The normal equations require Xᵀe = 0. We verify: column 1 dot e = 1(−0.3) + 1(−0.1) + 1(1.1) + 1(−0.7) = 0 ✓. Column 2 dot e = 1(−0.3) + 2(−0.1) + 3(1.1) + 4(−0.7) = −0.3 − 0.2 + 3.3 − 2.8 = 0 ✓. The residual sum of squares is SSE = (−0.3)² + (−0.1)² + (1.1)² + (−0.7)² = 0.09 + 0.01 + 1.21 + 0.49 = 1.80. With n − p = 4 − 2 = 2 degrees of freedom, s² = 1.80/2 = 0.90.

s² = 0.90, s = 0.949. The standard errors of β̂ are obtained from √(s² × diag((XᵀX)⁻¹)): SE(β̂₀) = √(0.9 × 1.5) = 1.16, SE(β̂₁) = √(0.9 × 0.2) = 0.424.

SECTION 7

Strengths, Limitations & Computational Considerations

Ordinary least squares via the normal equations is remarkably powerful, but no estimator is universally optimal. Understanding its strengths and weaknesses guides practitioners in selecting appropriate methods for their data.

OLS via Normal Equations: Strengths vs. Limitations

Strengths	Limitations
Closed-form solution β̂ = (XᵀX)⁻¹Xᵀy — no iterative optimization required.	Requires XᵀX to be invertible (full column rank). Fails when p > n or under perfect multicollinearity.
BLUE under Gauss–Markov conditions — optimal among all linear unbiased estimators.	Sensitive to outliers: squared loss magnifies the effect of extreme residuals.
Geometric transparency — projection interpretation provides deep intuition.	Forming XᵀX can worsen conditioning: cond(XᵀX) = cond(X)². QR or SVD approaches are numerically superior.
Under normality of ε, OLS = MLE, and exact t- and F-tests are available.	Assumes homoscedasticity and uncorrelated errors. Violations require GLS, robust, or sandwich methods.
Foundation for extensions: WLS, GLS, ridge, LASSO, and GLMs all build on the OLS framework.	Unbiasedness comes at the cost of potentially high variance when predictors are nearly collinear — the bias-variance tradeoff favors regularized alternatives.

⚙ COMPUTATIONAL NOTE

In practice, statisticians rarely invert XᵀX directly. Instead, the design matrix X is factored as X = QR (QR decomposition) or X = UΣVᵀ (singular value decomposition). From the QR factorization, the normal equations simplify to Rβ̂ = Qᵀy, which is solved via back-substitution. This avoids squaring the condition number and is the standard approach in software such as R's lm() function and Python's numpy.linalg.lstsq().

SECTION 8

Connection to Generalized Least Squares and GLMs

Ordinary least squares is the simplest member of a family of estimation methods. When the assumption Cov(ε) = σ²I is violated — due to heteroscedasticity or correlated errors — the generalized least squares (GLS) estimator replaces the identity with a known positive definite covariance matrix Ω, yielding β̂_GLS = (XᵀΩ⁻¹X)⁻¹XᵀΩ⁻¹y. This can be understood as OLS applied to a transformed model where both sides are premultiplied by Ω^(−1/2). Moving further from linearity, generalized linear models (GLMs) extend the framework by allowing non-Gaussian response distributions and a nonlinear link function g(μ) = Xβ. In GLMs, the parameter estimation is performed via iteratively reweighted least squares (IRLS), where each iteration solves a weighted normal equation with a weight matrix that depends on the current parameter estimates.

OLS → GLS → GLM: A Progressive Generalization

Feature	OLS (Normal Equations)	GLS	GLM (IRLS)
Error distribution	ε ~ (0, σ²I) — any distribution with constant variance	ε ~ (0, σ²Ω) — known covariance structure	y ~ Exponential family (Bernoulli, Poisson, Gamma, etc.)
Normal equations	XᵀXβ̂ = Xᵀy	XᵀΩ⁻¹Xβ̂ = XᵀΩ⁻¹y	XᵀW⁽ᵏ⁾Xβ⁽ᵏ⁺¹⁾ = XᵀW⁽ᵏ⁾z⁽ᵏ⁾ (iterative)
Link function	Identity: E[y] = Xβ	Identity: E[y] = Xβ	g(E[y]) = Xβ (logit, log, inverse, etc.)
Solution type	Closed-form	Closed-form (if Ω known)	Iterative convergence to MLE
Optimality	BLUE (Gauss–Markov)	BLUE under general Ω (Aitken theorem)	Asymptotically efficient (MLE properties)

The thread connecting all three frameworks is the idea of projection. OLS projects onto C(X) using the standard inner product. GLS uses a weighted inner product defined by Ω⁻¹. IRLS projects onto a linearized version of the mean function at each iteration, converging to the maximum likelihood solution. Mastering the OLS normal equations is therefore not merely an exercise in introductory regression—it builds the conceptual foundation for the entire generalized linear modeling framework that dominates modern applied statistics.

SECTION 9

Practice Problems

PROBLEM 1 — CONCEPTUAL

Explain, using the geometric interpretation of least squares, why the residual vector e = y − ŷ must be orthogonal to every column of the design matrix X. What would happen to the fit if this orthogonality condition were violated?

PROBLEM 2 — BASIC CALCULATION

Given X = [[1, 2], [1, 3], [1, 5]] and y = [3, 5, 8]ᵀ, compute XᵀX, Xᵀy, and the OLS estimate β̂ = (XᵀX)⁻¹Xᵀy.

PROBLEM 3 — INTERMEDIATE

For a simple linear regression model y = β₀ + β₁x + ε, show that the normal equations reduce to the familiar scalar formulas β̂₁ = Sxy/Sxx and β̂₀ = ȳ − β̂₁x̄, where Sxy = ∑(xᵢ − x̄)(yᵢ − ȳ) and Sxx = ∑(xᵢ − x̄)².

PROBLEM 4 — APPLIED

A materials scientist models tensile strength (y, in MPa) as a function of heat treatment temperature (x₁, in °C) and alloy composition (x₂, weight %). With n = 20 observations, the computed matrices are XᵀX = [[20, 1000, 50], [1000, 52000, 2600], [50, 2600, 140]] and Xᵀy = [800, 41500, 2150]ᵀ. If a colleague reports that (XᵀX)⁻¹ has been computed, describe how you would obtain β̂, the standard error of each coefficient, and the 95% confidence interval for the slope β̂₁. What additional quantity do you need?

PROBLEM 5 — CRITICAL THINKING

Prove that the hat matrix H = X(XᵀX)⁻¹Xᵀ is idempotent (H² = H) and symmetric (Hᵀ = H). Then explain why these properties imply that all eigenvalues of H are either 0 or 1, and argue that tr(H) = p.

SUMMARY

Lesson Summary

The method of least squares estimates the parameter vector β in the linear model y = Xβ + ε by minimizing the residual sum of squares S(β) = ‖y − Xβ‖². Setting the gradient to zero produces the normal equations XᵀXβ̂ = Xᵀy, which have the unique closed-form solution β̂ = (XᵀX)⁻¹Xᵀy when X has full column rank. Geometrically, the fitted vector ŷ = Xβ̂ is the orthogonal projection of y onto the column space of X, and the residual e = y − ŷ is perpendicular to every column of X — the very condition encoded by the normal equations.

Under the Gauss–Markov assumptions (E[ε] = 0, Cov(ε) = σ²I), the OLS estimator is BLUE — the Best Linear Unbiased Estimator. Its covariance matrix is Cov(β̂) = σ²(XᵀX)⁻¹, and the unbiased estimate of σ² is s² = eᵀe/(n − p). The hat matrix H = X(XᵀX)⁻¹Xᵀ (idempotent, symmetric, with trace p) drives the ANOVA decomposition SST = SSR + SSE. This OLS framework generalizes to GLS (weighted inner products) and GLMs (iteratively reweighted least squares), making it the conceptual bedrock of modern regression analysis.