Linear Regression (OLS)
Squared-error loss, the normal equations, and the Gauss-Markov theorem.
On this page
Beginner: you have a scatter of points — square footage vs. price, hours studied vs. exam score, whatever — and you want to draw the single straight line that best summarizes the trend. "Best" here has a precise meaning: for every point, measure the vertical gap between the point and the line, square that gap, and add up the squares over every point. Ordinary least squares (OLS) is the specific line — really, in more than one input dimension, a hyperplane — that makes that total as small as possible. Nothing about it is iterative or approximate in spirit: there is one line that minimizes total squared vertical distance, and OLS finds it exactly.
Intermediate: why square the gaps instead of, say, just adding up their absolute values? Two reasons, one practical and one statistical. Practically, the squared objective is smooth and differentiable everywhere — including exactly at a residual of zero — which is what lets you set a derivative to zero and solve for the optimum in one shot, in closed form, with no iteration at all (the derivation below does exactly this). Absolute error has a kink at zero and no such closed form. Statistically, squaring a residual penalizes it much more harshly as it grows — a point that misses by 4 units contributes 16 to the total, one that misses by 1 contributes only 1 — so squared error treats a single large miss as far more costly than several small ones, which is exactly the right sensitivity if you believe your noise is roughly Gaussian (made precise in Part 3 of the derivation below) and exactly the wrong sensitivity if your data has occasional wild outliers (the subject of Robust Regression, 2.3.10, later in this module).
Advanced — the geometric view. Stack every row of predictors into a matrix X (rows = data points, columns = features, typically with a column of 1s for the intercept) and every target into a vector y. If you could choose θ freely, Xθ can only ever land somewhere in the column space of X — the set of all vectors reachable as some linear combination of X's columns, exactly the "span" from Linear Algebra for ML, section 1.9. With noisy real data, the actual target vector y essentially never lies exactly in that column space — there's no θ that fits every point perfectly. The best you can do is find the point inside the column space that is closest to y in ordinary Euclidean distance, and that closest point is, by definition, the orthogonal projection of y onto the column space of X. OLS is exactly that projection: it doesn't search or iterate toward it, it computes it directly. This is the same projection machinery built up across section 1.9 (span, basis, rank), section 1.10 (the SVD, which is the general tool for describing what a matrix does to space), and the numerically stable way to actually compute a projection via section 1.14 (QR decomposition) — OLS is the first place in this course where that abstract machinery gets pointed at a concrete, named algorithm.
Worked example. Three points: (1, 2), (2, 3), (3, 7). Averages are x̄ = 2, ȳ = 4. The least-squares slope is the ratio of how x and y co-vary to how much x varies on its own: slope = Σ(x-x̄)(y-ȳ) / Σ(x-x̄)² = [(-1)(-2)+(0)(-1)+(1)(3)] / [1+0+1] = 5/2 = 2.5, and the intercept follows from forcing the line through the averages: intercept = ȳ - slope·x̄ = 4 - 2.5·2 = -1. The fitted line ŷ = -1 + 2.5x predicts 1.5, 4, and 6.5 at the three x-values — residuals of 0.5, -1, and 0.5, which do not all vanish, because three points don't lie on any single line. But among every possible line, this exact one minimizes the sum of squared residuals, and the derivation below is precisely the general-case version of the two formulas just used by hand.
θ̂ is the vector of fitted coefficients (intercept plus one weight per feature). H = X(XᵀX)⁻¹Xᵀ, sometimes called the hat matrix because it puts the hat on y, is exactly the orthogonal projection matrix onto the column space of X described above — it satisfies H² = H and Hᵀ = H, the defining algebraic properties of any orthogonal projection.
Part 1 — deriving the normal equations. The objective is total squared error, written as a single vector norm:
Expand the product term by term:
(the middle two terms are transposes of each other and both scalars, hence equal, hence they combine). Differentiate with respect to the vector θ using two standard matrix-calculus identities — ∇_θ(b^\top θ) = b and, for symmetric A, ∇_θ(θ^\top A θ) = 2Aθ (both covered in Linear Algebra for ML, section 1.17):
Set the gradient to the zero vector — the first-order condition for a minimum:
These are the normal equations. When X has full column rank (no exact linear dependence among its columns — see the multicollinearity pitfall below), XᵀX is invertible and this solves directly to the formula above:
One more check confirms this is really a minimum and not just a stationary point: the Hessian of J is 2XᵀX, which is always positive semi-definite (it's 2AᵀA in disguise, and vᵀAᵀAv = ‖Av‖² ≥ 0 for every v) and strictly positive definite exactly when X has full column rank — the same convexity argument developed in section 1.16 guarantees the normal equations' solution is the unique global minimum of J, not a saddle point.
Part 2 — the Gauss-Markov theorem. Assume the true relationship is linear in the parameters, y = Xθ + ε, under four conditions: (i) linearity as just stated; (ii) exogeneity, E[ε | X] = 0 — the noise carries no systematic relationship to the predictors; (iii) no perfect multicollinearity, i.e. X has full column rank; and (iv) homoscedastic, uncorrelated errors, Cov(ε | X) = σ²I — every error has the same variance and no two errors are correlated with each other. Under exactly these four assumptions, the Gauss-Markov theorem says OLS is BLUE — the Best Linear Unbiased Estimator: among every estimator that is (a) linear in y and (b) unbiased, OLS has the smallest variance.
Proof. Let θ̃ = Cy be any linear estimator, for some matrix C that may depend on X but not on y. Substituting y = Xθ + ε:
using exogeneity. Unbiasedness means this must equal θ for every possible θ, which forces the constraint CX = I. Now write C as the OLS-implied matrix plus a deviation:
Plugging into CX = I shows the deviation D must satisfy DX = 0 — any valid unbiased linear estimator is OLS's matrix plus something that annihilates X. Now compute the variance of θ̃ = Cy using Cov(ε|X) = σ²I:
Expand CCᵀ using C = (XᵀX)⁻¹Xᵀ + D. The two cross terms both vanish, because DX = 0 implies (X^\top X)^{-1}X^\top D^\top = [(X^\top X)^{-1}(DX)^\top] is a product containing DX = 0 (and symmetrically for the other cross term), leaving only the two "diagonal" pieces:
so that
DDᵀ is positive semi-definite for any matrix D (it'sAAᵀ in disguise again), so Var(θ̃) − Var(θ̂_OLS) ⪰ 0 in the matrix sense — every valid linear unbiased estimator has variance at least as large as OLS's, in every direction, with equality exactly when D = 0, i.e. when θ̃ = θ̂_OLS. That is the Gauss-Markov theorem: OLS is the unique minimum-variance member of the entire class of linear unbiased estimators.
Part 3 — the Gaussian-noise MLE view. Module 1's Maximum Likelihood & MAP lesson (2.1.6) proved the loss ⟺ likelihood correspondence in the abstract; here it's redone specifically for regression. Strengthen assumption (iv) above to a full distribution rather than just a mean and variance — assume the noise is exactly Gaussian, ε_i \sim \mathcal N(0, \sigma^2), independent across data points. Then y_i | x_i, \theta is Gaussian with mean x_i^\top\theta:
Assuming independence across points, sum the log of this density over the whole dataset:
The first term has no θ in it at all — a constant as far as the maximization goes — and 1/(2σ²) is a positive constant multiplier, so:
— exactly the OLS objective from Part 1. Maximizing the likelihood of the data under a linear-Gaussian-noise model is algebraically identical to minimizing squared error; they are the same optimization problem, not two separately-motivated ideas that happen to agree. In particular, this means the normal-equations solution θ̂ from Part 1 is the maximum-likelihood estimator whenever the noise really is Gaussian — squared error was never an arbitrary choice, it's the log-likelihood of the single most common noise assumption in all of statistics.
Where this is used: every method in this Regression module is a deliberate edit to one of these three results. Ridge (2.3.2) and Lasso (2.3.3) add a penalty term to the Part 1 objective — the module overview's L(y,Xθ) + λR(θ) template with R now nonzero — which breaks unbiasedness on purpose, trading some of the Part 2 guarantee away for lower variance. Weighted/generalized least squares relax assumption (iv)'s homoscedasticity requirement and reweight the normal equations accordingly. Generalized Linear Models (2.3.9) redo Part 3 with a non-Gaussian exponential- family distribution in place of the Gaussian, and logistic regression is the Bernoulli special case. Understanding exactly which assumption each later method changes — and which three it leaves alone — is the fastest way to keep fifteen regression methods from feeling like fifteen unrelated formulas.
On load, the blue line animates from a deliberately bad guess to the true OLS fit computed live from the normal equations, while the sum of squared residuals (dashed gray segments, squared and summed) shrinks to its minimum. Drag any teal point afterward and the line refits instantly.
Nested ellipses are the exact level sets of the sum-of-squared-residuals surface J(θ₀,θ₁), traced via a closed-form eigendecomposition of its quadratic form -- their being closed loops (not an open trough) is exactly the positive-definite-Hessian fact from Part 1 of the derivation. The intro animation runs plain gradient descent from a bad guess and zig-zags slowly down the shallow valley (intercept and slope are correlated because x isn't centered), never quite catching the red closed-form minimum in the time it's given -- unlike the normal equations, which land there in a single step. Drag the blue probe afterward to read off any point's SSE.
All three tabs fit the same ten points and land on the same intercept, slope, and sum of squared residuals. The from-scratch tabs solve the normal equations directly (a general linear solve in Python, an explicit 2x2 Cramer's-rule solve in C++); the library tab confirms scikit-learn is computing the identical mathematical object, just with a production-grade solver underneath.
#include <cstdio>
#include <cmath>
#include <vector>
int main() {
std::vector<double> xs = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
std::vector<double> ys = {2.1, 3.9, 6.2, 7.8, 10.1, 11.9, 14.2, 15.8, 18.1, 19.9};
int n = static_cast<int>(xs.size());
// Normal equations for simple linear regression collapse to a 2x2 system:
// [ n Sx ] [b0] [ Sy ]
// [ Sx Sxx ] [b1] = [ Sxy ]
double Sx = 0.0, Sy = 0.0, Sxx = 0.0, Sxy = 0.0;
for (int i = 0; i < n; ++i) {
Sx += xs[i];
Sy += ys[i];
Sxx += xs[i] * xs[i];
Sxy += xs[i] * ys[i];
}
// Explicit 2x2 solve via Cramer's rule -- no linear-algebra library at all.
double det = n * Sxx - Sx * Sx;
double b0 = (Sxx * Sy - Sx * Sxy) / det; // intercept
double b1 = (n * Sxy - Sx * Sy) / det; // slope
std::printf("OLS via normal equations: intercept=%.4f slope=%.4f\n", b0, b1);
double sse = 0.0;
for (int i = 0; i < n; ++i) {
double resid = ys[i] - (b0 + b1 * xs[i]);
sse += resid * resid;
}
std::printf("Sum of squared residuals: %.4f\n", sse);
return 0;
}- Hedonic pricing models — regressing home sale price on square footage, bedroom count, and location indicators to estimate how much each attribute contributes to price. The fitted coefficient on "one extra bedroom" is a genuinely useful, directly interpretable number precisely because OLS is linear in the parameters, not just a black-box prediction.
- Econometric wage equations — the classic Mincer equation regresses log-wage on years of education and years of experience. This is a textbook Gauss-Markov setting: economists rely explicitly on the BLUE guarantee to argue the estimated "return to a year of education" is the most precise unbiased estimate obtainable from a linear model, given the stated assumptions.
- A/B test effect estimation — regressing an outcome metric on a treatment indicator (0/1) directly estimates the average treatment effect as that indicator's coefficient; adding pre-treatment covariates to the same regression (CUPED-style variance reduction) keeps the estimate unbiased while shrinking its variance, exactly the trade-off Gauss-Markov formalizes.
- Sensor and instrument calibration — plotting known reference values against raw instrument readings and fitting a line gives a calibration curve (slope and offset) used to correct every future raw reading from that instrument back to true units.
- Fitting physical laws from measurements — Hooke's law (spring extension vs. applied force) or Ohm's law (current vs. voltage) are both linear relationships whose slope is the physical constant of interest (spring stiffness, resistance); OLS on repeated noisy measurements is the standard way to estimate that constant along with its uncertainty.
- Marketing mix modeling — regressing sales on spend across several advertising channels to estimate each channel's incremental contribution. This is also the textbook setting where the multicollinearity pitfall below bites hardest: ad spend across channels is frequently launched in correlated bursts, which is exactly the situation the next lesson's ridge penalty is designed to fix.
- Multicollinearity. If two or more columns of
Xare nearly (not even exactly) linearly dependent,XᵀXbecomes nearly singular — its determinant approaches zero, its inverse blows up, and the fitted coefficients become wildly unstable, with huge variance and swings in sign from tiny changes in the data, even though Gauss-Markov still calls the estimator unbiased. Exact linear dependence makes(XᵀX)⁻¹outright undefined. This single failure mode is the entire motivation for the very next lesson, Ridge Regression (2.3.2), which trades a small, controlled amount of bias for a large reduction in this variance. - Extrapolating beyond the training range. The fitted line is only ever justified by data inside the range (more precisely, the convex hull) of the observed predictors. Nothing about the model prevents you from evaluating it far outside that range, and nothing guarantees the true relationship stays linear out there — a calibration line fit between 0°C and 40°C says nothing reliable about behavior at 200°C.
- Treating a high R² as proof of a good model. R² measures only how well the fitted line matches this particular training sample — it says nothing about whether the Gauss-Markov assumptions actually hold (autocorrelated residuals and omitted relevant variables can both inflate it misleadingly), and it mechanically never decreases when you add more features to the same training data, regardless of whether those features have any real predictive value. A high R² and a well-specified, generalizing model are related but very much not the same claim.
Going deeper
Production numerical code essentially never computes θ̂ the way the formula literally reads — by forming XᵀX, inverting it, and multiplying. Forming XᵀX squares the condition number of the problem (cond(XᵀX) = cond(X)²), which is precisely why multicollinearity above is such a numerically dangerous pitfall — a moderately ill-conditioned X becomes a severely ill-conditioned XᵀX. Instead, real solvers factor X directly: QR decomposition (Linear Algebra for ML, section 1.14) turns the least- squares problem into a simple triangular back-substitution without ever forming XᵀX, and an SVD-based solve (section 1.10) goes further, producing the minimum-norm solution even when X doesn't have full column rank — a case where (XᵀX)⁻¹ doesn't exist at all but a sensible "best" answer still does. This is exactly what a call like NumPy's lstsq or scikit-learn's LinearRegression does under the hood; the normal-equations formula above is the right way to understand OLS, and the wrong way to actually compute it at scale.
Gauss-Markov says OLS is BLUE — the Best Linear Unbiased Estimator. Does that guarantee no other estimator, linear or not, could ever have lower variance than OLS?
No. BLUE is a best-in-class result restricted specifically to estimators that are both linear in y and unbiased — it says nothing about estimators outside that class. A biased estimator can have strictly lower variance, and sometimes lower total mean squared error (variance plus squared bias) as a result — this is exactly the trade-off ridge regression (2.3.2) makes deliberately, giving up unbiasedness to shrink variance enough that total error drops. Nonlinear estimators can also beat OLS if they exploit information Gauss-Markov's assumptions don't use — e.g. the true maximum-likelihood estimator under a known non-Gaussian noise distribution generally outperforms OLS in that specific setting, even though OLS remains the best available option among estimators restricted to being linear and unbiased.
OLS solves one clean problem exactly: minimize total squared vertical distance, which is the same as projecting y orthogonally onto the column space of X, which is the same as differentiating ‖y−Xθ‖² to zero to get the normal equations θ̂ = (XᵀX)⁻¹Xᵀy. Under four standard assumptions, Gauss-Markov guarantees that solution has the lowest possible variance of any linear unbiased estimator; under the stronger assumption of Gaussian noise, that same solution is exactly the maximum- likelihood estimate. Every later lesson in this Regression module is a deliberate, motivated departure from one part of this picture — a penalty added to the objective, an assumption relaxed, a distribution generalized — so the clearest way to understand any of them is to ask exactly which piece of this lesson they're changing, and why.