The Mathematics of Linear Regression

Almost every model you will ever train is fitted by taking small steps downhill: guess some parameters, measure the error, nudge, repeat, and hope it settles somewhere good. Linear regression is the rare exception. For this one model you can write down the exact optimal parameters in closed form, in a single line of algebra, with no iteration and no learning rate and no guessing.

That is worth understanding properly, and not because you will implement it by hand — a library call does it in one line. It is worth understanding because the derivation is where the pieces you use everywhere else fit together visibly: a loss function, a gradient, a critical point, a matrix inverse, and a geometric picture that turns out to explain what "fitting" means in general. It is also the cleanest place to see why the mathematically obvious way to compute something can be the numerically wrong way.

Setting the problem up

You have $n$ observations and $d$ features. Stack them into a design matrix $\mathbf{X}$ with one row per observation and one column per feature:

The model is

Each prediction is a weighted sum of that row's features. Multiplying the whole matrix at once produces all $n$ predictions in one operation.

The intercept — the constant offset that lets the fitted line avoid passing through the origin — is handled by a trick: add a column of ones to $\mathbf{X}$. Then the corresponding weight is multiplied by 1 for every observation, so it acts as a constant, and no special case is needed anywhere in the algebra.

The objective is mean squared error:

In words: for each observation, take how far the prediction missed, square it so that overshooting and undershooting both count as error, and average. Squaring also makes the function smooth and differentiable everywhere, which is what lets the calculus below work. The $1/n$ scales the result to a per-example average; it makes no difference to which $\mathbf{w}$ is best, since scaling a function by a positive constant does not move its minimum.

Deriving the closed-form solution

The minimum is where the gradient is zero. So expand the loss, differentiate it, set that to zero, and solve.

Step 1: expand the squared norm

A squared norm is a dot product of a vector with itself, so $\|\mathbf{v}\|^2 = \mathbf{v}^T\mathbf{v}$. Applying that and multiplying out:

The two cross terms $\mathbf{y}^T\mathbf{X}\mathbf{w}$ and $\mathbf{w}^T\mathbf{X}^T\mathbf{y}$ combined into one because each is a single number, and a number equals its own transpose — so they are equal and add to twice one of them.

Three terms, each with a role. The first does not contain $\mathbf{w}$ at all, so it is a constant as far as the optimisation is concerned. The second is linear in $\mathbf{w}$. The third is quadratic in $\mathbf{w}$. The whole loss is a bowl.

Step 2: two identities you need

Differentiating with respect to a vector means collecting the partial derivative with respect to each component into a vector of the same shape. Two patterns cover everything here.

The symmetry condition is satisfied: $\mathbf{X}^T\mathbf{X}$ is symmetric for any $\mathbf{X}$, because transposing it gives $\mathbf{X}^T(\mathbf{X}^T)^T = \mathbf{X}^T\mathbf{X}$ back.

Step 3: take the gradient

Dropping the constant $1/n$ and differentiating term by term:

This vector has one entry per feature, and each entry says how much the total squared error would grow if you nudged that feature's weight upward.

Step 4: set it to zero and solve

That middle form is called the normal equations, and it is the form you should actually remember, for reasons that become clear later. Assuming $\mathbf{X}^T\mathbf{X}$ is invertible:

No iteration, no learning rate, no convergence criterion. Feed in the data and the optimal weights come out.

Confirming it is a minimum

A zero gradient marks a minimum, a maximum, or a saddle. To tell which, look at the second derivative — here the Hessian matrix:

For any vector $\mathbf{z}$, $\mathbf{z}^T\mathbf{X}^T\mathbf{X}\mathbf{z} = (\mathbf{X}\mathbf{z})^T(\mathbf{X}\mathbf{z}) = \|\mathbf{X}\mathbf{z}\|^2 \ge 0$. A squared length can never be negative, so the Hessian is positive semi-definite: the surface curves upward, or is flat, in every direction. Never downward.

The loss surface of linear regression is a bowl. There is one minimum, it is global, and the formula finds it exactly.

If $\mathbf{X}$ has full column rank — no feature is an exact linear combination of the others — then $\|\mathbf{X}\mathbf{z}\| > 0$ for all non-zero $\mathbf{z}$, the Hessian is strictly positive definite, and the solution is unique. If two features are perfectly collinear, the bowl has a flat trough along one direction, infinitely many weight vectors achieve the same minimum loss, and $\mathbf{X}^T\mathbf{X}$ is singular so the inverse does not exist.

All the way through, with numbers

Three observations, one feature $x$, plus an intercept.

With the ones column first:

Compute $\mathbf{X}^T\mathbf{X}$. Its entries are $n$, $\sum x$, and $\sum x^2$:

Compute $\mathbf{X}^T\mathbf{y}$. Its entries are $\sum y = 10$ and $\sum xy = 2 + 6 + 15 = 23$:

Invert. The determinant is $3(14) - 6(6) = 42 - 36 = 6$, and for a $2\times2$ you swap the diagonal and negate the off-diagonal:

