Regularisation: Ridge, LASSO, and Priors

1. \(\frac{\lambda}{\lambda+1}\) gives \(\frac{100}{101}=0.9901\), \(\frac{10}{11}=0.9091\), \(\frac12=0.5\), \(\frac{0.01}{1.01}=0.0099\).

2. \(0.9901+0.9091+0.5+0.0099 = 2.409\) of 4 parameters. The model behaves as though it had about two and a half free parameters instead of four.

3. Unregularised \(\kappa=100/0.01=10{,}000\). With ridge, \(\kappa=\frac{101}{1.01}=100\) — a hundredfold improvement, and by Lecture 03 that is a hundredfold reduction in the gradient-descent iteration count. Regularisation and preconditioning really are one operation.

4. The \(\lambda=0.01\) direction, which retains only 1% of its OLS value. That is exactly right. A small eigenvalue means the data barely varies along \(\vq_i\), so \(\vq_i^\top\hat{\vw}_0 = \frac{\vq_i^\top\mX^\top\vy}{\lambda_i}\) divides by a tiny number and the OLS estimate in that direction is dominated by noise. Ridge suppresses the directions whose estimates were least trustworthy and leaves the well-determined ones essentially alone — which is why it reduces variance far more than it adds bias.

LASSO. Minimise \(g(w)=\tfrac12(w-z)^2+\alpha|w|\), strictly convex, so \(w^\star\) is the unique point with \(0\in\partial g(w^\star)=(w^\star-z)+\alpha\,\partial|w^\star|\).

Case \(w^\star>0\): \(\partial|w^\star|=\{1\}\), so \(w^\star=z-\alpha\), consistent only if \(z>\alpha\).
Case \(w^\star<0\): \(w^\star=z+\alpha\), consistent only if \(z<-\alpha\).
Case \(w^\star=0\): need \(0\in -z+\alpha[-1,1]\), i.e. \(z\in[-\alpha,\alpha]\).

The three cases are exhaustive and mutually exclusive, giving \(w^\star = \operatorname{sign}(z)\max(|z|-\alpha,0)\). ∎

