Regularised Regression

this page collects the closed-form solutions to regularised regression (where they exist) and the iterative approximations we fall back on (where they don’t). 𐃏 along the way we will see that regularisation is not an ad-hoc hack but a perfectly sensible artefact of estimation: it drops straight out of MAP (maximum a posteriori) inference once you put a prior on the coefficients.

motivation#

ordinary least squares is the best unbiased linear estimator — but nobody said unbiased was the place to be.

• when columns of \(X\) are correlated, or \(p\) approaches \(n\), the matrix \(X^\top X\) is near-singular: tiny wobbles in \(y\) produce wild swings in \(\hat\beta\). variance explodes while bias stays at zero — a terrible trade (see the bias-variance decomposition).
• regularisation deliberately buys a little bias to sell a lot of variance: shrink the coefficients toward zero, and the estimator stops chasing noise.
• a second, independent payoff: if the penalty can zero coefficients out exactly (the lasso can), you get variable selection and an interpretable model for free.

two pieces of housekeeping before any penalty touches anything (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009):

• standardise the inputs first. penalised solutions are not equivariant under feature scaling — a coefficient on “metres” and one on “millimetres” would be punished unequally.
• never penalise the intercept. shrinking the mean of \(y\) toward zero is not a belief anyone holds. centre \(y\) and the columns of \(X\), fit without an intercept, add the means back at the end.

ridge#

closed form#

ridge regression penalises the squared \(\ell_2\) norm:

\begin{equation} \hat\beta^{\mathrm{ridge}} = \arg\min_{\beta}\; \lVert y - X\beta \rVert_2^2 + \lambda \lVert \beta \rVert_2^2, \qquad \lambda \ge 0. \end{equation}

the objective is a quadratic in \(\beta\); setting the gradient \(-2X^\top(y - X\beta) + 2\lambda\beta\) to zero gives

\begin{equation} \boxed{\;\hat\beta^{\mathrm{ridge}} = (X^\top X + \lambda I)^{-1} X^\top y\;} \end{equation}

• the inverse always exists for \(\lambda > 0\): \(X^\top X\) is positive semi-definite, so every eigenvalue of \(X^\top X + \lambda I\) is at least \(\lambda\). this was the original selling point — ridge as a numerical fix for singular normal equations, the statistics arriving later.¹
• at \(\lambda = 0\) we recover OLS (ordinary least squares); as \(\lambda \to \infty\) the coefficients slide monotonically to zero. no coefficient ever becomes exactly zero at finite \(\lambda\) — shrinkage, not selection.

the svd view: shrinkage factors#

take the thin SVD \(X = U D V^\top\), with singular values \(d_1 \ge d_2 \ge \dots \ge d_p \ge 0\). substituting into the closed form, 𐃏

\begin{align*} X\hat\beta^{\mathrm{ridge}} &= X (X^\top X + \lambda I)^{-1} X^\top y \\ &= U D V^\top (V D^2 V^\top + \lambda V V^\top)^{-1} V D U^\top y \\ &= U D (D^2 + \lambda I)^{-1} D\, U^\top y \\ &= \sum_{j=1}^{p} u_j\, \frac{d_j^2}{d_j^2 + \lambda}\, u_j^\top y. \end{align*}

compare with OLS, which is \(\sum_j u_j\, u_j^\top y\): ridge computes the same coordinates of \(y\) in the basis \(\{u_j\}\), then shrinks the \(j\)-th coordinate by the factor

\begin{equation} s_j(\lambda) = \frac{d_j^2}{d_j^2 + \lambda} \in [0, 1). \end{equation}

• directions with large \(d_j\) (high-variance principal directions of the column space) are barely touched; directions with small \(d_j\) — precisely the ones where OLS variance blows up — are crushed (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).
• ridge therefore assumes the signal lives in the high-variance directions of the inputs. usually reasonable; occasionally spectacularly wrong.

