Linear regression

Ordinary least squares#

The model is linear in its parameters, which is what makes the closed form possible.

y = X * beta + epsilon

• y — response vector, shape (n, 1)
• X — design matrix, shape (n, p), first column all ones for the intercept
• beta — parameter vector, shape (p, 1)
• epsilon — error term, shape (n, 1), assumed mean zero
• n — number of observations, p — number of parameters including the intercept

OLS minimises the residual sum of squares:

RSS(beta) = || y - X*beta ||^2 = (y - X*beta)^T (y - X*beta)

Setting the gradient to zero yields the normal equations:

X^T X beta_hat = X^T y

and, when X^T X is invertible:

beta_hat = (X^T X)^{-1} X^T y

• X^T — transpose of the design matrix
• (X^T X)^{-1} — inverse of the p x p Gram matrix
• beta_hat — the estimated parameters

Geometrically, X * beta_hat is the orthogonal projection of y onto the column space of X. The residual vector is orthogonal to every column of X, which is the entire content of the normal equations.

Do not invert the matrix#

(X^T X)^{-1} is correct algebra and poor numerics. Forming X^T X squares the condition number of X:

kappa(X^T X) = kappa(X)^2

• kappa — condition number, the ratio of largest to smallest singular value

Production solvers use QR decomposition or SVD on X directly. scikit-learn's LinearRegression calls scipy.linalg.lstsq, which uses the SVD-based gelsd driver. When X is rank-deficient it returns the minimum-norm solution rather than raising an error.

Normal equations:  O(n p^2 + p^3)
Householder QR:    O(n p^2 - p^3/3)
SVD:               O(n p^2) with a larger constant, most robust
Gradient descent:  O(n p) per iteration, preferred when n or p is very large

Gauss-Markov#

The theorem requires four assumptions:

1. Linearity in the parameters.
2. Strict exogeneity, E[epsilon | X] = 0.
3. Homoscedasticity and no autocorrelation, Var(epsilon | X) = sigma^2 * I.
4. Full column rank of X.

Under these, beta_hat is the Best Linear Unbiased Estimator. It has minimum variance among all linear unbiased estimators.

Note what is not assumed: normality of errors. Normality is needed for exact t and F inference in finite samples, not for the Gauss-Markov result itself.

The variance of the estimator is:

Var(beta_hat) = sigma^2 * (X^T X)^{-1}

• sigma^2 — error variance, estimated by RSS / (n - p)

Note that "best unbiased" is a weaker claim than "best". Biased estimators can achieve lower mean squared error, which is exactly the opening that regularisation exploits.

Coefficient of determination#

R^2 = 1 - RSS / TSS,     TSS = SUM_i ( y_i - y_bar )^2

• RSS — residual sum of squares
• TSS — total sum of squares
• y_bar — mean of the observed response

R^2 never decreases when a predictor is added, including a column of pure noise. Adjusted R^2 penalises p:

R^2_adj = 1 - ( (1 - R^2) * (n - 1) / (n - p) )

On out-of-sample data R^2 can be negative, meaning the model is worse than predicting the training mean. That is a valid and informative result, not a bug.

Regularisation#

When p approaches or exceeds n, or predictors are collinear, OLS variance explodes. Adding a penalty trades bias for variance.

Ridge:  beta_hat = argmin  || y - X*beta ||^2 + lambda * || beta ||_2^2
Lasso:  beta_hat = argmin  || y - X*beta ||^2 + lambda * || beta ||_1

• lambda — regularisation strength, chosen by cross-validation
• || . ||_2^2 — squared L2 norm, the sum of squared coefficients
• || . ||_1 — L1 norm, the sum of absolute coefficients

Ridge has the closed form beta_hat = (X^T X + lambda*I)^{-1} X^T y. That matrix is always invertible for lambda > 0.

Lasso has no closed form. Its L1 penalty is non-differentiable at zero, and that is exactly why it drives coefficients to exact zero. The result is genuine variable selection. It is solved by coordinate descent or LARS.

Always standardise features before penalising. The penalty is not scale-invariant, so unstandardised features are penalised in proportion to their arbitrary units.

Diagnostics worth running#

For inference — standard errors, p-values, confidence intervals — use statsmodels rather than scikit-learn. scikit-learn is built for prediction and deliberately omits inferential output.

Key references#

• Legendre, A. M. (1805). Nouvelles méthodes pour la détermination des orbites des comètes. First publication of least squares.
• Gauss, C. F. (1809). Theoria Motus Corporum Coelestium. Claims prior use from 1795 and adds the probabilistic justification.
• Hoerl, A. & Kennard, R. (1970). Ridge Regression: Biased Estimation for Nonorthogonal Problems. Technometrics 12(1).
• Tibshirani, R. (1996). Regression Shrinkage and Selection via the Lasso. JRSS-B 58(1).
• Zou, H. & Hastie, T. (2005). Regularization and Variable Selection via the Elastic Net. JRSS-B 67(2).
• Golub, G. & Van Loan, C. (2013). Matrix Computations, 4th ed. Johns Hopkins. The reference for the numerical arguments above.

Current state#

Linear regression is not superseded, and treating it as a beginner's toy is a mistake. It remains the default in econometrics and epidemiology. It wins wherever a defensible coefficient matters more than a marginal accuracy gain.

Generalised linear models (Nelder & Wedderburn, 1972) extend it to non-Gaussian responses through a link function. That is where logistic and Poisson regression come from. Generalised additive models (Hastie & Tibshirani, 1986) relax linearity per feature while keeping additive interpretability.

It is also the correct baseline for any tabular regression problem. Report it alongside your gradient-boosted model. If the gap is small, the extra complexity is not earning its keep.

What to learn next#

• Logistic regression — the GLM with a logit link.
• Overfitting and underfitting — the bias-variance mechanics behind ridge and lasso.
• Optimization — gradient methods for when the closed form is too expensive.