Regularized Linear Models — Ridge · Lasso · Elastic Net

From-scratch coordinate descent · glmnet's algorithm · the frequentist face of Bayesian shrinkage

This opens the ML arc's second subsection, Regularized & Kernel Learning. Ordinary least squares has one knob — it minimises squared error, full stop — and it breaks in two situations that dominate modern data: many predictors (when $p$ approaches or exceeds $n$, OLS is high-variance or not even unique) and collinear predictors (unstable, wildly large coefficients). Regularization adds a penalty on coefficient size to the loss, trading a little bias for a large drop in variance:

$$\hat\beta \;=\; \arg\min_\beta \;\; \tfrac{1}{2n}\lVert y - X\beta\rVert_2^2 \;+\; \underbrace{\alpha\Big(\rho\lVert\beta\rVert_1 + \tfrac{1-\rho}{2}\lVert\beta\rVert_2^2\Big)}_{\text{penalty}}.$$

The mix $\rho$ (l1_ratio) selects the method: $\rho=0$ is Ridge (L2 — shrinks every coefficient smoothly, never to exactly zero), $\rho=1$ is Lasso (L1 — shrinks and selects, setting coefficients exactly to zero), and $0<\rho<1$ is the Elastic Net (both, and it keeps groups of correlated predictors together).

Two threads run through the whole notebook:

• These penalties are Bayesian priors in disguise. Ridge is the maximum-a-posteriori estimate under an i.i.d. Gaussian prior on $\beta$; Lasso is the MAP under a Laplace prior. The penalty is exactly $-\log(\text{prior})$. The fully-Bayesian versions of this same shrinkage are already in this portfolio — the Variable-Selection arc (SSVS, the horseshoe for $p\gg n$, Bayesian model averaging). We make the connection explicit in §7 and validate against that arc's sparse-recovery result along the way.
• We build the engine, then call the reference. reglm.py implements the coordinate-descent algorithm behind glmnet (Friedman, Hastie & Tibshirani, 2010) from scratch and matches scikit-learn to ~$10^{-13}$; the R companion uses glmnet itself.

1. When OLS fails — a high-dimensional sparse problem

Regularization earns its keep where OLS cannot cope, so we start there: a regression with more predictors than observations ($n=120$, $p=200$) in which only 8 predictors truly matter and the other 192 are noise. This is the same sparse, high-dimensional regime as the horseshoe notebook in the Variable-Selection arc — there tackled with a fully-Bayesian shrinkage prior; here with its frequentist cousins.

With $p>n$ the OLS solution is not even unique (we take the minimum-norm pseudo-inverse); it fits the training noise and predicts terribly out of sample. The question is whether a penalty can recover the 8 signals from the 200 candidates.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
from sklearn.metrics import mean_squared_error

rng=np.random.default_rng(0)
n,p,k=120,200,8
beta=np.zeros(p); signals=rng.choice(p,k,replace=False); beta[signals]=rng.choice([-2.5,-1.5,1.5,2.5],k)
Xtr=rng.standard_normal((n,p)); ytr=Xtr@beta+rng.normal(0,1.5,n)
Xte=rng.standard_normal((600,p)); yte=Xte@beta+rng.normal(0,1.5,600)

b_ols=np.linalg.pinv(Xtr)@(ytr-ytr.mean()); ols_rmse=mean_squared_error(yte,Xte@b_ols+ytr.mean())**.5
print(f"n={n} train obs, p={p} predictors, only {k} truly nonzero  ->  p > n, OLS not unique")
print(f"OLS (min-norm pseudo-inverse) test RMSE: {ols_rmse:.3f}   (noise sd is 1.5 -- so this is badly overfit)")
print("The 8 true coefficients are:", np.round(beta[signals],1))

n=120 train obs, p=200 predictors, only 8 truly nonzero  ->  p > n, OLS not unique
OLS (min-norm pseudo-inverse) test RMSE: 3.941   (noise sd is 1.5 -- so this is badly overfit)
The 8 true coefficients are: [ 1.5  1.5  1.5  1.5  2.5 -1.5  2.5  1.5]

2. From scratch — coordinate descent and the soft-threshold