Multiply.

So the fitted line is $\hat y = 0.3333 + 1.5x$. Check it against the standard slope formula: $\frac{n\sum xy - \sum x\sum y}{n\sum x^2 - (\sum x)^2} = \frac{69 - 60}{42 - 36} = 1.5$. Agreed.

Look at the residuals. Predictions are 1.8333, 3.3333, 4.8333, giving residuals $+0.1667$, $-0.3333$, $+0.1667$. Two facts about them are not coincidences:

• They sum to zero: $0.1667 - 0.3333 + 0.1667 = 0$.
• Weighted by $x$ they also sum to zero: $1(0.1667) + 2(-0.3333) + 3(0.1667) = 0.1667 - 0.6667 + 0.5 = 0$.

These are exactly the normal equations written out: $\mathbf{X}^T(\mathbf{y} - \mathbf{X}\mathbf{w}) = \mathbf{0}$ says every column of $\mathbf{X}$ has zero dot product with the residual vector. The first column is the ones, giving the first fact; the second column is $x$, giving the second.

The geometry: least squares is a projection

This is the picture that makes everything above feel inevitable rather than algebraic.

$\mathbf{y}$ is a point in $n$-dimensional space — one dimension per observation. Your predictions $\mathbf{X}\mathbf{w}$ can only ever be linear combinations of the columns of $\mathbf{X}$, so as $\mathbf{w}$ varies over all possible values, $\mathbf{X}\mathbf{w}$ sweeps out a $d$-dimensional flat subspace: the column space of $\mathbf{X}$.

In our example, $\mathbf{y} = [2,3,5]^T$ lives in 3-dimensional space, and the achievable predictions form a 2-dimensional plane spanned by $[1,1,1]^T$ and $[1,2,3]^T$. Since $\mathbf{y}$ does not lie on that plane — no straight line passes through all three points — there is no exact solution.

So you settle for the closest achievable point. And the closest point of a plane to an outside point is its orthogonal projection: drop a perpendicular. That perpendicular is the residual vector, which is therefore at right angles to the plane — at right angles to every column of $\mathbf{X}$ — which is precisely $\mathbf{X}^T(\mathbf{y} - \mathbf{X}\mathbf{w}) = \mathbf{0}$.

The normal equations are not an algebraic accident. They are the statement that the residual is perpendicular to everything your features can express.

That orthogonality has a plain-English meaning worth holding on to: there is no signal left in the residual that your features could have captured. If the residual still correlated with any feature, you could adjust that feature's weight and reduce the error. At the optimum, every last drop of linearly-extractable information has been squeezed out. What remains is either noise or structure your features simply cannot represent.

The hat matrix

Substituting the solution back gives the fitted values directly from $\mathbf{y}$:

$\mathbf{H}$ is called the hat matrix because it puts the hat on $\mathbf{y}$. It is the projection operator onto the column space, and it has two telling properties. It is symmetric, and it is idempotent: $\mathbf{H}^2 = \mathbf{H}$. Projecting something that is already on the plane leaves it alone — exactly what you would expect from a shadow.

Its diagonal entries $h_{ii}$ measure how much observation $i$ influences its own fitted value. High-leverage points, with $h_{ii}$ close to 1, can drag the whole fit towards themselves. The trace of $\mathbf{H}$ equals $d$, the number of parameters, which is the formal sense in which "degrees of freedom used by the model" equals "number of features".

What has to be true for the estimate to be trustworthy

The formula always returns numbers. Whether those numbers mean anything depends on assumptions.

When the first five hold, the Gauss-Markov theorem says the least-squares estimator is BLUE: the Best Linear Unbiased Estimator. Among all estimators that are linear in $\mathbf{y}$ and unbiased, it has the smallest variance. Note the two qualifiers, because they are where the interesting alternatives live: drop "unbiased" and you can do better, which is the whole justification for shrinkage methods that accept a little bias in exchange for a large reduction in variance.

Unbiasedness, shown

Substitute the true model $\mathbf{y} = \mathbf{X}\mathbf{w}^* + \boldsymbol{\epsilon}$ into the formula:

