What it actually does
Draw the straight line (or hyperplane) that comes closest to every point at once, where "closest" means the sum of squared vertical distances is as small as it can be. Every coefficient answers one question: if this feature goes up by one unit and nothing else moves, how much does the prediction change?
- The relationship between features and target is linear in the parameters (features themselves may be transformed)
- Residuals are independent — violated by default on time series, which is why autocorrelation must be checked
- Residual variance is roughly constant across the range of predictions (homoscedasticity)
- Features are not near-perfectly collinear
Loss function
Average the squared gap between what the model predicts and what actually happened. Squaring makes errors in both directions count, and makes large errors count disproportionately — which is a modelling choice, not a neutral one: it is what makes least squares sensitive to outliers.
- number of training examples
- feature vector for example i (with a leading 1 for the intercept)
- observed target for example i
- the coefficient vector being fit
How it is minimized
Closed-form least squares (QR factorization), or batch gradient descent
- the n x d design matrix
- learning rate, for the iterative form
- the gradient: features weighted by their residuals
J is convex and quadratic in theta, so setting the gradient to zero has a unique closed-form solution — there are no local minima to worry about. In practice you never form the inverse: solving via QR factorization is the same answer with roughly half the condition-number sensitivity. Gradient descent is the right choice only when n is too large for the design matrix to fit in memory, or when the model is being fit online.
learning rate (iterative only)
Step size; too large diverges, too small crawls (1e-4 to 1e-1, after feature scaling)
epochs (iterative only)
Number of full passes over the data (100 to 10,000)
fit_intercept
Whether to add a constant column; omit only if features are already centred
The closed form is exact — one shot, no tuning. Gradient descent converges at a rate governed by the condition number of X-transpose-X, so unscaled features are the usual cause of slow or oscillating convergence. The failure mode worth naming is collinearity: when two features are nearly duplicated, X-transpose-X is near-singular, the coefficients become enormous and unstable with opposite signs, and the fit still looks fine on training data. That is the problem Ridge exists to solve.
Closed form: O(nd^2 + d^3) time, O(nd) memory. Gradient descent: O(nd) per epoch, O(d) memory beyond the data.
How, where, and why it is used
Time-Series Forecasting
ViableTurn the series into a supervised table: the target is the value at t+h, and the features are lagged values, rolling statistics, calendar dummies, and any exogenous drivers. The model then treats forecasting as ordinary regression — which is exactly why the lag design matters more than the estimator does.
- Trend estimation and detrending as a preprocessing step before ARIMA
- Short-horizon demand forecasts where drivers (price, promotion, weather) dominate autocorrelation
- The interpretable baseline in a forecasting bake-off, where coefficients on each driver are themselves the deliverable
Not because it is the most accurate forecaster — it usually is not — but because it is the one whose coefficients you can hand to a planner and defend. When the question is "how much of the lift came from the promotion", a linear model answers it and a gradient-boosted one does not.
- Lags at the seasonal period, not just t-1
- Rolling mean and standard deviation over trailing windows
- Fourier terms for smooth seasonality instead of one dummy per period
- Difference the series when the trend is stochastic rather than deterministic
Rolling-origin backtesting with an expanding window, scored with MASE against a seasonal-naive baseline. Never shuffled k-fold — it lets the model see the future.
- Look-ahead leakage: rolling statistics computed over the full series before splitting
- Extrapolating a fitted linear trend far past the training range, where it becomes fiction
- Residual autocorrelation invalidates the standard errors even when the point forecasts are fine
Anomaly Detection
AdaptedFit the model to normal behaviour, then score residuals: a point whose actual value sits many robust standard deviations from its prediction is flagged. The model supplies the expectation; the anomaly score is how badly reality missed it.
- Sensor and telemetry monitoring where one signal is a near-linear function of others
- Financial reconciliation, where a total should be a fixed linear combination of its parts
Only appropriate when normal behaviour genuinely is linear in the features. It is fast, explainable, and gives a residual you can attribute to a specific input — but it will not model a curved or multi-modal notion of normal, and forcing it to will produce a detector that fires constantly.
- Robust scaling, since the fit itself is sensitive to the outliers you are hunting
- Fit on a confirmed-clean window, or the anomalies contaminate the definition of normal
Precision@k against confirmed incidents, with the threshold set from the residual quantiles of a known-clean period.
- Training on contaminated data teaches the model to accept the anomalies
- A single large outlier drags the fitted line toward itself and masks its own detection
Optimization
AdaptedLeast squares is itself an unconstrained convex optimization problem, and its closed form is what a solver would converge to. In applied pipelines it more often plays the upstream role: producing the cost or demand estimate that a downstream optimizer takes as a fixed input.
- Cost-curve fitting that feeds a linear program
- Response-surface estimation for capacity planning
Worth studying here because it is the smallest complete example of the pattern the whole category runs on — write the objective, take the gradient, set it to zero, and check whether the solution is stable. Ridge, Lasso, and the SVM are all this same argument with a constraint added.
- Scale features before any iterative solve; conditioning is the whole story
- Encode the decision variables so the coefficients carry the units a planner expects
Compare the closed-form solution against an iterative solver on the same data; agreement to numerical tolerance is the correctness check.
- Near-collinear features make the solution unstable without making the fit look bad
- Optimizing a fitted surrogate past the range where the fit was estimated
Causal Inference & Uplift
Primary fitWith a treatment indicator and a credible identification strategy, the coefficient on that indicator is the estimated treatment effect, and its standard error is the uncertainty in it. Prediction is not the goal — the coefficient is.
- Difference-in-differences estimation on panel data
- Regression adjustment for covariates in a randomized experiment
- Regression discontinuity around a policy threshold
The default in causal work precisely because the estimand is a coefficient with an interval. A more accurate predictor is not a better estimator here — accuracy and identification are different objectives, and this is the clearest case where optimizing the wrong one misleads.
- Include confounders the design requires, not the features that improve fit
- Interaction terms where the effect is expected to vary by subgroup
Placebo tests and negative-control outcomes, plus sensitivity analysis for unmeasured confounding. Not R-squared.
- Controlling for a collider or a post-treatment variable induces bias rather than removing it
- Reading a coefficient as causal when the design does not support it — the arithmetic is identical either way
What it costs to run
Effectively free. A closed-form fit on a million rows and fifty features finishes in well under a second on one core.
A single dot product — sub-microsecond, trivially vectorized, and small enough to run anywhere including inside a SQL query or an edge device.
Cheap enough to refit on every batch. The usual cadence is driven by how fast the relationship drifts, not by compute.
- Track coefficient stability across refits — a sign flip is a louder signal than a small accuracy drop
- Monitor residual mean by segment; sustained one-sided error means the linear form no longer holds
- Watch the condition number of the design matrix, since new collinearity destabilizes coefficients silently
- Feature scaling must be persisted with the model; refitting a scaler at inference silently changes predictions
- Unseen categorical levels produce a column mismatch — pin the encoder, do not re-derive it
- The model extrapolates without complaint, so guard the input range explicitly if that matters
Pros and cons, in context
Coefficients are directly interpretable as marginal effects
Decisive in regulated, causal, or stakeholder-facing work, where the explanation IS the deliverable. Nearly worthless when the only thing anyone acts on is the prediction itself.
Convex objective with an exact closed-form solution
Means no tuning, no seeds, and reproducible fits — which matters far more in audited pipelines than the last point of accuracy.
Extremely cheap to train and to serve
Lets you refit continuously and run inference in places a larger model cannot go. Irrelevant when accuracy dominates and you have GPUs to spare.
Well-understood statistical theory for uncertainty
Standard errors and prediction intervals come for free under the assumptions — but only under them, which is exactly why the assumptions above are worth checking.
Cannot represent interactions or curvature unless you build them in
A real ceiling on tabular problems where interactions dominate, which is where gradient boosting wins outright. Not a limitation at all if the relationship genuinely is close to linear.
Squared error makes it sensitive to outliers
A serious problem in anomaly-adjacent work, where the outliers are the subject. Addressed by Huber loss or robust regression rather than by accepting it.
Collinearity destabilizes coefficients without hurting the fit
Dangerous specifically because the training metrics look fine. If interpretation matters, this is the failure to check for first — and the reason Ridge exists.
Extrapolates confidently outside the training range
Fine for interpolation-heavy scoring; genuinely hazardous for long-horizon forecasting, where a fitted trend keeps rising forever.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Least squares by batch gradient descent - the objective, transcribed.
Every loop below maps to one symbol in J(theta): the inner loop is the dot
product x_i . theta, err is the residual, and grad accumulates X^T (X theta - y).
"""
def fit(X, y, lr=0.01, epochs=1_000):
n = len(X)
d = len(X[0])
theta = [0.0] * d
bias = 0.0
for _ in range(epochs):
grad = [0.0] * d
grad_bias = 0.0
for i in range(n):
# prediction: x_i . theta + b
pred = bias
for j in range(d):
pred += theta[j] * X[i][j]
err = pred - y[i] # residual for example i
for j in range(d):
grad[j] += err * X[i][j] # accumulate X^T r
grad_bias += err
# theta <- theta - lr * (2/n) * gradient
scale = lr * 2.0 / n
for j in range(d):
theta[j] -= scale * grad[j]
bias -= scale * grad_bias
return theta, biasO(n*d) per epoch in pure Python — roughly 100x slower than NumPy for the same arithmetic.
How this language actually wants it written.
"""Least squares by gradient descent - typed, validated, vectorized."""
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
@dataclass(frozen=True)
class LinearModel:
"""A fitted model. Frozen so a fitted object cannot be mutated in place."""
weights: Vector
bias: float
def predict(self, X: Matrix) -> Vector:
if X.ndim != 2 or X.shape[1] != self.weights.size:
raise ValueError(
f"expected (n, {self.weights.size}) design matrix, got {X.shape}"
)
return X @ self.weights + self.bias
def fit(X: Matrix, y: Vector, lr: float = 0.01, epochs: int = 1_000) -> LinearModel:
"""Fit by batch gradient descent. Raises ValueError on malformed input."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if X.shape[0] != y.shape[0]:
raise ValueError(f"X has {X.shape[0]} rows but y has {y.shape[0]}")
if X.shape[0] == 0:
raise ValueError("cannot fit on an empty dataset")
n, d = X.shape
weights = np.zeros(d, dtype=np.float64)
bias = 0.0
scale = lr * 2.0 / n
for _ in range(epochs):
residual = X @ weights + bias - y # the whole inner loop, vectorized
weights -= scale * (X.T @ residual)
bias -= scale * float(residual.sum())
return LinearModel(weights=weights, bias=bias)What changed: The three nested Python loops collapse into two matrix products, and the function now fails fast on malformed input instead of raising an opaque IndexError deep inside a loop. The frozen dataclass makes a fitted model an immutable value rather than a loose tuple whose element order you have to remember.
O(n*d) per epoch, executed in BLAS rather than the interpreter.
Faster, and the justification for each change.
"""Least squares in closed form - one QR solve, no iteration."""
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def fit_exact(X: Matrix, y: Vector) -> tuple[Vector, float]:
"""Solve min ||X theta - y||^2 exactly.
np.linalg.lstsq factorizes the design matrix (QR/SVD) rather than forming
and inverting X^T X. Same solution, roughly half the condition-number
sensitivity - the difference shows up precisely when features are
near-collinear, which is when you most need the answer to be trustworthy.
"""
if X.ndim != 2 or X.shape[0] != y.shape[0]:
raise ValueError(f"shape mismatch: X={X.shape}, y={y.shape}")
# Prepend the intercept column so it is solved for jointly, not separately.
design = np.column_stack([np.ones(X.shape[0], dtype=X.dtype), X])
coefficients, *_ = np.linalg.lstsq(design, y, rcond=None)
return coefficients[1:], float(coefficients[0])What changed: The iterative solve disappears entirely. The objective is convex and quadratic, so the minimum has a closed form — thousands of gradient steps were approximating an answer that one factorization computes exactly, with no learning rate to tune and no convergence to check.
Replace a closed-form solve with a numerically stabler factorization
lstsq uses QR/SVD on X directly instead of forming X^T X, which squares the condition number and loses roughly half the available precision.
Trade-off: QR costs about twice the flops of the normal equations, and is worth it every time except in a tight inner loop on well-conditioned data.
Eliminate Python-level loops over samples
The entire epoch loop is gone; the work becomes a single LAPACK call that is blocked and multithreaded internally.
Trade-off: Requires the full design matrix in memory — O(n*d) — so beyond that point the iterative version is the only option.
Ensure contiguous memory layout and a single dtype
column_stack produces one contiguous float64 array, so LAPACK operates on it without an internal copy.
Trade-off: Materializes a second array of the input size; for very wide X the copy is measurable.
O(n*d^2 + d^3) once, no iteration. Illustrative, not a measured benchmark.