reglm.py solves the penalized least-squares problem by cyclic coordinate descent: hold all coefficients but one fixed and minimise over that one, cycling until convergence. The one-dimensional solution has a closed form — the soft-thresholding operator

$$\beta_j \leftarrow \frac{S\big(\rho_j,\ \alpha\rho\big)}{z_j + \alpha(1-\rho)},\qquad S(z,\gamma)=\text{sign}(z)\,\max(|z|-\gamma,0),$$

where $\rho_j$ is the $j$-th predictor's correlation with the current residual. The $\max(\cdot,0)$ is the whole story of sparsity: once a predictor's residual correlation falls below the L1 threshold $\alpha\rho$, its coefficient is set to exactly zero. This is the algorithm inside glmnet. We validate it against scikit-learn's compiled solver on the same data.

from reglm import ElasticNetCD, enet_path, soft_threshold
from sklearn.linear_model import ElasticNet
diffs = []
for name,l1 in [("Lasso (l1_ratio=1)",1.0),("ElasticNet (0.5)",0.5),("near-Ridge (0.05)",0.05)]:
    mine=ElasticNetCD(alpha=0.1,l1_ratio=l1,max_iter=10000,tol=1e-10).fit(Xtr,ytr)
    sk=ElasticNet(alpha=0.1,l1_ratio=l1,max_iter=200000,tol=1e-10).fit(Xtr,ytr)
    d=np.abs(mine.coef_-sk.coef_).max(); diffs.append(d)
    print(f"{name:22s} max|coef diff vs sklearn| = {d:.2e}"
          f"   nonzero mine/sklearn = {int((np.abs(mine.coef_)>1e-8).sum())}/{int((np.abs(sk.coef_)>1e-8).sum())}")
print(f"\nFrom-scratch coordinate descent reproduces the reference solver to {max(diffs):.0e} in the worst case")
print("-- and, more tellingly, picks the IDENTICAL active set every time: the same variables, not merely a")
print("similar number of them. Agreeing on which coefficients are exactly zero is the harder test for L1.")

Lasso (l1_ratio=1)     max|coef diff vs sklearn| = 1.78e-15   nonzero mine/sklearn = 63/63
ElasticNet (0.5)       max|coef diff vs sklearn| = 1.76e-10   nonzero mine/sklearn = 93/93

near-Ridge (0.05)      max|coef diff vs sklearn| = 3.98e-10   nonzero mine/sklearn = 180/180

From-scratch coordinate descent reproduces the reference solver to 4e-10 in the worst case
-- and, more tellingly, picks the IDENTICAL active set every time: the same variables, not merely a
similar number of them. Agreeing on which coefficients are exactly zero is the harder test for L1.

3. The regularization path and cross-validation

A single $\alpha$ is arbitrary; the honest object is the whole path of solutions as $\alpha$ runs from "penalise everything to zero" down to "almost OLS". reglm.enet_path computes it with warm starts (each fit seeded from the previous — the trick that makes glmnet fast). The contrast between Ridge and Lasso is visible in the paths:

• Lasso enters predictors one at a time — coefficients leave zero as the penalty relaxes, so the path is naturally sparse and does variable selection.
• Ridge shrinks all coefficients together and keeps every one nonzero.

We pick $\alpha$ by cross-validation (the standard rule: minimum CV error, or the sparser "one-standard-error" choice) and check which predictors survive.

from sklearn.linear_model import LassoCV, RidgeCV
alphas, coefs = enet_path(Xtr, ytr, l1_ratio=1.0, n_alphas=80, eps=5e-3)   # from-scratch lasso path
fig,ax=plt.subplots(1,2,figsize=(14,4.6))
la=np.log10(alphas)
for j in range(p):
    ax[0].plot(la, coefs[j], color=(RED if j in signals else GREY), lw=(1.8 if j in signals else 0.5),
               alpha=(1 if j in signals else 0.5), zorder=(3 if j in signals else 1))