The first term collapsed because $(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{X} = \mathbf{I}$. Taking expectations and using $E[\boldsymbol{\epsilon}] = \mathbf{0}$, the second term vanishes and $E[\hat{\mathbf{w}}] = \mathbf{w}^*$. On average, across many datasets, you land on the truth.

Variance, and why collinearity hurts

The diagonal entries are the variances of the individual coefficient estimates. For our worked example, $(\mathbf{X}^T\mathbf{X})^{-1} = \begin{bmatrix} 2.333 & -1 \\ -1 & 0.5\end{bmatrix}$, so the slope estimate has variance $0.5\sigma^2$ and the intercept $2.333\sigma^2$.

Now consider what happens as two features become nearly collinear. $\mathbf{X}^T\mathbf{X}$ approaches singularity, its determinant approaches zero, and its inverse blows up. The variances explode. Concretely: the model cannot tell which of the two near-identical features deserves the credit, so it distributes weight between them almost arbitrarily — you might get $+800$ and $-798$ on features that differ by a rounding error. Refit on slightly different data and those numbers change completely, even though the predictions barely move.

This is the practical reason to care about the determinant: it is not an abstract quantity. It sits in the denominator of the inverse, so as it shrinks towards zero the coefficient variances grow in proportion to $1/\det$.

Why nobody computes that inverse

The formula says $(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$. Writing that literally in code is a mistake, for two reasons.

Cost. Inverting a $d \times d$ matrix takes $O(d^3)$ operations, and you then throw the inverse away after one multiplication. Solving the linear system $\mathbf{X}^T\mathbf{X}\mathbf{w} = \mathbf{X}^T\mathbf{y}$ directly costs the same order but with a smaller constant and much better numerical behaviour.

Conditioning. This is the serious one. The condition number of a matrix measures how much it amplifies input errors; a condition number of $10^k$ means you can lose about $k$ digits of accuracy. And the key fact is:

Forming $\mathbf{X}^T\mathbf{X}$ squares the conditioning problem. If $\mathbf{X}$ has a moderate condition number of $10^4$ — entirely ordinary with features on different scales — then $\mathbf{X}^T\mathbf{X}$ has $10^8$. Double precision carries about 16 significant digits, so you have thrown away half of them before doing any useful work.

The normal equations are correct mathematics and poor numerics. Solve the least-squares problem directly from $\mathbf{X}$, without ever forming $\mathbf{X}^T\mathbf{X}$.

Seeing the difference in code

1import numpy as np23rng = np.random.default_rng(0)45# A deliberately ill-conditioned design: the third column is almost6# a copy of the second, so the columns are nearly linearly dependent.7n = 2008x1 = rng.normal(size=n)9x2 = x1 + 1e-8 * rng.normal(size=n)10X = np.column_stack([np.ones(n), x1, x2])11w_true = np.array([1.0, 2.0, -1.0])12y = X @ w_true + 0.01 * rng.normal(size=n)1314print("cond(X)      =", f"{np.linalg.cond(X):.3e}")15print("cond(X.T@X)  =", f"{np.linalg.cond(X.T @ X):.3e}")   # roughly the square1617# 1. Explicit inverse -- mathematically right, numerically worst.18w_inv = np.linalg.inv(X.T @ X) @ X.T @ y1920# 2. Solve the normal equations without forming the inverse -- better.21w_solve = np.linalg.solve(X.T @ X, X.T @ y)2223# 3. SVD-based least squares -- never forms X.T @ X at all.24w_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)2526# 4. QR decomposition -- the classical stable approach.27Q, R = np.linalg.qr(X)28w_qr = np.linalg.solve(R, Q.T @ y)2930for name, w in [("inverse", w_inv), ("solve", w_solve),31                ("lstsq ", w_lstsq), ("qr    ", w_qr)]:32    resid = np.linalg.norm(y - X @ w)33    print(f"{name}  w = {np.round(w, 4)}   ||residual|| = {resid:.6f}")

Two things to notice when you run this. First, cond(X.T @ X) is roughly the square of cond(X), exactly as the theory predicts. Second, the second and third coefficients come out in the tens of thousands with opposite signs, and at different values for different methods, yet solve, lstsq and qr reach almost the same residual norm. That is the signature of an ill-conditioned problem: many very different weight vectors fit the data almost equally well, and which one you get depends on your numerical method rather than on the data. The explicit inverse is the odd one out: its residual is usually far larger (on one run, 5.2 against 0.14), which is what losing half your digits before you start looks like.

The individual coefficients are meaningless in that situation. If someone hands you a regression with a coefficient of $+800$ on one feature and $-798$ on a near-duplicate, and interprets those numbers as effects, they are interpreting floating-point noise.

What to do with this when you fit a real model

Call np.linalg.lstsq or a library's LinearRegression, which uses a stable decomposition internally. Do not write the inverse formula, however satisfying it looks.

Before fitting, check the condition number of your design matrix. Below about $10^3$, relax. Between $10^3$ and $10^6$, be careful about interpreting individual coefficients. Above $10^6$, the coefficients are not to be trusted at all — drop or combine the redundant features, or add a small penalty on the weights, which lifts the smallest eigenvalues away from zero and stabilises everything.

Standardise your features first when you intend to compare coefficient magnitudes. A coefficient's size depends entirely on its feature's units; a weight of 0.0001 on a variable measured in millimetres and a weight of 100 on the same variable in kilometres describe identical relationships.

And always inspect the residuals rather than only the error number. The orthogonality result says a correct fit leaves residuals with no linear relationship to any feature. So plot residuals against each feature and against the fitted values. A pattern — a curve, a fan shape, a trend — means an assumption has failed: nonlinearity, or non-constant variance, or a missing feature. A well-fitted model leaves a formless cloud, and that formlessness is precisely the geometric statement that the projection was done correctly.