Linear Regression, in depth: least squares from first principles, a worked example, QR versus normal equations, ridge, lasso and diagnostics

The model and the loss

Given n examples with p features, stack the features into a matrix X of shape n by p, with a column of ones for the intercept, and the targets into a vector y. The model predicts y_hat = X w for a weight vector w. The residual for example i is r_i = y_i - y_hat_i.

Ordinary least squares chooses w to minimise the sum of squared residuals, L(w) = ||y - X w||^2. Squaring has three justifications. It makes the loss smooth and convex, so there is one minimum and calculus finds it. It penalises big misses more than small ones. And if the noise is Gaussian with constant variance, minimising squared error is exactly maximum likelihood estimation, so OLS is the statistically efficient answer under that assumption.

Setting the gradient to zero gives the solution. The gradient of L is -2 X^T (y - X w); setting it to zero gives the normal equations X^T X w = X^T y. Geometrically, the residual vector must be orthogonal to every column of X: the prediction is the projection of y onto the space the columns span. If the columns are linearly independent, X^T X is invertible and the solution is unique.

A worked example by hand

Take five points: x = 1, 2, 3, 4, 5 and y = 2, 4, 5, 4, 5. With one feature the solution simplifies to slope = Sxy / Sxx and intercept = mean(y) - slope x mean(x), where Sxy and Sxx are sums of centred products.

The means are 3 and 4. The centred x values are -2, -1, 0, 1, 2 and the centred y values are -2, 0, 1, 0, 1. Sxy = (-2)(-2) + (-1)(0) + (0)(1) + (1)(0) + (2)(1) = 6. Sxx = 4 + 1 + 0 + 1 + 4 = 10. The slope is 0.6 and the intercept is 4 - 0.6 x 3 = 2.2.

Predictions are 2.8, 3.4, 4.0, 4.6 and 5.2, so the residuals are -0.8, 0.6, 1.0, -0.6 and -0.2. They sum to zero. The sum of squared errors is 0.64 + 0.36 + 1.00 + 0.36 + 0.04 = 2.4. The total sum of squares around the mean is 4 + 0 + 1 + 0 + 1 = 6, so R squared = 1 - 2.4 / 6 = 0.6: the line explains 60 percent of the variance in y. Check it in code:

import numpy as np
x = np.array([1, 2, 3, 4, 5.0]); y = np.array([2, 4, 5, 4, 5.0])
X = np.column_stack([np.ones_like(x), x])
w, *_ = np.linalg.lstsq(X, y, rcond=None)
r = y - X @ w
print(w)                                   # [2.2 0.6]
print((r**2).sum(), 1 - (r**2).sum() / ((y - y.mean())**2).sum())   # 2.4 0.6

Solving it well: QR, SVD, Cholesky, gradient descent

The formula w = inv(X^T X) X^T y is correct algebra and poor numerics. Forming X^T X squares the condition number of the problem: if X has condition number 10^6, X^T X has 10^12, and in double precision you can lose about 12 of your 16 significant digits. Explicitly inverting is worse than solving. Three better routes exist.

• QR decomposition. Factor X = QR with Q orthonormal and R upper triangular, then solve R w = Q^T y by back substitution. It works on X directly, so conditioning is not squared. Cost is about 2np^2 operations. This is the safe default.
• SVD. Factor X = U S V^T and set w = V diag(1/s_i) U^T y, dropping singular values below a tolerance. Slower, but it handles rank-deficient X gracefully and the singular values tell you exactly which directions are poorly determined. np.linalg.lstsq uses an SVD-based LAPACK routine.
• Cholesky of the normal equations. Fastest, about np^2 operations to form X^T X plus p^3/3 to factor, and fine when features are standardised and not collinear. Use it when p is small, n is huge and you have checked conditioning.

When n or p is too large for any factorisation, or data arrives as a stream, use iterative methods. Gradient descent on the mean squared error converges at a rate governed by the condition number of X^T X, which is why standardising features first speeds it up dramatically; gradient descent in depth covers step sizes and the optimisation landscape explains why conditioning controls convergence.

def fit_gd(X, y, lr=0.1, steps=2000):
    mu, sd = X[:, 1:].mean(0), X[:, 1:].std(0)
    Xs = X.copy(); Xs[:, 1:] = (X[:, 1:] - mu) / sd      # keep intercept column
    w = np.zeros(X.shape[1]); n = len(y)
    for _ in range(steps):
        grad = -2 / n * Xs.T @ (y - Xs @ w)
        w -= lr * grad
    w[1:] /= sd; w[0] -= mu @ w[1:]                       # back to original units
    return w

Assumptions and diagnostics

Fitting needs no assumptions; interpreting does. The coefficients are unbiased estimates of the true linear effect if the errors have mean zero given X. Standard errors and p-values additionally assume errors with constant variance that are independent across examples, and exact small-sample intervals assume Gaussian errors. Each assumption has a diagnostic.