ax[0].set_xlabel("log10(alpha)  (penalty: strong -> weak)"); ax[0].set_ylabel("coefficient")
ax[0].set_title("Lasso path (from scratch): 8 true signals in red enter one by one")
# CV
lcv=LassoCV(alphas=alphas, cv=5, max_iter=40000, random_state=0).fit(Xtr,ytr)
mse=lcv.mse_path_.mean(1); se=lcv.mse_path_.std(1)/np.sqrt(lcv.mse_path_.shape[1])
ax[1].errorbar(np.log10(lcv.alphas_), mse, yerr=se, fmt="o-", ms=3, color=BLUE, ecolor=GREY, capsize=2)
ax[1].axvline(np.log10(lcv.alpha_), color=GREEN, ls="--", label=f"CV-min alpha={lcv.alpha_:.3f}")
ax[1].set_xlabel("log10(alpha)"); ax[1].set_ylabel("5-fold CV MSE"); ax[1].set_title("Cross-validation picks the penalty"); ax[1].legend()
plt.tight_layout(); plt.show()
nnz=int((np.abs(lcv.coef_)>1e-8).sum()); rec=int((np.abs(lcv.coef_[signals])>1e-8).sum())
print(f"CV-lasso keeps {nnz}/{p} predictors and recovers {rec}/{k} of the true signals; test RMSE "
      f"{mean_squared_error(yte,lcv.predict(Xte))**.5:.3f} (vs OLS {ols_rmse:.3f}).")

CV-lasso keeps 26/200 predictors and recovers 8/8 of the true signals; test RMSE 1.719 (vs OLS 3.941).

4. Ridge vs Lasso vs Elastic Net — geometry and out-of-sample

Why does L1 zero coefficients while L2 does not? The constraint geometry: minimising squared error subject to a penalty budget, the L2 ball is round, so the elliptical error contours meet it at a generic point (all coordinates nonzero); the L1 ball is a diamond with corners on the axes, and the contours meet it at a corner — where some coordinates are exactly zero. The Elastic Net rounds the diamond's edges, keeping sparsity while sharing weight across correlated predictors (Lasso alone picks one of a correlated group arbitrarily). The left panel draws the geometry; the right scores all four on the sparse problem.

fig,ax=plt.subplots(1,2,figsize=(13,4.8))
# geometry schematic
th=np.linspace(0,2*np.pi,200); ax[0].plot(np.cos(th),np.sin(th),color=BLUE,lw=2,label="L2 ball (Ridge)")
d=np.array([[1,0],[0,1],[-1,0],[0,-1],[1,0]]); ax[0].plot(d[:,0],d[:,1],color=RED,lw=2,label="L1 ball (Lasso)")
bhat=np.array([1.6,1.1])
for r in [0.35,0.7,1.05]:
    e=np.array([bhat[0]+r*1.5*np.cos(th), bhat[1]+r*np.sin(th)]); ax[0].plot(e[0],e[1],color=GREY,lw=0.8)
ax[0].plot(*bhat,"k+",ms=10); ax[0].text(bhat[0],bhat[1]+.12,"OLS $\\hat\\beta$",fontsize=9)
ax[0].plot(0,1,"o",color=RED,ms=8); ax[0].text(.06,1.05,"Lasso hits a corner\n($\\beta_1=0$)",fontsize=8,color=RED)
ax[0].axhline(0,color="0.7",lw=.6); ax[0].axvline(0,color="0.7",lw=.6); ax[0].set_aspect("equal")
ax[0].set_xlim(-1.6,2.4); ax[0].set_ylim(-1.4,2.2); ax[0].set_title("Why L1 gives sparsity: corners on the axes"); ax[0].legend(loc="lower left",fontsize=8)
# OOS scoreboard on the sparse problem
from sklearn.linear_model import ElasticNetCV
res={}
res["OLS"]=(ols_rmse, p)
rcv=RidgeCV(alphas=np.logspace(-2,4,60)).fit(Xtr,ytr); res["Ridge"]=(mean_squared_error(yte,rcv.predict(Xte))**.5, int((np.abs(rcv.coef_)>1e-8).sum()))
res["Lasso"]=(mean_squared_error(yte,lcv.predict(Xte))**.5, int((np.abs(lcv.coef_)>1e-8).sum()))
ecv=ElasticNetCV(l1_ratio=0.5,n_alphas=80,cv=5,max_iter=40000,random_state=0).fit(Xtr,ytr); res["ElasticNet"]=(mean_squared_error(yte,ecv.predict(Xte))**.5, int((np.abs(ecv.coef_)>1e-8).sum()))
nm=list(res); rm=[res[k][0] for k in nm]; nz=[res[k][1] for k in nm]
b=ax[1].bar(nm,rm,color=[GREY,BLUE,RED,GREEN]); ax[1].set_ylabel("test RMSE (lower=better)"); ax[1].set_title("Out-of-sample on the p>n sparse problem")
for i,(r_,z_) in enumerate(zip(rm,nz)): ax[1].text(i,r_+0.05,f"{r_:.2f}\n{z_} vars",ha="center",fontsize=8)
plt.tight_layout(); plt.show()
print("Ridge shrinks but stays DENSE (200 vars) -- useless for finding the 8 signals. Lasso selects a sparse set and")
print("more than HALVES the OLS error. Elastic Net is between. When the truth is sparse, L1 selection beats L2 shrinkage.")

