Ridge Regression and L2 Regularization

Ordinary least squares (OLS) chooses weights w to minimize the residual sum of squares ``` L(w) = ‖y − Xw‖² ``` Ridge regression modifies this objective by adding a penalty proportional to the squared L2 norm of the weight vector: ``` L_ridge(w) = ‖y − Xw‖² + λ‖w‖² ``` Here λ ≥ 0 is the regularization strength, a hyperparameter that is not learned from the training loss directly — it is chosen by validation. When λ = 0, ridge regression reduces exactly to OLS. As λ → ∞, every coefficient is forced toward zero and the model collapses to predicting a constant. The penalty term ‖w‖² = w₁² + w₂² + ... + wₚ² is called the L2 penalty or Tikhonov regularization; the resulting shrinkage of coefficients toward zero is called ridge shrinkage. Crucially, the bias term (intercept) is conventionally excluded from the penalty — we want to shrink the slopes that describe feature effects, not force predictions toward zero regardless of the data's mean. In practice this is handled by centering y and each column of X before fitting, so the intercept can be recovered separately and the penalty applies only to the centered coefficients. The name "ridge" comes from the geometry of the solution: adding λI to XᵀX is literally adding a ridge along the diagonal of the matrix before inverting it, which pushes the matrix's eigenvalues away from zero and away from singularity. This single algebraic modification is what makes ridge regression well-defined even when OLS is not — for example, when the number of features p exceeds the number of examples n, or when features are perfectly collinear.

Step 1 — Center (and optionally scale) the data. Subtract the mean of y and the mean of each column of X. If features are on different scales, standardize each column to unit variance — otherwise the penalty λ‖w‖² unfairly punishes large-scale features more than small-scale ones, since a coefficient's contribution to ‖w‖² depends on the units of its feature.

Step 2 — Form the penalized normal equations. Differentiate the ridge loss with respect to w and set the gradient to zero: ``` ∇_w L_ridge = -2Xᵀ(y - Xw) + 2λw = 0 ⟹ XᵀXw + λw = Xᵀy ⟹ (XᵀX + λI)w = Xᵀy ``` Step 3 — Solve the linear system. ``` w* = (XᵀX + λI)⁻¹ Xᵀy ``` This is the ridge estimator. Unlike the OLS normal equations, (XᵀX + λI) is guaranteed invertible for any λ > 0, because adding λI shifts every eigenvalue of the positive semi-definite matrix XᵀX up by λ, making all eigenvalues strictly positive.

Step 4 — Select λ by cross-validation. Fit the model for a grid of λ values, evaluate held-out error (e.g., via k-fold cross-validation), and pick the λ that minimizes validation error. This step is essential: λ is a hyperparameter controlling the bias-variance tradeoff, and it cannot be tuned by minimizing training loss (training loss is always minimized at λ = 0). ``` function FIT_RIDGE(X, y, lambda):

# X: n-by-p design matrix (rows = examples), y: n-vector of targets x_mean <- column_means(X) y_mean <- mean(y) X_c <- X - x_mean # center features y_c <- y - y_mean # center target I <- identity_matrix(p) A <- transpose(X_c) X_c + lambda I # p-by-p, SPD for lambda > 0 b <- transpose(X_c) * y_c # p-vector w <- SOLVE_LINEAR_SYSTEM(A, b) # e.g. Cholesky decomposition of A intercept <- y_mean - dot(x_mean, w) return (intercept, w) function PREDICT_RIDGE(intercept, w, x_new): return intercept + dot(w, x_new) function CHOOSE_LAMBDA(X, y, lambda_grid, k): best_lambda <- None best_score <- +infinity for lambda in lambda_grid: score <- K_FOLD_CV_ERROR(X, y, lambda, k) # average squared error over folds if score < best_score: best_score <- score best_lambda <- lambda return best_lambda ``` The linear system `A w = b` is solved with Cholesky decomposition rather than by explicitly forming the matrix inverse, both for numerical stability and speed — A is symmetric positive definite whenever λ > 0, which is exactly the setting Cholesky is designed for.

Time complexity. Forming XᵀX costs O(np²) (multiplying an n×p matrix by its transpose). Adding λI is O(p²) and negligible. Solving the p×p system A w = b via Cholesky decomposition costs O(p³). Total training cost is ``` O(np² + p³) ``` which is identical to OLS via the normal equations — the ridge penalty adds no asymptotic overhead, only a diagonal shift. If p ≫ n, an equivalent dual formulation solves an n×n system instead, costing O(n²p + n³); this dual form is exactly the bridge to Kernels and the Kernel Trick, since it expresses the solution purely in terms of inner products XXᵀ. Prediction on a new point costs O(p) per query.

Space complexity. Storing X costs O(np); storing XᵀX (or the Cholesky factor) costs O(p²). Total space is O(np + p²).

Correctness — why (XᵀX + λI) is always invertible. XᵀX is symmetric and positive semi-definite: for any vector v, vᵀ(XᵀX)v = ‖Xv‖² ≥ 0, so all its eigenvalues are ≥ 0. If M has eigenvalue μ ≥ 0 with eigenvector v, then M + λI has eigenvalue μ + λ with the same eigenvector. Since μ ≥ 0 and λ > 0, every eigenvalue of XᵀX + λI is strictly positive, so the matrix is positive definite and therefore invertible. This holds regardless of whether X has full column rank — ridge regression is well-posed even when p > n or columns are exactly collinear, precisely the situations where plain OLS fails.