Collinearity deserves emphasis because it hides. Two features that move together, such as floor area and number of rooms, leave the prediction stable but make the individual coefficients swing wildly, even flip sign, between samples. The variance inflation factor for feature j is 1 / (1 - R_j^2), where R_j^2 comes from regressing feature j on the others; values above about 10 mean that coefficient should not be interpreted.

Ridge and lasso

Regularisation adds a penalty on the weights. It trades a little bias for less variance, fixes ill-conditioning, and with the right penalty selects features. Standardise features first, because the penalty is applied to raw coefficient sizes, and never penalise the intercept.

Ridge minimises ||y - X w||^2 + lambda ||w||^2. The closed form is w = (X^T X + lambda I)^-1 X^T y: adding lambda to the diagonal lifts every eigenvalue, so the system is always solvable and better conditioned. On the worked example with centred data, the slope becomes Sxy / (Sxx + lambda). With lambda = 2 that is 6 / 12 = 0.5, and the intercept becomes 4 - 0.5 x 3 = 2.5. The slope shrank towards zero; training error rose slightly; predictions on new data often improve when features are noisy or correlated.

Lasso minimises (1/2) ||y - X w||^2 + lambda ||w||_1. The absolute value has a corner at zero, so some weights land exactly at zero, which is feature selection. For one feature the solution is soft thresholding: slope = sign(Sxy) max(|Sxy| - lambda, 0) / Sxx. With lambda = 2, the slope is (6 - 2) / 10 = 0.4; with lambda = 6 or more it is exactly zero. With many features there is no closed form, and coordinate descent applies that one-feature update to each weight in turn:

def soft(z, t):
    return np.sign(z) * np.maximum(np.abs(z) - t, 0.0)

def lasso_cd(X, y, lam, iters=200):
    """X standardised, y centred; objective 0.5*||y - Xw||^2 + lam*||w||_1."""
    n, p = X.shape
    w = np.zeros(p); r = y.copy()
    col_sq = (X ** 2).sum(0)
    for _ in range(iters):
        for j in range(p):
            r += X[:, j] * w[j]                 # remove feature j's contribution
            w[j] = soft(X[:, j] @ r, lam) / col_sq[j]
            r -= X[:, j] * w[j]
    return w

Elastic net mixes both penalties and behaves better than lasso when correlated features come in groups, since lasso tends to pick one arbitrarily. Choose lambda by cross-validation over a log-spaced grid; scikit-learn's RidgeCV, LassoCV and ElasticNetCV do this, using the same objective scaling conventions documented for each estimator, so read them before comparing lambda values across libraries.

Failure modes

Most bad regressions fail in one of these ways.

• Leakage. A feature that is computed from the target, or from the future, gives a superb R squared and useless predictions. Audit how and when each feature is produced.
• Extrapolation. A line fitted on x from 1 to 5 says nothing reliable about x = 50. Record the training range and flag predictions outside it.
• Reading coefficients as causes. A coefficient is the association holding the other included features fixed. Omitted confounders bias it, and adding or removing a correlated feature changes it.
• Scaling mistakes. Regularising unstandardised features penalises them by unit choice, and fitting the scaler on all data before splitting leaks test information.
• R squared worship. Training R squared always rises with more features. Judge on held-out error, and on whether residual plots look like noise.
• Heavy-tailed targets. Prices, latencies and incomes are skewed; a few large values dominate squared loss. Model log(y), or use Huber loss, and transform predictions back carefully.

Trade-offs and where it fits

Linear regression's limits are also its value. It cannot model interactions or curvature you do not engineer into the features, and it will underfit problems where tree ensembles or neural networks shine. In exchange it trains in milliseconds, needs little data, extrapolates in a predictable (if wrong) way, gives coefficients a domain expert can argue with, and has exact uncertainty estimates when assumptions hold. A useful habit is to fit it first on every regression task: if a complex model cannot beat it on held-out data, the complex model is not learning anything real. It also reappears inside larger systems, as the final layer of a neural network, as linear probes on embeddings, and as the local model in explanation methods. Expected value in algorithms is a good companion for the probability behind its error estimates, and double descent shows what happens to least squares when p grows past n.

What to do next

1. Re-derive the five-point example by hand, then confirm it with the NumPy snippet.
2. On a real dataset, split train and test first, standardise using training statistics only, and fit OLS with np.linalg.lstsq or scikit-learn's LinearRegression.
3. Plot residuals against fitted values and against each feature; compute variance inflation factors.
4. Fit ridge and lasso across a log-spaced lambda grid with cross-validation and compare held-out error and which coefficients survive.
5. Implement the gradient descent and coordinate descent functions above and check they match the library answers to several decimal places.
6. Write down the training range of each feature and add an out-of-range warning to anything that serves predictions.