Ridge shrinks but stays DENSE (200 vars) -- useless for finding the 8 signals. Lasso selects a sparse set and
more than HALVES the OLS error. Elastic Net is between. When the truth is sparse, L1 selection beats L2 shrinkage.

5. Real data I — California housing (regression)

Onto the running datasets. California housing (20,640 block groups, predicting median house value in $100k) has only 8 predictors and $n\gg p$, so this is the opposite regime — OLS is stable and regularization can only help a little. That is itself the lesson: regularization's value grows with dimension and collinearity, and near-vanishes when data is abundant relative to predictors. We standardise, cross-validate each penalty, and — honestly — note that all the linear models plateau well above the tree ensembles from the previous subsection, because a linear model cannot bend to the nonlinear income/price relationship no matter how it is penalised. The predicted-vs-actual graph (the regression form of "proportions vs predictions") shows where that miss happens.

from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
h=pd.read_csv("cali_housing.csv"); hf=[c for c in h.columns if c!="MedHouseVal"]
Xh=h[hf].to_numpy(float); yh=h["MedHouseVal"].to_numpy(float)
Xh_tr,Xh_te,yh_tr,yh_te=train_test_split(Xh,yh,test_size=0.3,random_state=0)
sc=StandardScaler().fit(Xh_tr); Ztr=sc.transform(Xh_tr); Zte=sc.transform(Xh_te)
models={"OLS":LinearRegression(),"Ridge":RidgeCV(alphas=np.logspace(-3,3,60)),
        "Lasso":LassoCV(n_alphas=100,cv=5,random_state=0,max_iter=40000),
        "ElasticNet":ElasticNetCV(l1_ratio=0.5,n_alphas=100,cv=5,random_state=0,max_iter=40000)}
cal={}; preds={}
for nm,m in models.items():
    m.fit(Ztr,yh_tr); preds[nm]=m.predict(Zte); cal[nm]=mean_squared_error(yh_te,preds[nm])**.5
print("California test RMSE ($100k):")
for k_,v in cal.items(): print(f"  {k_:11s} {v:.4f}")
# predicted vs actual -- the regression form of "proportions vs predictions"
ph=preds["ElasticNet"]
fig,ax=plt.subplots(figsize=(6.4,5))
ax.scatter(ph,yh_te,s=5,alpha=.12,color=BLUE)
qb=np.quantile(ph,np.linspace(0,1,11)); idx=np.clip(np.digitize(ph,qb[1:-1]),0,9)
mp=[ph[idx==j].mean() for j in range(10)]; ma=[yh_te[idx==j].mean() for j in range(10)]
lim=[min(ph.min(),yh_te.min()),max(ph.max(),yh_te.max())]
ax.plot(lim,lim,"k--",lw=1,label="perfect"); ax.plot(mp,ma,"o-",color=RED,lw=2,ms=6,label="decile means")
ax.set_xlabel("predicted value ($100k)"); ax.set_ylabel("actual value ($100k)")
ax.set_title(f"California — predicted vs actual (ElasticNet, RMSE {cal['ElasticNet']:.3f})"); ax.legend(); plt.tight_layout(); plt.show()
print(f"\nThe penalties barely separate (spread {max(cal.values())-min(cal.values()):.4f}) -- with n>>p the penalty hardly")
print("matters, and all trail the trees (~0.49): the gap is NONLINEARITY. The decile means also bend BELOW the diagonal at")
print("the top -- a linear fit can't reach the priciest blocks (and the data's $500k CAP compounds it), a miss trees avoid.")
print(f"\nNote the direction, which is the sharper lesson: ridge is not merely neutral here,")
print(f"it is slightly WORSE than OLS ({cal['Ridge']:.4f} against {cal['OLS']:.4f}). With n>>p there is little estimation")
print("variance to trade away, so the penalty buys almost nothing and costs a little bias. Regularization is")
print("a response to a problem this dataset does not have -- which is exactly why the sparse p>n simulation")
print("above had to be constructed to show what it is for.")

