How Linear Regression Works: Least Squares, Normal Equations and Gradient Descent

Linear regression is usually the first model people learn and the one they keep using. It is also easy to use without understanding, so this post walks through what the fitting procedure actually does: the cost function, the normal equations and gradient descent. It was written with Bohdan Boiprav.

A note on the old title: the original version of this post referred to the “statistical mechanics” of linear regression. That is a different field, a branch of physics. What we mean here is the statistics and linear algebra of least squares.

The model

Linear regression models a numeric target y as a weighted sum of inputs plus an intercept:

y_hat = b0 + b1*x1 + b2*x2 + ... + bp*xp

With one input, this is a line. With two it is a plane, and with more it is a hyperplane. Fitting the model means choosing the coefficients b0 ... bp. In investment data work, a familiar example is regressing a stock’s returns on market returns: the slope is the stock’s beta and the intercept is its alpha over the sample.

Residuals and the cost function

A residual is the vertical gap between an observed value and the fitted line: e_i = y_i - y_hat_i. Ordinary least squares (OLS) picks the coefficients that make the sum of squared residuals as small as possible:

J(b) = sum_i (y_i - y_hat_i)^2

In machine-learning terms, J is the cost (or loss) function. Squaring does two things. It stops positive and negative errors cancelling out, and it makes the function smooth and convex, so there is a single minimum and no local traps. The price is that large errors count disproportionately, which makes OLS sensitive to outliers.

The normal equations

Write the data as a matrix X with one row per observation and a leading column of ones for the intercept. Setting the partial derivative of J with respect to every coefficient to zero gives a set of linear equations, known as the normal equations:

(X'X) b = X'y        so        b = (X'X)^-1 X'y

The earlier draft called these the “equations of normalisation”, which is not a standard term. The closed-form solution exists as long as X'X is invertible, which fails when one input is a perfect linear combination of others (perfect multicollinearity). In practice you rarely invert X'X directly. Solvers use a QR or SVD factorisation, which is more numerically stable when inputs are nearly collinear. That is what numpy.linalg.lstsq does, and scikit-learn’s LinearRegression relies on the same least-squares machinery.

Gradient descent

With a very large number of rows or columns, solving the normal equations can be expensive, and many other models have no closed form at all. Gradient descent is the general alternative. It starts from some coefficients and repeatedly steps in the direction in which the cost falls fastest, which is the negative gradient:

b := b - eta * grad J(b)       where   grad J(b) = (2/n) * X'(Xb - y)  for the mean squared error

Here eta is the learning rate. If it is too large the updates overshoot and diverge. If it is too small, convergence is slow. Gradient descent is also sensitive to the scale of the inputs, so standardise them first. The earlier draft described the update as moving “in the direction of maximum descent on the vector”, which is a garbled version of this idea.

All three routes give the same answer

The snippet below fits the same synthetic data (true coefficients 1.5, 2.0 and -0.5, plus noise) three ways.

import numpy as np

rng = np.random.default_rng(0)
n = 200
x1, x2 = rng.normal(size=n), rng.normal(size=n)
y = 1.5 + 2.0 * x1 - 0.5 * x2 + rng.normal(scale=0.3, size=n)
X = np.column_stack([np.ones(n), x1, x2])        # first column = intercept

# 1. Normal equations: solve (X'X) beta = X'y
beta_ne = np.linalg.solve(X.T @ X, X.T @ y)

# 2. Library least squares (SVD-based, numerically safer)
beta_ls, *_ = np.linalg.lstsq(X, y, rcond=None)

# 3. Gradient descent on the mean squared error
beta = np.zeros(3)
lr = 0.1
for _ in range(500):
    grad = (2 / n) * X.T @ (X @ beta - y)
    beta -= lr * grad

print(beta_ne.round(3), beta_ls.round(3), beta.round(3))

All three print [ 1.503 2.02 -0.482], close to the true values. The small differences from 1.5, 2.0 and -0.5 come from the noise in the sample, not from the method.

What the fit does not guarantee

Getting coefficients is the easy part. For them to carry meaning, and for standard errors and confidence intervals to be trustworthy, the usual assumptions need to hold reasonably well:

• The relationship is linear in the coefficients (you can still include transformed inputs such as x^2 or log(x)).
• Residuals have roughly constant variance. Financial returns often violate this, and robust standard errors are a common fix.
• Residuals are not correlated with each other, which is a frequent issue in time series.
• Inputs are not almost perfectly correlated. Otherwise the coefficients become unstable even though predictions may stay fine.

A plot of residuals against fitted values is the quickest check. A random cloud around zero is what you want. A curve, a funnel or a run of same-signed residuals each point to a specific problem. For inference rather than prediction, a library such as statsmodels reports standard errors, t-statistics and diagnostic tests directly.

Summary

Least squares minimises the sum of squared residuals. The normal equations give the exact answer in one step, and gradient descent reaches it iteratively, which is useful at scale or when no closed form exists. The fitting is just linear algebra. The harder work is checking that the assumptions behind the numbers are reasonable.