What it actually does
Ordinary least squares will happily hand two nearly identical features enormous coefficients of opposite sign, because only their sum is pinned down by the data. Put a price on coefficient size and that stops. An L2 price (Ridge) is smooth, so it shrinks correlated features toward each other and keeps all of them; an L1 price (Lasso) has a corner at zero, so it pushes most coefficients exactly to zero and hands back a subset. Elastic net charges both. In every case you are buying a large reduction in variance with a small amount of bias — deliberately fitting the training data worse in order to predict new data better.
- Everything ordinary least squares assumes about linearity and residual structure still holds — the penalty changes the estimator, not the model
- Features are standardized, because the penalty is charged in feature units and is not scale-invariant
- For Lasso to recover the correct support, the true model must actually be sparse and the irrelevant features must not be too correlated with the relevant ones (the irrepresentable condition)
- The induced bias is acceptable — this is a prediction estimator by construction, and its coefficients are shrunk on purpose
Loss function
Squared error, plus a fee for every unit of coefficient you spend. Lambda sets the price; alpha decides whether the fee is charged on absolute size or on squared size. The two penalties look almost identical written down and behave completely differently, and the difference is entirely the kink: an absolute value is not differentiable at zero, and that non-differentiability is what produces exact zeros rather than merely small numbers.
- penalty strength — how expensive a unit of coefficient is
- mixing parameter: 1 is pure Lasso, 0 is pure Ridge, in between is elastic net
- sum of absolute coefficients — the penalty with a corner at zero
- sum of squared coefficients — smooth everywhere, so it never reaches zero
- number of training examples; the scaling keeps lambda comparable across dataset sizes
How it is minimized
Cyclic coordinate descent with soft-thresholding; Ridge alone also has a closed form
- partial residual: the residual with feature j’s own contribution added back in
- the soft-threshold — move z toward zero by gamma, and stop at zero
- column j of the design matrix
- the L2 share, which appears in the denominator: Ridge shrinks by dividing, Lasso by subtracting
The L1 term is convex but not differentiable at zero, and zero is exactly where the interesting behaviour lives — so a gradient method has nothing to work with. Coordinate descent sidesteps it: hold every other coefficient fixed and the remaining one-dimensional problem has an exact solution in closed form, kink included. That is the soft-threshold above, and reading it tells you the whole story — the numerator subtracts and clips at zero (selection), the denominator divides (shrinkage). Ridge alone is smooth and has the closed form (X-transpose-X + n·lambda·I)^-1 X-transpose-y, which is the normal equations with lambda added to the diagonal: the numerical fix for a near-singular matrix and the statistical fix for unstable coefficients are literally the same operation.
lambda
Penalty strength. The one hyperparameter that matters, and it must be cross-validated, not guessed (a geometric path of ~100 points from lambda_max (the smallest value that zeroes everything) down to 1e-3 x lambda_max)
alpha (l1_ratio)
L1 share of the penalty: selection versus shrinkage (0 (Ridge), 1 (Lasso), 0.5 to 0.9 for elastic net on correlated features)
standardization
Not optional. The penalty is charged in feature units, so unscaled features get arbitrary amounts of shrinkage
max sweeps / tolerance
Coordinate-descent stopping rule; convergence is fast, so these rarely bind (100 sweeps, tol 1e-7 on the largest coefficient change)
The objective is convex, so coordinate descent reaches the global optimum, and because each coordinate step is exact rather than approximate it converges quickly and monotonically. The failure modes are about identification, not convergence. With a group of correlated features, Lasso keeps one essentially arbitrarily and zeroes the rest: predictions barely move, but the selected support changes from one bootstrap resample to the next — which makes it a good predictor and a bad answer to "which variables matter". Elastic net blunts this and does not remove it. The second failure is silent: the penalty is not scale-invariant, so if features are not standardized the model is quietly encoding your choice of units as a modelling decision.
Coordinate descent with residual updates: O(nd) per sweep. With a precomputed Gram matrix: O(d^2) per sweep independent of n, which is the right trade when n >> d and the wrong one when d >> n. Ridge closed form: O(nd^2 + d^3). A full 100-point regularization path with warm starts typically costs only a few times a single fit.
How, where, and why it is used
Time-Series Forecasting
ViableThe same lag-and-calendar design matrix OLS would use, but built deliberately wide — every lag out to the seasonal period, every rolling statistic, every exogenous driver — and then the penalty decides which survive. This inverts the usual workflow: instead of choosing lags by squinting at an ACF plot, you offer all of them and read back the support.
- Demand forecasting with many correlated drivers (price, promotion, weather, competitor activity)
- Load and energy forecasting, where calendar and weather features are numerous and heavily overlapping
- Automatic lag selection on high-frequency series where hand-picking lags does not scale across thousands of series
- The regularized baseline in a forecasting bake-off, chosen because its coefficients stay stable across refits
Lag features are collinear almost by construction — lag 1 and lag 2 of a smooth series are nearly the same column — and that is precisely where OLS coefficients blow up. Ridge stabilizes them; Lasso additionally returns a short, defensible list of drivers. It is the wrong choice when the dynamics themselves are the story: ARIMA and state-space models encode autocorrelation structurally, whereas this only ever sees the lags you thought to include.
- Standardize before penalizing, using statistics computed only on the training window
- Fourier terms for seasonality rather than one dummy per period, so the penalty is spent on a handful of harmonics
- Difference or detrend first, or a large share of the penalty budget goes on fighting the trend
- Group the dummies of one categorical so the penalty keeps or drops the category as a unit
Rolling-origin cross-validation to pick lambda — never shuffled k-fold, which lets the model see the future. Score MASE against seasonal-naive, and separately check that the selected support is stable across folds: a support that changes completely between folds means the selection is noise, even when the accuracy looks fine.
- Choosing lambda by in-sample fit, which always prefers lambda = 0 and defeats the entire method
- Computing standardization statistics over the whole series before splitting — the most common leak in this design
- Reading the Lasso support as the set of causal drivers; it is a prediction subset chosen under correlation
- Large lambda collapses the model toward the intercept, which looks deceptively well-behaved on stable series
Anomaly Detection
AdaptedFit on a clean window and score residuals, as with OLS — but the penalty is what makes it workable when the monitored signals are many and redundant, which in telemetry they always are. The sharper variant regresses each signal on all the others and treats the Lasso support as a learned dependency graph, then alarms when a residual breaks a relationship that graph asserts.
- Multivariate process monitoring in manufacturing, where dozens of sensors measure overlapping physics
- Infrastructure telemetry where signals are redundant by design (replicas, mirrored feeds)
- Reconciliation across correlated financial series, where a total should be a fixed combination of parts
OLS residuals on collinear sensors are unstable: the fit can shuffle weight between two near-duplicate sensors between refits, and the residual moves with it, generating alerts that are artefacts of the estimator rather than events in the world. Ridge pins the residual down. It remains the wrong tool when normal behaviour is nonlinear or multi-modal — forcing a linear notion of normal onto that produces a detector that fires constantly and gets muted.
- Standardize on a confirmed-clean window and persist those statistics with the model
- Robust scaling, since the fit is sensitive to exactly the outliers being hunted
- Exclude the target sensor from its own predictor set, or the model learns the identity map
Precision@k against confirmed incidents, with the alarm threshold set from residual quantiles over a known-clean period and lambda cross-validated on clean data only.
- A contaminated training window teaches the model to treat the anomalies as normal
- Large lambda shrinks predictions toward the mean, inflating every residual uniformly and invalidating a threshold calibrated at a smaller lambda
- Lasso dropping a sensor entirely makes anomalies confined to that sensor invisible — sparsity and coverage are in direct tension here
Optimization
ViableThe Lasso is itself a constrained convex program — minimize squared error subject to an L1 ball of radius t — and the objective above is its Lagrangian form. What earns it a place under this heading is that the constraint is non-smooth, so ordinary gradient machinery fails and you have to reach for proximal or coordinate methods. Downstream, the sparse model it returns is also what makes a fitted surrogate cheap enough to embed inside a larger optimizer.
- Compressed sensing and sparse signal recovery, where the sparse solution is the answer rather than a convenience
- Sparse surrogate models embedded in the objective of an MPC or LP
- Feature-budget selection when the cost of scoring a feature at serving time is a hard constraint
This is the smallest problem where the geometry of the constraint set changes the answer. The L1 ball has corners on the axes; a convex objective touching a cornered set generically touches it at a corner, which is precisely why Lasso produces exact zeros and Ridge — whose ball is round — never does. Every later argument about non-smooth optimization is a variation on that picture.
- Standardize so the constraint set is isotropic; an unscaled feature is effectively given a larger budget
- Use per-coefficient penalty weights when features genuinely differ in acquisition cost, rather than a single lambda
Check the KKT conditions at the returned solution: every zeroed coefficient must have a correlation with the residual inside the threshold. That is an exact, cheap correctness test — unlike "the objective stopped moving", which a coordinate method can satisfy while still far from optimal on a coordinate it has not revisited.
- Declaring convergence from a stalled objective without checking the KKT conditions on the inactive set
- Proximal gradient with a step size above 1/L, which diverges quietly rather than failing loudly
- Optimizing against a fitted sparse surrogate outside the range where it was estimated
Risk & Fraud
Primary fitAttach the same penalties to logistic log loss instead of squared error. L2-penalized logistic regression is the default credit and fraud scorecard, and L1 is how a candidate pool of several thousand engineered features is cut down to the twenty a model risk committee will actually review.
- Credit scorecards and probability-of-default models
- Card-present and card-not-present transaction fraud scoring under heavy class imbalance
- Insurance underwriting and anti-money-laundering triage models
Regulated scoring wants a low-dimensional, inspectable model whose coefficients can be defended one at a time and which stays stable across quarterly refits. The penalty delivers both directly. Gradient boosting usually wins on AUC and routinely loses the deployment argument, because there is no comparable review story — which is a governance fact, not a modelling one, and it decides the choice anyway.
- Weight-of-evidence or monotone binned encodings, so each coefficient has a sign a reviewer can sanity-check
- Standardize before penalizing; leave the intercept unpenalized
- Group the dummies of one categorical under a single penalty so the variable is kept or dropped as a unit
Cross-validated AUC and KS for ranking, plus a calibration check, with lambda selected on out-of-time folds rather than random ones — the deployment gap in this domain is temporal, and random folds systematically over-state performance.
- Penalizing the intercept, which biases the base rate and quietly mis-calibrates every score
- Class imbalance interacting with the penalty so the model shrinks toward predicting the majority class
- Treating the selected support as a stable feature list when correlated features swap places between refits
Causal Inference & Uplift
AdaptedLasso is a natural way to choose which of many candidate confounders to adjust for — and the naive version of that is exactly the error that motivated modern causal ML. Regularizing the treatment coefficient shrinks the estimate toward zero, and selecting confounders from the outcome equation alone omits variables that matter mainly through treatment assignment. The correct form is post-double-selection, or the orthogonalized residual-on-residual construction that double machine learning uses.
- High-dimensional confounder selection in observational studies with more covariates than units
- The nuisance-function fits inside double machine learning and related orthogonal estimators
- Post-double-selection inference on panel data with many fixed effects
It is the honest answer to "I have 500 covariates and 2,000 units" — but only inside a design that keeps the penalty away from the parameter you care about. Used naively it yields intervals with essentially no coverage, which is strictly worse than a wide interval, because a confidently wrong estimate is acted upon.
- Never penalize the treatment indicator; it enters the model unshrunk by construction
- Standardize covariates but not the treatment, so the penalty does not depend on treatment coding
- Take the union of the variables selected in the treatment and outcome equations (double selection)
Interval coverage under simulation with a known effect, not predictive fit — a Lasso that predicts the outcome better may have selected away the confounders that mattered. Follow with sensitivity analysis for the variables the selection dropped.
- Reporting a shrunk treatment coefficient as an effect estimate
- Post-selection inference that ignores that the model itself was chosen from the data
- Assuming sparsity when the true confounding is dense — the guarantees are conditional on a sparsity that nothing verifies
What it costs to run
Effectively free at a single lambda, and a full 100-point path with warm starts typically costs only a few times one fit — which is why cross-validating lambda is cheap enough that there is no excuse for hand-picking it.
A dot product, and after Lasso a shorter one: the selected support is often a tenth of the candidate features. That matters less for the multiplications than for the feature fetches it removes from the serving path.
Cheap enough to refit continuously. The real constraint is not compute but that the selected support changes between refits, so anything downstream that keyed on "the model uses these twenty features" breaks quietly.
- Track the selected support across refits — churn in the support is the earliest signal that features are correlated and the selection is noise
- Watch the cross-validated lambda itself; a lambda drifting upward means the signal-to-noise in the features is deteriorating
- Monitor the standardization statistics separately from the model, since the penalty is scale-dependent and a scaler drift silently changes how much each feature is shrunk
- Compare coefficient magnitudes against the previous refit; under a fixed lambda a large jump means the conditioning of the design changed
- The standardization is part of the model, not preprocessing — persist the exact means and scales, and never refit them at inference
- Do not penalize the intercept; established libraries handle this, hand-rolled implementations routinely do not
- Lasso coefficients are biased toward zero by construction, so they are shrunk estimates and not marginal effects; reporting them as effects is a category error
- A feature Lasso dropped still has to exist in the serving pipeline if lambda is ever re-tuned — deleting it upstream makes the next refit silently different
Pros and cons, in context
Trades a controlled amount of bias for a large reduction in variance
Decisive when d approaches or exceeds n, or when features are collinear — the regime where OLS has no unique solution at all. Close to pointless when n >> d and features are near-orthogonal, where OLS is already the best unbiased estimator and the penalty only adds bias.
Lasso performs feature selection inside the fit rather than as a separate step
One convex problem instead of a greedy stepwise search, and the result is reproducible. But the support is unstable under correlated features, which makes it an excellent prediction tool and a poor answer to "which variables matter".
Ridge has a closed form and stays defined even when X-transpose-X is singular
Adding lambda to the diagonal is simultaneously the numerical fix for an ill-conditioned matrix and the statistical fix for unstable coefficients. That the two are the same operation is the most useful thing in the method.
The whole regularization path is cheap, so lambda can be cross-validated honestly
Warm starts make a hundred fits cost a few times one. Since there is no reason to guess lambda, a guessed lambda is a reliable sign that the rest of the pipeline was not validated either.
Coefficients are biased and cannot be read as effect estimates
Fatal in causal or regulatory work where the coefficient is the deliverable. It is exactly why double machine learning orthogonalizes rather than simply regularizing — the penalty that helps prediction is the same penalty that destroys the estimate.
Lasso chooses essentially arbitrarily among correlated features
Two near-duplicate columns: one is kept, one is zeroed, and which is which can change with the seed of your cross-validation split. Elastic net blunts this by keeping correlated groups together, but does not make the selection identifiable.
Requires standardization, which becomes state that ships with the model
Not a modelling flaw, but the flaw that actually bites in production: predictions drift because a scaler was refit somewhere upstream, and nothing in the model surfaces it.
Still a linear model
Regularization buys stability, never expressiveness. If the ceiling is that interactions dominate, no amount of penalty tuning moves it — that is what trees and boosting are for.
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.
"""Elastic-net coordinate descent - the objective, transcribed.
Each sweep visits one coefficient at a time and solves the resulting
one-dimensional problem exactly. The partial residual is the residual with
feature j's own contribution added back in, which is what turns the univariate
solution into a soft-threshold rather than a full re-fit.
"""
def soft_threshold(z, gamma):
# S(z, gamma) = sign(z) * max(|z| - gamma, 0)
if z > gamma:
return z - gamma
if z < -gamma:
return z + gamma
return 0.0
def fit(X, y, lam=0.1, alpha=1.0, sweeps=100):
"""alpha=1 is Lasso, alpha=0 is Ridge, in between is elastic net."""
n = len(X)
d = len(X[0])
beta = [0.0] * d
# norms[j] = (1/n) * sum_i x_ij^2 - the curvature of the 1-D subproblem
norms = [0.0] * d
for j in range(d):
total = 0.0
for i in range(n):
total += X[i][j] * X[i][j]
norms[j] = total / n
for _ in range(sweeps):
for j in range(d):
# rho_j = (1/n) sum_i x_ij * (y_i - sum_{k != j} x_ik beta_k)
rho = 0.0
for i in range(n):
partial = 0.0
for k in range(d):
if k != j:
partial += X[i][k] * beta[k]
rho += X[i][j] * (y[i] - partial)
rho /= n
# numerator subtracts and clips (selection); denominator divides
# (shrinkage). Both halves of the update are visible here.
beta[j] = soft_threshold(rho, lam * alpha) / (norms[j] + lam * (1.0 - alpha))
return betaO(sweeps * d^2 * n) — the partial residual is recomputed from scratch for every coordinate, which is the obvious thing and the wrong thing.
How this language actually wants it written.
"""Elastic net by coordinate descent - typed, validated, residual-tracked."""
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 ElasticNetModel:
"""A fitted penalized linear model, in the caller's original feature units."""
coefficients: Vector
intercept: float
lam: float
alpha: float
@property
def support(self) -> Vector:
"""Indices of the features the penalty actually kept."""
return np.flatnonzero(self.coefficients)
def predict(self, X: Matrix) -> Vector:
if X.ndim != 2 or X.shape[1] != self.coefficients.size:
raise ValueError(
f"expected (n, {self.coefficients.size}) design matrix, got {X.shape}"
)
return X @ self.coefficients + self.intercept
def _soft_threshold(z: float, gamma: float) -> float:
return float(np.sign(z) * max(abs(z) - gamma, 0.0))
def fit(
X: Matrix,
y: Vector,
lam: float = 0.1,
alpha: float = 1.0,
sweeps: int = 100,
tol: float = 1e-7,
) -> ElasticNetModel:
"""Fit by cyclic coordinate 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 not 0.0 <= alpha <= 1.0:
raise ValueError(f"alpha must lie in [0, 1], got {alpha}")
if lam < 0.0:
raise ValueError(f"lam must be non-negative, got {lam}")
# Standardization is part of the model, not a preprocessing nicety the
# caller may skip: the penalty is charged in feature units.
centre = X.mean(axis=0)
scale = X.std(axis=0)
scale[scale == 0.0] = 1.0
Xs = (X - centre) / scale
y_mean = float(y.mean())
n, d = Xs.shape
beta = np.zeros(d, dtype=np.float64)
residual = y - y_mean # y - X beta, kept current in place
norms = np.einsum("ij,ij->j", Xs, Xs) / n
for _ in range(sweeps):
largest_step = 0.0
for j in range(d):
column = Xs[:, j]
# Add feature j back in, solve, then push the change through the
# residual - O(n) per coordinate instead of O(n*d).
rho = float(column @ residual) / n + norms[j] * beta[j]
updated = _soft_threshold(rho, lam * alpha) / (norms[j] + lam * (1.0 - alpha))
step = updated - beta[j]
if step == 0.0:
continue
residual -= step * column
beta[j] = updated
largest_step = max(largest_step, abs(step))
if largest_step < tol:
break
coefficients = beta / scale # back into the caller's units
return ElasticNetModel(
coefficients=coefficients,
intercept=y_mean - float(centre @ coefficients),
lam=lam,
alpha=alpha,
)What changed: Two changes, one algorithmic and one structural. The inner recomputation of the partial residual is replaced by an incrementally maintained residual vector, which drops a whole factor of d from every sweep — the largest single win in the progression, and it is a change of algorithm rather than of language. Structurally, standardization moves inside the fit (the penalty is meaningless without it), input is validated before any work happens, and the fitted model becomes an immutable dataclass that knows its own support instead of a bare coefficient array.
O(n*d) per sweep, with the per-coordinate work executed in BLAS rather than the interpreter.
Faster, and the justification for each change.
"""Elastic net over a full lambda path - covariance updates, warm starts."""
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
def fit_path(
X: Matrix,
y: Vector,
alpha: float = 1.0,
n_lambdas: int = 100,
eps: float = 1e-3,
sweeps: int = 100,
tol: float = 1e-7,
) -> tuple[Vector, Matrix]:
"""Fit the entire regularization path in one call.
Two structural changes carry the speedup. Covariance updates: the only
quantities coordinate descent needs are X^T X and X^T y, so both are formed
once as single BLAS calls and every later sweep reads from the d x d Gram
matrix without touching the data again - the cost stops depending on n.
Warm starts: neighbouring solutions along the path are close, so each lambda
begins from the previous solution and converges in a handful of sweeps.
"""
if X.ndim != 2 or X.shape[0] != y.shape[0]:
raise ValueError(f"shape mismatch: X={X.shape}, y={y.shape}")
design = np.ascontiguousarray(X, dtype=np.float64)
scale = design.std(axis=0)
scale[scale == 0.0] = 1.0
design = (design - design.mean(axis=0)) / scale
centred_y = np.ascontiguousarray(y, dtype=np.float64) - y.mean()
n, d = design.shape
gram = (design.T @ design) / n # one GEMM, reused by every lambda
correlation = (design.T @ centred_y) / n
# Path from lambda_max - the smallest penalty that zeroes every coefficient
# - geometrically down to eps * lambda_max.
lam_max = float(np.abs(correlation).max()) / max(alpha, 1e-3)
lambdas = np.geomspace(lam_max, lam_max * eps, n_lambdas)
beta = np.zeros(d, dtype=np.float64)
coefficients = np.empty((n_lambdas, d), dtype=np.float64) # allocated once
for path_index, lam in enumerate(lambdas):
l1, l2 = lam * alpha, lam * (1.0 - alpha)
for _ in range(sweeps):
largest_step = 0.0
for j in range(d):
# Partial correlation straight from the Gram row.
rho = correlation[j] - gram[j] @ beta + gram[j, j] * beta[j]
updated = np.sign(rho) * max(abs(rho) - l1, 0.0) / (gram[j, j] + l2)
step = updated - beta[j]
if step == 0.0:
continue
beta[j] = updated
largest_step = max(largest_step, abs(step))
if largest_step < tol:
break
coefficients[path_index] = beta # warm start carries into the next lambda
return lambdas, coefficients / scaleWhat changed: The data is touched twice in total, to form the Gram matrix and the correlation vector, and after that the solver never sees it again — every sweep at every lambda is O(d^2) on cached quantities rather than O(nd) over the samples. That is only a win when n is much larger than d, which is the judgement call the docstring states outright. On top of it, the whole lambda path is solved in one call with warm starts, so cross-validating the penalty costs a few times a single fit instead of a hundred times.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
X^T X and X^T y become two GEMM calls, and each coordinate update then reads a cached Gram row instead of walking n samples.
Trade-off: Forming the Gram costs O(n*d^2) once and keeps a d x d matrix resident. When d >> n — the exact regime Lasso is usually reached for — this is the wrong trade and the residual-update version wins.
Batch work to amortize interpreter overhead
Solving the whole path in one call with warm starts means each subsequent lambda converges in a few sweeps, so 100 fits cost a small multiple of one.
Trade-off: The caller pays for the entire path even when a single lambda is wanted, and warm starting couples the solutions — a numerical problem at one lambda propagates down the rest of the path.
Ensure contiguous memory layout and a single dtype
ascontiguousarray with an explicit float64 means the GEMM operates on the array directly rather than on an internal copy, and the Gram rows read sequentially in the inner loop.
Trade-off: Materializes a copy of the input when the caller passed a Fortran-ordered or float32 array, which for very wide X is measurable memory pressure.
O(n*d^2) once for the Gram, then O(d^2) per sweep independent of n. Illustrative, not a measured benchmark.