effective degrees of freedom#

ridge is a linear smoother, \(\hat y = H_\lambda y\) with \(H_\lambda = X(X^\top X + \lambda I)^{-1}X^\top\). its effective degrees of freedom is the trace:

\begin{equation} \mathrm{df}(\lambda) = \operatorname{tr}(H_\lambda) = \sum_{j=1}^{p} \frac{d_j^2}{d_j^2 + \lambda}, \end{equation}

a smooth dial from \(\mathrm{df}(0) = p\) (full OLS) down to \(\mathrm{df}(\infty) = 0\) (constant model). it is the honest way to compare a ridge fit against a subset-selection fit “with \(k\) variables” — both are spending degrees of freedom, ridge just spends them in fractional coins (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

lasso#

swap the squared norm for an absolute one and the character of the problem changes completely:

\begin{equation} \hat\beta^{\mathrm{lasso}} = \arg\min_{\beta}\; \frac{1}{2n}\lVert y - X\beta \rVert_2^2 + \lambda \lVert \beta \rVert_1. \end{equation}

the objective is convex but no longer differentiable at \(\beta_j = 0\) — and that kink is exactly what manufactures sparsity.²

orthonormal design: soft-thresholding#

when \(X^\top X = I\) the problem separates across coordinates and each penalty has a closed form in terms of the OLS estimate \(\hat\beta_j\) (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009):

the lasso formula is the soft-thresholding operator \(S_\lambda(z) = \operatorname{sign}(z)(\lvert z\rvert - \lambda)_+\): translate every coefficient \(\lambda\) toward zero, and everything that reaches zero stays there. ridge never zeroes; hard thresholding zeroes discontinuously; soft thresholding zeroes continuously.

why corners give sparsity#

write each penalised problem in its equivalent constrained form: minimise the residual sum of squares subject to \(\lVert\beta\rVert_1 \le t\) (lasso) or \(\lVert\beta\rVert_2^2 \le t^2\) (ridge). the solution is the first point of the constraint set touched by the expanding elliptical contours of the RSS.

• the \(\ell_1\) ball is a cross-polytope: its corners sit on the coordinate axes, where some coefficients are exactly zero. an expanding ellipse has positive probability of first touching a corner (or an edge, zeroing several coefficients at once).
• the \(\ell_2\) ball is rotationally symmetric — no privileged directions, no corners, no zeros.
• in \(p\) dimensions the effect intensifies: the \(\ell_1\) ball is nearly all corner. 𐃏

coordinate descent#

no closed form in general, but the problem is tailor-made for coordinate descent (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009). fix all coordinates except \(\beta_j\) and define the partial residual \(r^{(j)} = y - \sum_{k \ne j} x_k \beta_k\). the one-dimensional problem in \(\beta_j\) is a quadratic plus an absolute value, minimised exactly by a soft-threshold:

\begin{equation} \beta_j \leftarrow \frac{S_{n\lambda}\!\big(x_j^\top r^{(j)}\big)}{\lVert x_j \rVert_2^2}. \end{equation}

• input standardised \(X\), centred \(y\), penalty \(\lambda\), initial \(\beta = 0\).
• loop cycle \(j = 1, \dots, p\): form the partial residual cheaply (add \(x_j\beta_j\) back into the running residual), soft-threshold, subtract the new contribution out.
• stop when no coordinate moves more than a tolerance.

convergence is guaranteed here because the non-smooth part of the objective is separable across coordinates — the one situation where coordinate descent on a non-differentiable convex function is safe. each sweep costs \(O(np)\). in practice one solves a whole grid of \(\lambda\) values from large to small, warm-starting each solve at the previous solution — the “pathwise” strategy that makes a full lasso path nearly as cheap as a single fit.

elastic net#

correlated features embarrass the lasso: among a group of near-duplicates it arbitrarily picks one and zeroes the rest, and the choice flips under resampling. the elastic net compromises:³

\begin{equation} \hat\beta = \arg\min_\beta\; \frac{1}{2n}\lVert y - X\beta\rVert_2^2 + \lambda\left(\alpha \lVert\beta\rVert_1 + \frac{1-\alpha}{2}\lVert\beta\rVert_2^2\right). \end{equation}

• \(\alpha = 1\) is the lasso, \(\alpha = 0\) is ridge; anything strictly between keeps the corners (sparsity survives) while the quadratic term makes the objective strictly convex — unique solutions even when \(p > n\), and a grouping effect: strongly correlated features receive similar coefficients instead of a winner-takes-all.
• coordinate descent barely changes; the update becomes

\begin{equation} \beta_j \leftarrow \frac{S_{n\lambda\alpha}\!\big(x_j^\top r^{(j)}\big)}{\lVert x_j\rVert_2^2 + n\lambda(1-\alpha)}, \end{equation}

a soft-threshold followed by a ridge-style proportional shrink.

the bayesian view: penalties are priors#

the cleanest justification for all of the above. model the data as gaussian, \(y \mid \beta \sim \mathcal{N}(X\beta, \sigma^2 I)\), and put a prior on \(\beta\). the MAP estimate maximises the log posterior:

\begin{equation} \hat\beta^{\mathrm{map}} = \arg\max_\beta\; \log p(y \mid \beta) + \log p(\beta). \end{equation}

• gaussian prior \(\beta \sim \mathcal{N}(0, \tau^2 I)\): the log prior is \(-\frac{1}{2\tau^2}\lVert\beta\rVert_2^2\) plus a constant, so

\begin{equation} \hat\beta^{\mathrm{map}} = \arg\min_\beta\; \lVert y - X\beta\rVert_2^2 + \frac{\sigma^2}{\tau^2}\lVert\beta\rVert_2^2, \end{equation}

which is ridge with \(\lambda = \sigma^2/\tau^2\) (Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon, 2020). tight prior (small \(\tau\)) means heavy shrinkage; the dial \(\lambda\) is literally your prior confidence that coefficients are small. for the gaussian–gaussian pair the posterior is itself gaussian, so the MAP estimate coincides with the posterior mean.

• laplace prior \(p(\beta_j) \propto \exp(-\lvert\beta_j\rvert/b)\): the log prior is \(-\frac{1}{b}\lVert\beta\rVert_1\) plus a constant, and the MAP estimate is the lasso. the laplace density’s spike at zero is the probabilistic shadow of the \(\ell_1\) ball’s corner.
• one honest caveat: the full laplace posterior puts zero mass on the event \(\beta_j = 0\) — sparsity is a property of the mode, not of the posterior. a bayesian who reports a posterior mean gets no zeros at all. 𐃏

in code#

ridge from scratch, with the svd check#

closed form, the SVD route, and the effective degrees of freedom — all on a synthetic problem with 3 real signals and 5 correlated decoys:

import numpy as np

rng = np.random.default_rng(42)
n, p = 60, 8
# correlated design: 3 real signals, 5 decoys correlated with them
Z = rng.normal(size=(n, 3))
X = np.hstack([Z, Z @ rng.normal(size=(3, 5)) * 0.5 + rng.normal(size=(n, 5))])
beta_true = np.array([4.0, -2.0, 1.0, 0, 0, 0, 0, 0])
y = X @ beta_true + rng.normal(0, 1.0, n)

# standardise columns, centre y (never penalise the intercept)
X = (X - X.mean(0)) / X.std(0)
y = y - y.mean()

def ridge(X, y, lam):
    p = X.shape[1]
    return np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

# svd view: identical fit via shrunken singular values
U, d, Vt = np.linalg.svd(X, full_matrices=False)
def ridge_svd(lam):
    return Vt.T @ ((d / (d**2 + lam)) * (U.T @ y))

lams = [0.0, 1.0, 10.0, 100.0, 1000.0]
print(f"{'lam':>7} {'df':>6}  coefficients (first 5)")
for lam in lams:
    b = ridge(X, y, lam)
    assert np.allclose(b, ridge_svd(lam)), "svd route disagrees"
    df = np.sum(d**2 / (d**2 + lam))
    print(f"{lam:>7.0f} {df:>6.2f}  " + " ".join(f"{v:>7.3f}" for v in b[:5]))

 lam     df  coefficients (first 5)
   0   8.00    3.211  -2.021   0.848   0.218   0.248
   1   7.74    3.118  -1.923   0.771   0.244   0.247
  10   6.20    2.542  -1.404   0.396   0.350   0.248
 100   2.61    1.095  -0.547   0.056   0.246   0.188
1000   0.44    0.183  -0.099   0.008   0.042   0.039

• the assert passes: \((X^\top X + \lambda I)^{-1}X^\top y\) and \(V(D^2+\lambda I)^{-1}D\,U^\top y\) are the same estimator.
• \(\mathrm{df}\) slides smoothly from \(8\) down to under half a degree of freedom.
• every coefficient shrinks, none dies: at \(\lambda = 1000\) the decoys are small but stubbornly nonzero. that is the coefficient path of a shrinker, not a selector.

lasso via coordinate descent#

same data, soft-thresholding in the loop:

def soft(z, g):                      # soft-threshold operator
    return np.sign(z) * np.maximum(np.abs(z) - g, 0.0)

def lasso_cd(X, y, lam, iters=500, tol=1e-8):
    n, p = X.shape
    b = np.zeros(p)
    col_ss = (X**2).sum(0)           # ||x_j||^2
    r = y.copy()                     # residual y - Xb
    for _ in range(iters):
        b_old = b.copy()
        for j in range(p):
            r += X[:, j] * b[j]      # remove j's contribution
            b[j] = soft(X[:, j] @ r, n * lam) / col_ss[j]
            r -= X[:, j] * b[j]      # add it back
        if np.max(np.abs(b - b_old)) < tol:
            break
    return b

print(f"{'lam':>6} {'nonzero':>8}  coefficients")
for lam in [0.001, 0.05, 0.2, 0.5, 1.0]:
    b = lasso_cd(X, y, lam)
    nz = int((np.abs(b) > 1e-10).sum())
    print(f"{lam:>6.3f} {nz:>8}  " + " ".join(f"{v:>7.3f}" for v in b))

  lam  nonzero  coefficients
0.001        8    3.209  -2.018   0.846   0.217   0.247   0.319  -0.138  -0.030
0.050        7    3.131  -1.879   0.743   0.164   0.188   0.301  -0.000  -0.032
0.200        5    3.073  -1.603   0.554   0.000   0.000   0.292   0.000  -0.027
0.500        5    2.768  -1.226   0.132   0.000   0.000   0.221   0.000  -0.021
1.000        3    2.302  -0.814   0.000   0.000   0.000   0.003   0.000  -0.000

• exact zeros appear and multiply as \(\lambda\) grows — the corners at work.
• the decoys die first; the two strong signals are the last coefficients standing (the weak third signal is eventually mistaken for a decoy — sparsity is not free, it is bias by another name).
• note the survivors are also shrunk: soft thresholding both selects and biases. this is the lasso’s standing critique, and the usual patch is to refit OLS on the selected support.

see also#

• ordinary least squares — the \(\lambda = 0\) endpoint of everything here
• bias-variance decomposition — the currency this page trades in
• cross validation — how \(\lambda\) actually gets chosen
• kernel methods — ridge in feature space is kernel ridge regression
• logistic regression — the same penalties bolt straight onto the logistic loss

References#

Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon (2020). Mathematics for Machine Learning, Cambridge University Press.

Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome (2009). The Elements of Statistical Learning, Springer.