Regularization: Ridge, Lasso, and the Geometry of Sparsity

Ridge, lasso, and elastic net penalize coefficient size while fitting the data. They are useful when several coefficient vectors fit similarly well: the penalty helps choose among them. Ridge uses squared coefficients, lasso uses absolute values, and elastic net combines the two. We will examine how these choices affect small coefficients and correlated features, then tune the penalties on housing data.

Shrinkage and exact zeros

Let \(X\) be the feature matrix, \(y\) the target vector, and \(\beta\) the coefficients. Start with centered data and orthonormal columns, meaning \(X^\top X=I\): each column has length one and different columns have zero dot product. No intercept is needed in this centered calculation.

Define ridge by minimizing \(\tfrac12\|y-X\beta\|_2^2+\tfrac\lambda2\|\beta\|_2^2\), and lasso by minimizing \(\tfrac12\|y-X\beta\|_2^2+\lambda\|\beta\|_1\), where \(\lambda\geq0\) controls the penalty. Here \(\|\beta\|_2^2=\sum_j\beta_j^2\) and \(\|\beta\|_1=\sum_j|\beta_j|\). With these conventions, ridge divides each least-squares coefficient by \(1+\lambda\). Lasso applies the soft threshold \(\operatorname{sign}(b)\max(|b|-\lambda,0)\) to each least-squares coefficient \(b\).

lam = 1.0
for b_ols in (-2.0, -0.5, 0.0, 0.3, 1.0, 3.0):
    ridge = b_ols / (1 + lam)
    lasso = (1 if b_ols > 0 else -1) * max(abs(b_ols) - lam, 0)
    print(f"OLS {b_ols:+5.2f}  ridge {ridge:+7.4f}  lasso {lasso:+7.4f}")
# OLS -2.00  ridge -1.0000  lasso -1.0000
# OLS -0.50  ridge -0.2500  lasso +0.0000
# OLS +0.00  ridge +0.0000  lasso +0.0000
# OLS +0.30  ridge +0.1500  lasso +0.0000
# OLS +1.00  ridge +0.5000  lasso +0.0000
# OLS +3.00  ridge +1.5000  lasso +2.0000

Here \(\lambda=1\), so ridge halves each coefficient. For an OLS coefficient of 0.3, ridge gives 0.15 and lasso gives zero; for 3, they give 1.5 and 2. Lasso sets coefficients with magnitude at most \(\lambda\) to zero. Above that threshold it reduces their magnitude by \(\lambda\). These coordinate formulas rely on orthonormal columns. With correlated columns, individual ridge coefficients can cross zero as the penalty changes, even though ridge does not perform lasso-style thresholding.

For one coordinate, the loss above becomes \(\tfrac12(\beta-b)^2\), up to a constant. Ridge sets its derivative to zero: \(\beta-b+\lambda\beta=0\), giving \(\beta=b/(1+\lambda)\). For a positive lasso solution the equation is \(\beta-b+\lambda=0\); for a negative solution it is \(\beta-b-\lambda=0\). At zero, the absolute-value penalty has a range of supporting slopes, from \(-\lambda\) to \(\lambda\), called its subgradient. Zero is optimal when that range can balance the loss slope, which happens when \(|b|\leq\lambda\).

This explains why an L1 penalty can leave a coefficient exactly zero over a range of data values. Non-differentiability alone is not a general guarantee of sparse coefficients; the particular penalty and loss determine what becomes zero. Elastic net retains this thresholding behavior when its L1 component is positive, while its L2 component also encourages smaller, more evenly distributed coefficients.

What each does to correlated features

The library conventions differ from the first calculation. Ridge minimizes the sum of squared residuals plus alpha times the squared coefficient norm. Lasso divides that residual sum by \(2n\), then adds alpha times the absolute coefficient sum. Here \(n\) is the number of training rows. ElasticNet uses the same averaged loss as Lasso, with L1 weight \(\alpha r\) and squared-L2 weight \(\alpha(1-r)/2\), where \(r\) is l1_ratio. Equal alpha values across these estimators do not mean equal penalty strength. Their intercepts are unpenalized.

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet

rng = np.random.default_rng(0)
n = 200
z = rng.normal(size=n)
X = np.column_stack([z + 0.01 * rng.normal(size=n),     # two near-copies of z
                     z + 0.01 * rng.normal(size=n),
                     rng.normal(size=n)])               # and one noise column
