Regularization Theory

Regularization modifies the empirical risk minimization objective

by adding a term that penalizes complex parameter vectors $\theta$. The penalized (Tikhonov / structural) form is

where $R(\theta)$ is a penalty function (commonly a norm of $\theta$) and $\lambda \geq 0$ is the regularization strength. The constrained form instead solves

restricting the search to a norm-ball of radius $C$. These are not two different ideas — they are the same idea viewed through primal and dual eyes.

Lagrangian duality connects them. Introduce a Lagrange multiplier $\lambda \geq 0$ for the constraint $R(\theta) \leq C$; the KKT stationarity condition of the constrained problem is exactly the stationarity condition of the penalized problem with that same $\lambda$. For convex, differentiable loss and convex penalty $R$, strong duality holds (Slater's condition is satisfied whenever the constraint set has nonempty interior), so for every $C$ there exists a $\lambda(C) \geq 0$ such that $\hat{\theta}_\lambda = \hat{\theta}_C$, and this correspondence is monotone: larger $C$ (looser constraint) corresponds to smaller $\lambda$ (weaker penalty). This is why practitioners freely move between "constrain the norm" and "penalize the norm" language — they describe the same solution path traced by one scalar knob. The two canonical penalties are the squared L2 norm, $R(\theta) = \|\theta\|_2^2 = \sum_j \theta_j^2$, giving ridge regression, and the L1 norm, $R(\theta) = \|\theta\|_1 = \sum_j |\theta_j|$, giving the LASSO (Least Absolute Shrinkage and Selection Operator). Both are convex, so the penalized objective remains a convex optimization problem whenever the base loss is convex, preserving global optimality guarantees. The qualitative difference between them — L2 shrinks, L1 selects — comes entirely from the shape of their unit balls, not from any difference in convexity or smoothness of the loss term itself.

Deriving the ridge solution

Take linear regression with squared loss. The penalized objective is

Differentiating and setting the gradient to zero:

Because $X^\top X + \lambda I$ is symmetric positive definite for any $\lambda > 0$ (even if $X^\top X$ is singular, e.g. when features are collinear or $n < p$), the solution is unique:

This is the single biggest practical reason ridge regression is used even absent any belief in sparsity: it repairs the invertibility problem that breaks ordinary least squares whenever features are correlated or the design matrix is rank-deficient.

Deriving the LASSO update (coordinate descent)

The LASSO objective, $\|y - X\theta\|_2^2 + \lambda\|\theta\|_1$, has no closed form for the full vector because $\|\theta\|_1$ is not differentiable at $\theta_j = 0$. But it is separable coordinate-wise given the other coordinates fixed, which makes coordinate descent exact and efficient. Fixing all $\theta_k$ for $k \neq j$, the one-dimensional subproblem in $\theta_j$ reduces (after algebra using the partial residual $r_j = y - \sum_{k\neq j} X_k \theta_k$) to the soft-thresholding operator:

The soft-threshold function is the proximal operator of the L1 norm; it shrinks its argument toward zero and sets it exactly to zero once $|z| \leq \gamma$ — this exact zeroing is the mechanism of sparsity. ``` Algorithm: Coordinate Descent for LASSO Input: X (n×p design matrix, columns standardized), y (n-vector), λ ≥ 0, tolerance ε, max_iters Output: θ (p-vector), the LASSO solution 1. Initialize θ ← 0 # or warm start from previous λ 2. Precompute column norms: c_j ← X_j^T X_j for all j 3. repeat 4. θ_old ← θ 5. for j = 1 to p: 6. r_j ← y − X θ + X_j θ_j # partial residual excluding feature j 7. z_j ← X_j^T r_j 8. θ_j ← soft_threshold(z_j, λ / 2) / c_j 9. until ‖θ − θ_old‖ < ε or iters ≥ max_iters 10. return θ function soft_threshold(z, γ): if z > γ: return z − γ if z < −γ: return z + γ else: return 0 ``` Each outer sweep over all $p$ coordinates costs $O(np)$ (recomputing residuals incrementally, not from scratch each time), and coordinate descent for LASSO is known to converge because the objective is convex and each coordinate update solves its 1-D subproblem exactly — a block coordinate minimization scheme, which for convex separable-nonsmooth-plus-smooth objectives is guaranteed to converge to the global optimum.

Choosing λ

Neither ridge nor LASSO tells you what $\lambda$ should be — that is determined empirically. The standard procedure is k-fold cross-validation over a geometric grid of $\lambda$ values, selecting the $\lambda$ that minimizes held-out error (or the largest $\lambda$ within one standard error of the minimum, the "one-standard-error rule," which prefers simpler models when performance is statistically indistinguishable). ``` Algorithm: Selecting λ via k-fold Cross-Validation Input: data (X, y), candidate grid Λ = {λ_1 > λ_2 > ... > λ_m}, number of folds k Output: λ, the selected regularization strength 1. Split data into k folds F_1, ..., F_k 2. for each λ in Λ: 3. for i = 1 to k: 4. train on all folds except F_i, using penalty λ 5. evaluate validation error on F_i 6. CV_error(λ) ← average validation error over the k folds 7. λ ← argmin_λ CV_error(λ) 8. refit final model on all data using λ 9. return λ ``` The LASSO path — the sequence of solutions as $\lambda$ decreases from $\lambda_{\max}$ (at which all coefficients are zero) to 0 (ordinary least squares) — can be computed efficiently with warm starts: initialize the coordinate descent for the next, smaller $\lambda$ from the solution at the previous, larger $\lambda$, which is both a good initial guess and dramatically speeds convergence.

The figure shows why the two penalties behave so differently. The elliptical contours are level sets of the loss $\|y - X\theta\|_2^2$; the solution to the constrained problem is the point where the smallest such ellipse touches the constraint region. The L2 ball is smooth everywhere, so the tangency point generically lands at an interior point of an edge — no coordinate is forced to zero. The L1 diamond has sharp corners lying exactly on the coordinate axes; because ellipses are far more likely to first touch a corner than a smooth edge (corners are "attractive" in a measure-theoretic sense — a positive-probability set of ellipse orientations touches each corner, while touching a specific interior edge point is a measure-zero event), the LASSO solution frequently has some coordinates exactly zero. This is the entire geometric explanation for feature selection, foreshadowing the deeper connection to Feature Selection.

Ridge regression. Computing $(X^\top X + \lambda I)^{-1} X^\top y$ requires forming $X^\top X$ in $O(np^2)$ time and inverting the resulting $p \times p$ matrix in $O(p^3)$ (via Cholesky decomposition, since $X^\top X + \lambda I$ is symmetric positive definite). Total: $O(np^2 + p^3)$ time, $O(p^2)$ space for the matrix. This is the same asymptotic cost as ordinary least squares — regularization adds essentially no overhead here, only numerical stability, since $\lambda I$ guarantees the matrix being inverted is well-conditioned even when $X^\top X$ alone is singular.

LASSO via coordinate descent. Each full sweep over $p$ coordinates costs $O(np)$ if residuals are updated incrementally rather than recomputed from scratch (each single-coordinate update touches only one column of $X$, an $O(n)$ operation, times $p$ coordinates). Convergence to a fixed tolerance typically requires a small constant number of sweeps in practice (often $10$–$100$), though the worst-case iteration count depends on the condition number of the problem; the overall practical cost is commonly quoted as $O(np)$ per sweep, times sweeps needed. Because the LASSO objective is convex (sum of a convex quadratic loss and a convex, though nonsmooth, L1 penalty), coordinate descent — which exactly minimizes each coordinate's one-dimensional convex subproblem in turn — provably converges to the global minimum: this is a special case of the general theorem that block coordinate descent converges for objectives that are the sum of a smooth convex function and a separable (coordinate-wise decomposable), possibly nonsmooth, convex function.

Why regularization improves generalization — the bias-variance argument. Decompose the expected test error of an estimator $\hat{\theta}$ at a fixed input as

Ordinary least squares is the unique unbiased linear estimator with minimum variance among unbiased estimators (the Gauss–Markov theorem), but "minimum variance among unbiased estimators" does not mean minimum variance overall. Ridge regression deliberately introduces bias — $\mathbb{E}[\hat\theta_{\text{ridge}}] \neq \theta^*$ in general — in exchange for a reduction in variance that can be large enough to lower total mean squared error. Concretely, for orthogonal design ($X^\top X = I$), the ridge estimate of each coefficient is $\hat\theta_j / (1+\lambda)$, a pure shrinkage; its variance is $\sigma^2/(1+\lambda)^2$, strictly smaller than the OLS variance $\sigma^2$ for any $\lambda > 0$, while its squared bias grows from zero. There exists a range of small $\lambda > 0$ for which the variance reduction outweighs the bias increase, so the ridge estimator has strictly lower MSE than OLS — this is a classical, provable fact, not a heuristic, and it directly foreshadows The Bias-Complexity Tradeoff.

Tradeoffs to be honest about. Ridge never produces exact zeros (its penalty is smooth away from the origin, so no soft-thresholding effect occurs), meaning it does not perform feature selection and all $p$ features remain in the model, which can hurt interpretability. LASSO's sparsity is a genuine benefit for interpretability and computational cost at prediction time, but it has known weaknesses: with strongly correlated features it tends to arbitrarily select one and zero out the rest rather than sharing weight, and it can select at most $n$ features when $p > n$. Elastic Net, $\lambda_1 \|\theta\|_1 + \lambda_2 \|\theta\|_2^2$, is the standard fix, combining both penalties to get sparsity plus grouped shrinkage among correlated features.

Consider a tiny regression with one feature standardized to unit variance, target $y$, and $n = 5$ observations, where the true relationship is $y = 2x + \text{noise}$ but a second, spurious feature $x_2$ is highly noisy and uncorrelated with $y$ in the population but weakly correlated by chance in this small sample.

(Step 6's LASSO arithmetic is schematic to show the mechanism: with $z_2 = X_2^\top r_2 = 1.3$ and threshold $\gamma = 0.5$ applied on the normalized scale, $z_2/c_2 = 0.26 < \gamma/c_2$, so soft-thresholding sends it to exactly zero.) Reading the table: OLS assigns a nonzero, noise-inflated coefficient to the spurious feature ($0.26$) simply because of sampling coincidence. Ridge shrinks both coefficients toward zero proportionally, reducing but not eliminating the spurious signal — it is still present, just smaller ($0.217$). LASSO's soft-thresholding step exactly zeros out the spurious coefficient because its correlation with the residual falls below the threshold $\gamma$, while the true signal survives (shrunk from $1.92$ to $1.42$, but nonzero). This concretely demonstrates the sparsity mechanism: it is not that LASSO "decides" a feature is useless via some search procedure, it is the direct algebraic consequence of the soft-threshold operator's flat region around zero.

Andrey Tikhonov developed what is now called Tikhonov regularization in the 1940s–60s to solve ill-posed inverse problems — problems (arising in geophysics, imaging, and integral equations) where a unique, stable solution does not exist without additional constraints. His insight, that adding a penalty term restores well-posedness by making the solution operator continuous in the data, is mathematically identical to what became ridge regression when Arthur Hoerl and Robert Kennard introduced it independently into statistics in 1970 to fix multicollinearity in linear regression. Robert Tibshirani introduced the LASSO in 1996, explicitly motivated by the desire for a penalty that would perform variable selection automatically rather than merely shrinking coefficients — recognizing that the geometry of the L1 ball, rather than the L2 ball, was the key structural ingredient.

Ridge regression and its generalizations underpin production systems wherever features are correlated or numerous relative to sample size: genomics (where the number of genes vastly exceeds the number of patients), econometric forecasting, and any linear model built on engineered feature sets with redundancy. LASSO and Elastic Net are the default tools for automatic feature selection in high-dimensional settings — recommender systems selecting a sparse set of predictive signals from thousands of candidates, and biostatistics pipelines identifying a handful of biomarkers from tens of thousands of gene expression measurements. Regularized logistic regression is a staple of interpretable risk-scoring models in credit and healthcare, where sparsity is prized because a model with 12 nonzero coefficients can be audited by a human regulator in a way a model with 10,000 cannot. In deep learning, the L2 penalty reappears as weight decay, added to virtually every training recipe for large neural networks to control the effective capacity of massively overparameterized models; L1 penalties appear in structured pruning objectives that seek to zero out entire neurons or channels for model compression. Support Vector Machines are, at their core, an L2-regularized hinge-loss minimization — the "maximum margin" geometric story and the "penalized empirical risk" algebraic story are the same convex program viewed two ways, another instance of the constrained/penalized duality developed here. Any framework offering `alpha`, `C`, or `weight_decay` hyperparameters — scikit-learn, PyTorch, XGBoost's regularization terms — is exposing exactly the $\lambda$ of this article to the user, making regularization theory one of the most universally applied ideas in applied machine learning.

Ask the tutor: How does the choice between L1 and L2 regularization change when the loss function itself is non-convex, such as in deep neural networks?