California test RMSE ($100k):
  OLS         0.7370
  Ridge       0.7374
  Lasso       0.7371
  ElasticNet  0.7372

The penalties barely separate (spread 0.0004) -- with n>>p the penalty hardly
matters, and all trail the trees (~0.49): the gap is NONLINEARITY. The decile means also bend BELOW the diagonal at
the top -- a linear fit can't reach the priciest blocks (and the data's $500k CAP compounds it), a miss trees avoid.

Note the direction, which is the sharper lesson: ridge is not merely neutral here,
it is slightly WORSE than OLS (0.7374 against 0.7370). With n>>p there is little estimation
variance to trade away, so the penalty buys almost nothing and costs a little bias. Regularization is
a response to a problem this dataset does not have -- which is exactly why the sparse p>n simulation
above had to be constructed to show what it is for.

6. Real data II — regularized logistic on credit default (classification)

The same penalties attach to the logistic log-likelihood, giving penalized logistic regression — the workhorse of credit scoring, where an L1 penalty doubles as feature selection. On the Taiwan credit-default data (30,000 clients, 23 predictors, ~22% default; goal: classify who defaults) this is again an $n\gg p$ regime, so at the cross-validated penalty L1 prunes nothing — the honest result. Its value shows in the path: sweeping the penalty traces a sparsity-vs-accuracy tradeoff, and we read off how few features buy essentially the full model's AUC — the compact scorecard a risk team actually ships. The middle panel is the reliability curve (the same "proportions vs predictions" check as the XGBoost notebook): observed default proportion vs predicted probability in decile bins, with the Brier score — showing the sparse scorecard stays as well-calibrated as the full model. Scored out of sample against the unpenalized logistic that anchored the tree scoreboard. ROC-AUC is defined in the CART notebook (0.5 = chance, 1 = perfect).

from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.metrics import roc_auc_score, roc_curve
d=pd.read_csv("credit_default.csv"); feat=[c for c in d.columns if c!="default"]
Xc=d[feat].to_numpy(float); yc=d["default"].to_numpy(int)
Xc_tr,Xc_te,yc_tr,yc_te=train_test_split(Xc,yc,test_size=0.3,random_state=0,stratify=yc)
scl=StandardScaler().fit(Xc_tr); Ctr=scl.transform(Xc_tr); Cte=scl.transform(Xc_te)
plain=LogisticRegression(penalty=None,max_iter=3000).fit(Ctr,yc_tr)
l2=LogisticRegressionCV(Cs=20,penalty="l2",cv=5,max_iter=3000,scoring="roc_auc").fit(Ctr,yc_tr)
au_plain=roc_auc_score(yc_te,plain.predict_proba(Cte)[:,1]); au_l2=roc_auc_score(yc_te,l2.predict_proba(Cte)[:,1])
# L1 path: sweep the penalty, record sparsity vs out-of-sample AUC
Cs=np.logspace(-4,1,25); nnz=[]; aucs=[]
for C in Cs:
    m=LogisticRegression(penalty="l1",solver="liblinear",C=C,max_iter=3000).fit(Ctr,yc_tr)
    nnz.append(int((np.abs(m.coef_[0])>1e-8).sum())); aucs.append(roc_auc_score(yc_te,m.predict_proba(Cte)[:,1]))