Ridge. \(h(w)=\tfrac12(w-z)^2+\tfrac{\alpha}{2}w^2\) is smooth, so \(h'(w)=(w-z)+\alpha w=0\), giving \(w^\star=\frac{z}{1+\alpha}\). ∎

The distinction. Ridge applies a multiplicative shrink of \(\frac{1}{1+\alpha}\) — the same proportion to every coefficient, so a coefficient of 10 loses ten times as much as a coefficient of 1, and no coefficient ever reaches zero at finite \(\alpha\). LASSO applies an additive shrink of \(\alpha\) with a floor at zero — every coefficient loses the same absolute amount, so small ones are wiped out entirely while large ones are barely affected in relative terms.

So "LASSO shrinks more" is wrong in both directions. For a large coefficient (\(z=10\), \(\alpha=1\)): LASSO gives 9, ridge gives 5 — ridge shrinks more. For a small one (\(z=0.5\)): LASSO gives 0, ridge gives 0.33 — LASSO shrinks more. They are qualitatively different operations, not the same operation at different strengths. LASSO makes a discrete decision about each coefficient (in or out); ridge makes a continuous adjustment to all of them. That is why LASSO's answer is unstable under resampling (small perturbations flip discrete decisions) while ridge's is smooth.

Proof. The OLS objective on the augmented data is

because the squared norm of a stacked vector is the sum of the squared norms of its blocks. That is exactly the ridge objective, so the minimisers coincide. ∎

Confirming via the normal equations: \(\tilde{\mX}^\top\tilde{\mX} = \mX^\top\mX+\alpha\mI\) and \(\tilde{\mX}^\top\tilde{\vy}=\mX^\top\vy\), giving precisely the ridge system.

Ridge as data augmentation. The \(d\) extra rows are \(d\) fictitious observations, each saying "feature \(j\) alone, scaled by \(\sqrt\alpha\), produces output 0". They are the prior, written as data. This makes the strength of a prior directly comparable to an amount of evidence, which is a genuinely useful mental model: \(\alpha\) large means "I have a lot of pseudo-observations telling me the weights are small", and it explains why the prior's influence fades as \(n\) grows — the real rows come to outnumber the fake ones.

It also generalises immediately. Replacing \(\sqrt\alpha\mI\) with \(\sqrt\alpha\mathbf{D}\) implements a different penalty \(\alpha\|\mathbf{D}\vw\|^2\), which is how smoothing splines penalise curvature and how graph-Laplacian regularisers encourage neighbouring coefficients to agree.

Why it is numerically better. The condition number argument from Phase A Lecture 10: forming \(\mX^\top\mX\) squares \(\kappa\), so a design with \(\kappa(\mX)=10^{8}\) gives \(\kappa(\mX^\top\mX)=10^{16}\), exhausting float64 before \(\alpha\) is even added. The augmented matrix \(\tilde{\mX}\) has \(\kappa(\tilde{\mX})=\sqrt{\kappa(\mX^\top\mX+\alpha\mI)}\) — the square root, not the square. Running lstsq on \(\tilde{\mX}\) therefore loses about half as many digits.

This is what production ridge implementations do (scikit-learn's solver='lsqr' and 'sparse_cg' work on the augmented system), and it is also the trick that makes ridge feasible for sparse \(\mX\), where forming the dense Gram matrix would be impossible.

import numpy as np
rng = np.random.default_rng(0)

def make(n=100, d=200, seed=None):
    r = np.random.default_rng(seed)
    g1, g2 = r.normal(size=(n,1)), r.normal(size=(n,1))
    X = r.normal(size=(n,d))
    X[:, :5]   = 0.9*g1 + 0.44*X[:, :5]        # correlated group A
    X[:, 5:10] = 0.9*g2 + 0.44*X[:, 5:10]      # correlated group B
    w = np.zeros(d); w[:10] = 1.5
    return X, X@w + r.normal(size=n), w

def enet(Z, y, a, rho, T=300):
    n_, d_ = Z.shape
    w = np.zeros(d_); G = Z.T@Z; b = Z.T@y; dg = np.diag(G).copy()
    for _ in range(T):
        for j in range(d_):
            r_ = b[j] - G[j]@w + G[j,j]*w[j]
            z = r_/(dg[j] + a*n_*(1-rho))
            t = a*n_*rho/(dg[j] + a*n_*(1-rho))
            w[j] = np.sign(z)*max(abs(z)-t, 0.0)
    return w

methods = {"ridge": lambda Z,y,a: enet(Z,y,a,0.0),
           "lasso": lambda Z,y,a: enet(Z,y,a,1.0),
           "enet":  lambda Z,y,a: enet(Z,y,a,0.5)}

sel = {k: np.zeros(200) for k in methods}
stats = {k: [] for k in methods}
for rep in range(50):
    X, y, wt = make(seed=rep); Xte, yte, _ = make(seed=1000+rep)
    Z = (X - X.mean(0))/(X.std(0)+1e-12); yc = y - y.mean()
    Zt = (Xte - X.mean(0))/(X.std(0)+1e-12)
    for k, f in methods.items():
        best = min(((np.mean((Zt@f(Z,yc,a) - (yte-y.mean()))**2), a)
                    for a in np.logspace(-3, 0, 12)))
        w = f(Z, yc, best[1])
        nz = np.abs(w) > 1e-8
        sel[k] += nz
        stats[k].append((best[0], nz.sum(), nz[:10].sum(), nz[10:].sum()))

for k in methods:
    s = np.array(stats[k])
    print(f"{k:6s} test MSE {s[:,0].mean():7.3f}  nonzero {s[:,1].mean():6.1f}"
          f"  true found {s[:,2].mean():5.2f}/10  false {s[:,3].mean():6.1f}")
    freq = sel[k]/50
    print(f"        selection frequency: true features {freq[:10].mean():.2f}, "
          f"noise features {freq[10:].mean():.3f}, "
          f"true features always chosen: {(freq[:10]==1).sum()}/10")

2. Ridge keeps all 200 coefficients (no zeros ever) and gets respectable test MSE. LASSO selects roughly 20–40 features, recovering perhaps 5–7 of the 10 true ones — note it typically picks only one or two from each correlated group of five, because once one is in, its correlated twins add little and are thresholded away. Elastic net selects more (40–70), recovers 8–10 of the true ones, and usually has the best or near-best test MSE.

3 & 4. Stability is the striking result. Across 50 resamples, LASSO's selection frequency for the true features is erratic — a given member of a correlated group might be chosen in 30% of runs, its twin in 25%, with essentially arbitrary allocation. Elastic net selects true features far more consistently (frequencies often above 0.8), because the \(\ell_2\) term rewards spreading weight across correlated features rather than concentrating it on one. Ridge is trivially "stable" only because it selects everything.

What this implies about reporting important features. A LASSO run gives you a sparse model that predicts well; it does not give you the set of relevant variables. If a colleague reruns your analysis on next quarter's data and gets a different feature set, that is expected behaviour, not a bug. Three legitimate responses:

• Stability selection: run LASSO on many bootstrap resamples and report features selected in, say, over 60% of them. This has genuine false-discovery guarantees and is what you should do if the feature list is the deliverable.
• Elastic net if correlated features should be reported as a group.
• Say what you mean. "These 12 features suffice for prediction" is a defensible claim. "These 12 features are the causally relevant ones" is not — and no regularisation path can make it so.