y = 3 * z + 0.5 * rng.normal(size=n)

for name, m in (("OLS", LinearRegression()),
                ("ridge a=1", Ridge(alpha=1.0)),
                ("lasso a=0.1", Lasso(alpha=0.1)),
                ("enet a=0.1 l1=0.5", ElasticNet(alpha=0.1, l1_ratio=0.5))):
    print(f"{name:18s} coef {np.round(m.fit(X, y).coef_, 4)}")
# OLS                coef [ 2.8705  0.1525 -0.0133]
# ridge a=1          coef [ 1.5306  1.4823 -0.0106]
# lasso a=0.1        coef [ 2.9139e+00  4.0000e-04 -0.0000e+00]
# enet a=0.1 l1=0.5  coef [ 1.4462  1.4421 -0.    ]

OLS assigns 2.87 and 0.15 to the near-copy features. Their individual coefficients are difficult to distinguish reliably when the columns contain almost the same information. Ridge assigns about 1.53 and 1.48. If two columns were exactly identical, keeping their coefficient sum fixed would preserve predictions, and the squared penalty would be smallest at an equal split: \(1.5^2+1.5^2=4.5\), compared with \(3^2+0^2=9\). With near-copies, the fit term also matters. Lasso strongly concentrates the weight in one column; the displayed small second coefficient is not an exact zero.

For prediction, compare the penalties on held-out data and examine how sensitive the fitted coefficients are to changes in the training sample. A sensor retained by lasso is a candidate to keep when reducing measurement costs. Before removing the others, check measurement quality, availability, and performance without them. A zero coefficient does not establish that a sensor is scientifically irrelevant, and an even ridge split does not identify separate causal effects.

Elastic net splits this pair nearly evenly and sets the noise coefficient to zero. Its L2 component can help retain correlated predictors together while L1 can set other coefficients to zero. Retaining a group still needs to be evaluated against the prediction or measurement goal.

Effective degrees of freedom

For fixed \(X\) and a fixed ridge penalty, predictions are a linear function of the targets: \(\hat y=Hy\), where \(H=X(X^\top X+\lambda I)^{-1}X^\top\). The trace, the sum of the diagonal entries, measures total sensitivity: \(\operatorname{tr}(H)=\sum_i\partial\hat y_i/\partial y_i\). This is ridge’s effective degrees of freedom. The following example fits no intercept and has full column rank, so the inverse also exists at \(\lambda=0\).

import numpy as np

rng = np.random.default_rng(1)
n, p = 100, 20
X = rng.normal(size=(n, p))

for lam in (0.0, 1.0, 10.0, 100.0, 1000.0):
    H = X @ np.linalg.solve(X.T @ X + lam * np.eye(p), X.T)
    print(f"lambda {lam:7.1f}  effective df {np.trace(H):7.3f}   (p = {p})")
# lambda     0.0  effective df  20.000   (p = 20)
# lambda     1.0  effective df  19.748   (p = 20)
# lambda    10.0  effective df  17.799   (p = 20)
# lambda   100.0  effective df   9.513   (p = 20)
# lambda  1000.0  effective df   1.807   (p = 20)

At \(\lambda=0\), effective degrees of freedom is \(p=20\). At \(\lambda=1000\), it is about 1.8 while all twenty coefficients remain in the model. This quantifies reduced sensitivity to the training targets; it does not make the fit equivalent to a model with two selected features.

If \(\sigma_j\) are the singular values of \(X\), the trace is \(\sum_j\sigma_j^2/(\sigma_j^2+\lambda)\). These fractions are the shrinkage factors along the corresponding fitted-value directions. A direction with \(\sigma_j^2=1\) keeps half its fitted component at \(\lambda=1\); one with \(\sigma_j^2=100\) keeps \(100/101\). Ridge suppresses weakly determined directions more strongly, but positive singular directions are not removed exactly at a finite penalty.

A trace of 9.5 is close to the value 10 for a full-rank, ten-column OLS fit without an intercept. The comparison concerns this sensitivity measure, not equal predictions or equal test error. A separately fitted, unpenalized intercept adds one degree of freedom to the centered-data calculation. Selecting the penalty using the targets adds adaptivity that the fixed-penalty trace does not include.