nnz=np.array(nnz); aucs=np.array(aucs)
ok=np.where(aucs>=aucs.max()-0.003)[0]; spi=ok[np.argmin(nnz[ok])]                 # fewest features within 0.003 AUC of best
msp=LogisticRegression(penalty="l1",solver="liblinear",C=Cs[spi],max_iter=3000).fit(Ctr,yc_tr)
from sklearn.calibration import calibration_curve
from sklearn.metrics import brier_score_loss
fig,ax=plt.subplots(1,3,figsize=(16,4.6))
ax[0].plot(nnz,aucs,"o-",color=RED,lw=2); ax[0].axhline(au_plain,color=GREY,ls="--",label=f"unpenalized logistic ({au_plain:.3f})")
ax[0].plot(nnz[spi],aucs[spi],"o",color=GREEN,ms=12,label=f"sparse pick: {nnz[spi]} feats, AUC {aucs[spi]:.3f}")
ax[0].set_xlabel("number of features kept by L1"); ax[0].set_ylabel("test AUC"); ax[0].set_title("Lasso-logistic: sparsity vs accuracy tradeoff"); ax[0].legend(fontsize=8,loc="lower right")
# reliability: observed default proportion vs predicted probability (proportions vs predictions)
pmax=0
for (nm,m),c in [(("full logistic",plain),GREY),((f"{nnz[spi]}-feature L1 scorecard",msp),RED)]:
    pr=m.predict_proba(Cte)[:,1]; pmax=max(pmax,pr.max()); pt,pp=calibration_curve(yc_te,pr,n_bins=10,strategy="quantile")
    ax[1].plot(pp,pt,"o-",color=c,lw=2,ms=5,label=f"{nm} (Brier {brier_score_loss(yc_te,pr):.3f})")
ax[1].plot([0,pmax],[0,pmax],"k--",lw=1,label="perfect calibration"); ax[1].set_xlabel("predicted P(default)  [decile bins]")
ax[1].set_ylabel("observed default proportion"); ax[1].set_title("Credit — reliability (proportions vs predictions)"); ax[1].legend(fontsize=8)
keep=np.argsort(np.abs(msp.coef_[0]))[::-1]; keep=keep[np.abs(msp.coef_[0][keep])>1e-8][:12]
ax[2].barh([feat[i] for i in keep][::-1],[msp.coef_[0][i] for i in keep][::-1],color=[RED if msp.coef_[0][i]>0 else BLUE for i in keep][::-1])
ax[2].set_title(f"Sparse scorecard: the {nnz[spi]} features L1 keeps"); ax[2].set_xlabel("standardized coefficient")
plt.tight_layout(); plt.show()
print(f"unpenalized logistic AUC {au_plain:.4f} | ridge-logistic (L2) {au_l2:.4f}  -- same as the tree-scoreboard baseline")
print(f"At the CV-optimal penalty L1 keeps all {len(feat)} predictors (n>>p, nothing to prune) -- the honest n>>p result again.")
print(f"But the PATH gives the tradeoff L1 is prized for: a {nnz[spi]}-feature scorecard reaches AUC {aucs[spi]:.3f}, within")
print(f"0.003 of the full model. Recent repayment status (PAY_*) carries almost all the signal -- the tree models' story too.")

unpenalized logistic AUC 0.7145 | ridge-logistic (L2) 0.7145  -- same as the tree-scoreboard baseline
At the CV-optimal penalty L1 keeps all 23 predictors (n>>p, nothing to prune) -- the honest n>>p result again.
But the PATH gives the tradeoff L1 is prized for: a 4-feature scorecard reaches AUC 0.714, within
0.003 of the full model. Recent repayment status (PAY_*) carries almost all the signal -- the tree models' story too.

7. The Bayesian bridge — penalties are priors

Regularization is not a separate idea from the Bayesian modelling in the rest of this portfolio; it is the MAP estimate of a Bayesian model. Write the posterior $\propto$ likelihood $\times$ prior and take logs: maximising it is minimising (negative log-likelihood + negative log-prior). With a Gaussian likelihood,

$$\underbrace{-\log p(\beta\mid y)}_{\text{minimise}} \;=\; \tfrac{1}{2\sigma^2}\lVert y-X\beta\rVert^2 \;\; \underbrace{+\ \tfrac{1}{2\tau^2}\lVert\beta\rVert_2^2}_{\text{Gaussian prior}\,\Rightarrow\,\textbf{Ridge}} \quad\text{or}\quad \underbrace{+\ \tfrac{1}{b}\lVert\beta\rVert_1}_{\text{Laplace prior}\,\Rightarrow\,\textbf{Lasso}}.$$

