How does Ordinary Least Squares derive the coefficient vector, and what is the closed-form solution?

How to think about it

The direct answer

Ordinary Least Squares, or OLS, chooses the coefficient vector that minimizes the sum of squared prediction errors. For a feature matrix X with full column rank, meaning no feature column is an exact linear combination of the others, the unique solution is:

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

The derivation comes from setting the gradient of the squared-error loss to zero. The important production nuance is that we usually do not calculate the inverse explicitly; we solve the least-squares problem with QR or SVD instead.

Why this is the right objective

Suppose there are n observations and p features:

• X is the n by p design matrix.
• y is the vector of n observed target values.
• β is the vector of p coefficients.
• Xβ is the vector of predictions.

If the model includes an intercept, the first column of X is usually all ones. For example, with one feature, the rows might look like [1, 4.2], [1, 5.1], and [1, 6.7]. The first coefficient is the intercept; the second is the feature slope.

The residual vector is the prediction error for every observation:

e = y - Xβ

OLS minimizes the squared length of that residual vector:

L(β) = ||y - Xβ||²

The square is not decorative. It makes positive and negative errors unable to cancel, gives larger errors more influence, and produces a smooth quadratic objective that can be differentiated. That last property is what gives us the closed form.

The derivation

Start by writing the squared norm as a matrix product:

L(β) = (y - Xβ)ᵀ(y - Xβ)

Expand it:

L(β) = yᵀy - 2βᵀXᵀy + βᵀXᵀXβ

The middle terms combine because each is a scalar and therefore equal to its transpose.

Now take the gradient, which is the vector of partial derivatives with respect to every element of β:

∇L(β) = -2Xᵀy + 2XᵀXβ

At the minimum, the gradient is zero:

-2Xᵀy + 2XᵀXβ = 0

Move one term to the other side and divide by 2:

XᵀXβ = Xᵀy

These are the normal equations, the linear equations whose solution minimizes the squared-error objective.

If XᵀX is invertible, multiply both sides by its inverse:

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

That is the familiar closed-form solution.

Why does a zero gradient give the global minimum rather than a random stationary point? Because the objective is a convex quadratic. Its curvature is governed by XᵀX, which is positive semidefinite. With full column rank, it is positive definite, so the quadratic has one unique bottom.

A concrete calculation

Consider a model predicting a delivery time from one numeric feature:

Feature value	Observed time
1	2
2	3
3	5

Include an intercept, so:

X = [[1, 1], [1, 2], [1, 3]]

and:

y = [2, 3, 5]

The matrix products are:

XᵀX = [[3, 6], [6, 14]]

Xᵀy = [10, 23]

So the normal equations are:

3β₀ + 6β₁ = 10

6β₀ + 14β₁ = 23

Solving them gives:

β₀ = 1/3

β₁ = 3/2

The fitted line is therefore:

ŷ = 1/3 + 1.5x

The predictions are approximately [1.833, 3.333, 4.833]. The residuals are approximately [0.167, -0.333, 0.167], and their squared sum is about 0.167.

The calculation can be reproduced with NumPy:

import numpy as np

X = np.array([[1., 1.],
              [1., 2.],
              [1., 3.]])
y = np.array([2., 3., 5.])

beta = np.linalg.lstsq(X, y, rcond=None)[0]
print(np.round(beta, 6))
# [0.333333 1.5     ]

The geometric meaning

The normal equations become clearer if we write the residual as e = y - Xβ:

Xᵀe = 0

Each column of X has a dot product of zero with the residual. In plain language, once OLS has chosen the best coefficients, the remaining error is orthogonal to every direction the model knows how to use.

The columns of X span a subspace of possible prediction vectors. OLS projects y onto that subspace and calls the projected vector ŷ. This is why OLS is often described as a projection.

The hat matrix, which maps observed targets to fitted targets, is:

H = X(XᵀX)⁻¹Xᵀ

It is called the hat matrix because:

ŷ = Hy