Correctness — the objective is strictly convex, so the closed form is the unique global minimum. The ridge loss L_ridge(w) = ‖y − Xw‖² + λ‖w‖² is a sum of a convex quadratic (the squared error, convex but only strictly convex along directions where X has rank) and a strictly convex quadratic λ‖w‖² (for λ > 0, strictly convex in every direction). The sum of a convex function and a strictly convex function is strictly convex, so L_ridge has a unique minimizer everywhere, and setting the gradient to zero (as derived above) finds it exactly.

Why shrinkage reduces variance — the bias-variance decomposition. Write the ridge estimator in terms of the eigendecomposition of XᵀX = QΛQᵀ, with eigenvalues d₁ ≥ ... ≥ dₚ ≥ 0. One can show ``` w_ridge = Q diag( d_j / (d_j + λ) ) Qᵀ w_OLS ``` Each OLS component along eigenvector qⱼ is shrunk by the factor dⱼ/(dⱼ + λ) ∈ (0, 1]. Directions with small eigenvalue dⱼ (high-variance, poorly determined directions, exactly where OLS is unstable) are shrunk the most, because dⱼ/(dⱼ + λ) → 0 as dⱼ → 0. Directions with large eigenvalue (well-determined by the data) are barely shrunk. This is exactly the desired behavior: ridge regression suppresses the noisy, ill-determined components of the OLS solution while preserving the well-determined ones, which lowers estimator variance far more than it increases bias whenever the design matrix is ill-conditioned.

Tradeoffs, stated honestly. Ridge never sets a coefficient to exactly zero (unlike Lasso's L1 penalty), so it does not perform feature selection — it shrinks but does not sparsify. Ridge also requires choosing λ, adding a hyperparameter search (though a single cheap trick, generalized cross-validation, can select λ without refitting from scratch for each candidate). And centering/scaling is not optional bookkeeping — skipping it silently distorts which features get penalized most.

Consider a tiny regression problem with one feature (already centered, so no intercept needed), n = 4 examples:

First, center x and y. Mean of x is 2.5, mean of y is 5.0.

Compute the sufficient statistics: ``` Sxx = Σ(xᶜ)² = 2.25 + 0.25 + 0.25 + 2.25 = 5.0 Sxy = Σ(xᶜ · yᶜ) = (−1.5)(−2.9) + (−0.5)(−1.1) + (0.5)(1.2) + (1.5)(2.8) = 4.35 + 0.55 + 0.60 + 4.20 = 9.70 ``` OLS solution (λ = 0): w_OLS = Sxy / Sxx = 9.70 / 5.0 = 1.94.

Ridge solution at λ = 1: ``` w_ridge = Sxy / (Sxx + λ) = 9.70 / (5.0 + 1) = 9.70 / 6.0 ≈ 1.6167 ``` Ridge solution at λ = 5: ``` w_ridge = 9.70 / (5.0 + 5) = 9.70 / 10.0 = 0.97 ``` Observe the shrinkage factor Sxx/(Sxx+λ) directly: at λ=1 it is 5/6 ≈ 0.833, at λ=5 it is 5/10 = 0.5. This matches the general formula dⱼ/(dⱼ+λ) derived above, with p=1 so d₁ = Sxx = 5.

As λ grows, the fitted slope shrinks toward zero and the intercept (which is unpenalized) drifts toward the unconditional mean of y, since the model relies less and less on x to explain y. In a real application, λ = 1 vs. λ = 5 vs. λ = 50 would be compared not by inspecting coefficients but by k-fold cross-validation error on held-out points — the "right" λ depends on how noisy the true relationship is, which the training data alone under-determines.

Ridge regression's mathematical machinery predates its name: adding a diagonal term to stabilize an ill-conditioned linear system is a special case of Tikhonov regularization, developed in the Soviet school of numerical analysis and inverse problems for stabilizing solutions to ill-posed integral equations. Statisticians independently rediscovered the same fix for collinear regression design matrices and gave it the name that stuck in machine learning.

Ridge regression (and its L2 penalty in general) is embedded throughout modern machine learning far beyond textbook linear regression. In neural networks, the identical penalty term appears as weight decay, added to the training loss to prevent weights from growing unboundedly and to control overfitting, especially important in deep networks with millions of parameters. In genomics and chemometrics, ridge regression is standard when the number of measured features (genes, spectral bands) vastly exceeds the number of samples, exactly the p ≫ n regime where OLS is singular. Support vector machines and kernel ridge regression use the L2-penalized objective directly, and the dual formulation of ridge regression — expressing the solution through the Gram matrix XXᵀ rather than XᵀX — is the algebraic seed that grows into Kernels and the Kernel Trick, allowing ridge regression to be lifted into high-dimensional or infinite-dimensional feature spaces without ever computing the features explicitly. Recommender systems and collaborative filtering use L2 penalties on latent factor matrices to prevent overfitting to sparse user-item interactions. And in Bayesian terms, ridge regression is exactly the maximum a posteriori estimate under a Gaussian prior on the weights, connecting it directly to Bayesian linear regression and to the broader theme of regularization-as-prior that recurs throughout probabilistic machine learning.

Ask the tutor: How does the shrinkage factor dⱼ/(dⱼ+λ) change if I use a different penalty, like λ‖w‖₁, instead of λ‖w‖²?