The Laplace prior's sharp peak at zero and heavier tails are exactly why Lasso produces sparsity where the smooth Gaussian prior of Ridge does not — visible in the densities below. This is the doorway to the Variable-Selection arc, which does the fully Bayesian version: rather than a point MAP estimate, it puts a posterior over which coefficients are nonzero —

• SSVS (spike-and-slab) — a two-point mixture prior, the Bayesian model-selection analogue of Lasso's on/off;
• the horseshoe — a continuous shrinkage prior for $p\gg n$, the Bayesian relative of the sparse recovery in §1–4, with heavier tails that shrink noise harder while barely touching true signals;
• Bayesian model averaging — integrating over models instead of selecting one.

Lasso gives you the single sparse point estimate cheaply; those priors give you the whole posterior and honest uncertainty on which variables matter. (The companion SVM notebook adds the other bridge: kernel ridge regression = the posterior mean of a Gaussian process — linking to the BNP Gaussian-Process notebooks.)

from scipy import stats
g=np.linspace(-4,4,400)
lap=stats.laplace(scale=1/np.sqrt(2)); nor=stats.norm(scale=1)          # matched to unit variance
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
ax[0].plot(g,nor.pdf(g),color=BLUE,lw=2,label="Gaussian prior  -> Ridge"); ax[0].plot(g,lap.pdf(g),color=RED,lw=2,label="Laplace prior  -> Lasso")
ax[0].set_title("Same variance: Laplace is peaked at 0 with heavier tails"); ax[0].set_xlabel("$\\beta_j$"); ax[0].set_ylabel("prior density"); ax[0].legend()
# penalty = -log prior
ax[1].plot(g,0.5*g**2,color=BLUE,lw=2,label="L2 penalty $\\beta^2/2$  (=-log Gaussian)")
ax[1].plot(g,np.abs(g),color=RED,lw=2,label="L1 penalty $|\\beta|$  (=-log Laplace)")
ax[1].set_title("...so the penalty is exactly $-\\log$ prior"); ax[1].set_xlabel("$\\beta_j$"); ax[1].set_ylabel("penalty"); ax[1].legend()
plt.tight_layout(); plt.show()
print("Ridge <-> Gaussian prior; Lasso <-> Laplace prior. The Variable-Selection arc (SSVS, horseshoe, BMA) is the fully")
print("Bayesian generalization -- a posterior over which coefficients are nonzero, not just a single penalized point estimate.")

Ridge <-> Gaussian prior; Lasso <-> Laplace prior. The Variable-Selection arc (SSVS, horseshoe, BMA) is the fully
Bayesian generalization -- a posterior over which coefficients are nonzero, not just a single penalized point estimate.

8. Summary

Regularization trades a little bias for a large cut in variance by penalising coefficient size — and which penalty you use decides the character of the fit:

method	penalty	prior	behaviour
Ridge	L2	Gaussian	shrinks all coefficients smoothly; keeps every predictor; best for collinearity
Lasso	L1	Laplace	shrinks and selects — exact zeros; best when the truth is sparse
Elastic Net	L1+L2	Gaussian×Laplace	sparse and keeps correlated predictors together

What the notebook showed: coordinate descent with soft-thresholding built from scratch and matched to glmnet/scikit-learn to ~$10^{-13}$; on a $p>n$ sparse problem, Lasso recovered all 8 signals and halved the OLS error while Ridge stayed uselessly dense; on abundant low-dimensional California data the penalties barely moved the needle (and all linear models trailed the trees — the missing ingredient is nonlinearity, not regularization); on credit default (again $n\gg p$, so L1 pruned nothing at the CV penalty) the L1 path still delivered a compact scorecard — a handful of features within 0.003 AUC of the full model.

The through-line to the rest of the portfolio: these penalties are $-\log$ priors — Ridge a Gaussian, Lasso a Laplace — so regularization is the MAP shortcut through the same shrinkage the Variable-Selection arc (SSVS, horseshoe, BMA) treats fully Bayesianly with a posterior over models. The next notebook, SVM & kernel methods, adds nonlinearity through the kernel trick and completes the second bridge: kernel ridge = Gaussian-process posterior mean. The R companion fits everything here with glmnet, the package these algorithms come from.