Do not confuse this with the coefficient formula. β contains the model coefficients; H maps the target vector to predictions. The hat matrix is useful for leverage and influence diagnostics, but it is usually not formed in production because it is an n by n matrix.

The senior-level caveat: the inverse is not the implementation

The formula is exact in algebra, but explicitly computing (XᵀX)⁻¹ is usually a poor numerical method.

If two feature columns are nearly duplicates, XᵀX becomes ill-conditioned, meaning small rounding errors in the input can cause large changes in the computed coefficients. Forming XᵀX also squares the condition number: in the usual two-norm, cond(XᵀX) = cond(X)². That can throw away meaningful precision.

A least-squares solver avoids this explicit inverse:

beta = np.linalg.lstsq(X, y, rcond=None)[0]

Such routines use a factorization-based method, commonly QR or SVD depending on the implementation. These methods are more stable because they work with the structure of X directly. QR is often a good general-purpose choice; SVD is especially useful when diagnosing rank deficiency.

The computational trade-off also matters:

• Forming XᵀX costs roughly O(np²), then solving the resulting system costs roughly O(p³).
• One gradient-descent pass costs roughly O(np), but many passes may be needed.
• Storing XᵀX requires O(p²) memory, which becomes painful when the feature count is large.

For a dataset with millions of rows and a few dozen dense features, a direct least-squares solver may be perfectly sensible. For hundreds of thousands of sparse features, iterative methods are usually more appropriate. “Closed form” does not automatically mean “fastest.”

What happens when the matrix is singular?

If one feature is exactly duplicated by another, or if there are more features than independent observations, X does not have full column rank. Then XᵀX is singular and the inverse does not exist.

The predictions may still be unique, but the individual coefficients are not. For example, if two columns are identical, increasing one coefficient by 2 and decreasing the other by 2 leaves every prediction unchanged. The model cannot decide how to split the contribution between them.

The general expression is:

β = X⁺y

where X⁺ is the Moore–Penrose pseudoinverse, a generalized inverse that returns the minimum-length least-squares solution.

In practice, inspect the reported rank and singular values, remove or combine redundant features, or use regularization. Ridge regression changes the objective by adding a penalty on coefficient size and gives:

βridge = (XᵀX + λI)⁻¹Xᵀy

for λ > 0. The penalty usually makes the system better conditioned, but it also introduces bias. That is a deliberate bias–variance trade-off, not a free numerical fix.

A common first symptom of collinearity is not a terrible validation score. It is wildly unstable coefficients: a tiny change in the training sample makes one coefficient jump from 12 to -40 while another moves in the opposite direction. Predictions can remain fairly stable because the combined contribution is what the model identifies reliably.

What they’ll ask next

Why not use gradient descent for every linear regression?

Gradient descent works because this objective is convex, but it needs a learning rate, stopping rule, and possibly feature scaling. For modest p, a factorization-based least-squares solve is deterministic and usually simpler. For very large or sparse problems, iterative methods can use less memory and avoid materializing XᵀX.

Do we need Gaussian errors for the derivation?

No. The derivation is purely an optimization argument. Gaussian errors make OLS equivalent to maximum likelihood and support familiar exact inference formulas. Unbiasedness instead depends on assumptions such as the errors having conditional mean zero given X. Heteroskedastic or correlated errors can leave the coefficient calculation unchanged while making ordinary standard errors unreliable.

What is the first thing you would check if OLS gives a singular-matrix error?

I would check duplicate or linearly dependent features, the intercept encoding, missing-value handling, and the matrix rank or singular values. I would not “fix” it by adding a tiny number to the diagonal without understanding the cause. That is regularization, and it changes the model.

Say this in the interview

“OLS minimizes ||y - Xβ||²; differentiating gives the normal equations XᵀXβ = Xᵀy, so under full column rank the closed form is β = (XᵀX)⁻¹Xᵀy, although in production I solve the system with QR or SVD rather than explicitly forming the inverse.”