Both penalties on the same real data

The California housing data include the correlated features AveRooms and AveBedrms. We reserve a test set, then compare candidate penalties using the same five training-data folds and mean squared error. The scaler is inside each pipeline, so every fold estimates its means and standard deviations from that fold’s training rows. After selection, each pipeline is refitted on all training rows and evaluated on the reserved test set. GridSearchCV tries each grid value; model__alpha names the alpha parameter of the pipeline’s model step. Its scoring convention maximizes a score, so neg_mean_squared_error selects the smallest mean squared error.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, KFold, GridSearchCV
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import root_mean_squared_error

X, y = fetch_california_housing(return_X_y=True, as_frame=True)
names = list(X.columns)
Xtr, Xte, ytr, yte = train_test_split(
    X.to_numpy(), y.to_numpy(), test_size=0.3, random_state=0)
folds = KFold(n_splits=5, shuffle=True, random_state=0)

def pipeline(model):
    return Pipeline([("scale", StandardScaler()), ("model", model)])

fits = {"OLS": pipeline(LinearRegression()).fit(Xtr, ytr)}
for name, model, grid in (
    ("ridge", Ridge(), np.logspace(-3, 3, 25)),
    ("lasso", Lasso(max_iter=50000, tol=1e-8), np.logspace(-3, 1, 25))
):
    search = GridSearchCV(pipeline(model), {"model__alpha": grid},
                          cv=folds, scoring="neg_mean_squared_error")
    search.fit(Xtr, ytr)
    fits[name] = search.best_estimator_
    print(f"{name} alpha {search.best_params_['model__alpha']:.4f}")

print(f"  {'feature':12s} {'OLS':>9s} {'ridge':>9s} {'lasso':>9s}")
coefs = [fits[name].named_steps["model"].coef_ for name in fits]
for i, name in enumerate(names):
    print(f"  {name:12s}" + "".join(f" {c[i]:9.4f}" for c in coefs))
for name, fit in fits.items():
    c = fit.named_steps["model"].coef_
    rmse = root_mean_squared_error(yte, fit.predict(Xte))
    print(f"{name:5s} test RMSE {rmse:.4f}  nonzero {np.count_nonzero(c)}/8")
# ridge alpha 31.6228
# lasso alpha 0.0022
#   feature            OLS     ridge     lasso
#   MedInc          0.8449    0.8439    0.8345
#   HouseAge        0.1157    0.1185    0.1177
#   AveRooms       -0.2702   -0.2656   -0.2446
#   AveBedrms       0.2908    0.2844    0.2653
#   Population     -0.0108   -0.0097   -0.0077
#   AveOccup       -0.0281   -0.0284   -0.0266
#   Latitude       -0.8753   -0.8487   -0.8522
#   Longitude      -0.8496   -0.8228   -0.8253
# OLS   test RMSE 0.7370  nonzero 8/8
# ridge test RMSE 0.7372  nonzero 8/8
# lasso test RMSE 0.7370  nonzero 8/8

Lasso selects alpha 0.0022 and retains all eight features here. That describes the fitted model; it does not prove that every feature contributes useful signal. Cross-validation chooses from a finite grid to reduce prediction error; it does not test which coefficients are truly nonzero. The coefficient count describes the resulting fit, and the test errors provide a separate comparison of predictive performance.

The coefficients in the table are for standardized training features: each describes a change in the prediction per training standard deviation of that feature, holding the others fixed. The target is measured in hundreds of thousands of dollars. Shrinking the opposite-signed room and bedroom coefficients toward zero is different from the equal positive split in the synthetic example. Correlation alone does not determine the signs or the shrinkage of individual coefficients.

The selected penalties and errors are specific to this split, these grids, and squared-error prediction. The test RMSE values are nearly equal: 0.7370 for OLS, 0.7372 for ridge, and 0.7370 for lasso, roughly $73,700 in the target’s units. These results do not establish that one method will consistently win. To choose a method, use training-data validation; repeatedly changing the analysis to improve this test score would turn the test set into another validation set.

Exercises

1. Walk the lasso path. Fit a lasso at increasing penalties on data with three real coefficients and five zeros, and report the number of nonzero coefficients and their values at each.

Compare which coefficients become zero in this run and how the surviving estimates differ from the known generating values.

import numpy as np
from sklearn.linear_model import Lasso

rng = np.random.default_rng(0)
n, p = 200, 8
X = rng.normal(size=(n, p))
y = X @ np.array([3.0, -2.0, 1.0, 0, 0, 0, 0, 0]) + 0.5 * rng.normal(size=n)

for a in (0.001, 0.05, 0.2, 0.5, 1.0, 2.0):
    c = Lasso(alpha=a).fit(X, y).coef_
    print(f"alpha {a:6.3f}  nonzero {np.sum(np.abs(c) > 1e-8)}  coef {np.round(c, 3)}")
# alpha  0.001  nonzero 8  coef [ 3.016 -1.963  0.955 -0.048  0.028 -0.018  0.021 -0.007]
# alpha  0.050  nonzero 4  coef [ 2.957 -1.91   0.893 -0.006  0.    -0.     0.     0.   ]
# alpha  0.200  nonzero 3  coef [ 2.776 -1.743  0.689 -0.     0.     0.     0.     0.   ]
# alpha  0.500  nonzero 3  coef [ 2.414 -1.408  0.28  -0.     0.     0.     0.     0.   ]
# alpha  1.000  nonzero 2  coef [ 1.844 -0.877  0.    -0.     0.     0.     0.     0.   ]
# alpha  2.000  nonzero 1  coef [ 0.743 -0.     0.     0.     0.     0.     0.     0.   ]

2. Changing units changes the penalty. Regularize two features with equal true effects, multiplying the numeric values of one feature by 1000, and compare the cross-validated score with and without standardization.

You should get: a scaled score that is identical across all three unit choices and a raw score that is not.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(1)
n = 500
a = rng.normal(size=n)
b = rng.normal(size=n)
y = 1.5 * a + 1.5 * b + 0.3 * rng.normal(size=n)        # both matter equally
alpha = 1e4

scaled = lambda: Pipeline([("s", StandardScaler()), ("m", Ridge(alpha=alpha))])
for label, X in (("a in units, b in units", np.column_stack([a, b])),
                 ("a in units, b x1000",    np.column_stack([a, b * 1000])),
                 ("a x1000, b in units",    np.column_stack([a * 1000, b]))):
    raw = cross_val_score(Ridge(alpha=alpha), X, y, cv=5).mean()
    std = cross_val_score(scaled(), X, y, cv=5).mean()
    print(f"{label:>24}  raw {raw:7.4f}   scaled {std:7.4f}"
          f"   raw coefs {np.round(Ridge(alpha=alpha).fit(X, y).coef_, 6)}")
#   a in units, b in units  raw  0.0579   scaled  0.0585   raw coefs [0.060436 0.080306]
#      a in units, b x1000  raw  0.5802   scaled  0.0585   raw coefs [0.060466 0.001526]
#      a x1000, b in units  raw  0.4478   scaled  0.0585   raw coefs [0.001511 0.080336]

3. The grouping effect. Give a lasso and an elastic net four near-identical copies of one predictor plus three noise columns, and compare how each allocates the coefficients.

Compare the allocation within the four-column group, the coefficient sum, and whether each noise coefficient is exactly zero.

import numpy as np
from sklearn.linear_model import Lasso, ElasticNet

rng = np.random.default_rng(2)
n = 300
z = rng.normal(size=n)
group = np.column_stack([z + 0.05 * rng.normal(size=n) for _ in range(4)])
X = np.column_stack([group, rng.normal(size=(n, 3))])
y = 2 * z + 0.4 * rng.normal(size=n)

for name, m in (("lasso", Lasso(alpha=0.15)),
                ("enet l1=0.2", ElasticNet(alpha=0.15, l1_ratio=0.2))):
    c = m.fit(X, y).coef_
    print(f"{name:12s} group {np.round(c[:4], 4)}  noise {np.round(c[4:], 4)}"
          f"  group sum {c[:4].sum():.4f}")
# lasso        group [0.3469 0.4223 1.066  0.    ]  noise [0. 0. 0.]  group sum 1.8352
# enet l1=0.2  group [0.4766 0.477  0.4895 0.4564]  noise [0.0059 0.     0.    ]  group sum 1.8995