03 Regularization — Exercises¶
Test your understanding of Ridge, Lasso, and regularization geometry.
Prerequisites. Read theory.md and work through first_principles.ipynb before attempting these.
import random
import matplotlib.pyplot as plt
import numpy as np
%load_ext autoreload
%autoreload 2
SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)
plt.rcParams["figure.dpi"] = 90
Exercise 1 — Hand Derivation: Ridge Closed Form (2×2)¶
Given 3 data points with a single feature (no intercept, already centered):
$$ X = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix}, \quad y = \begin{pmatrix} 2 \\ 5 \\ 4 \end{pmatrix} $$
The Ridge objective is $L(\theta) = \frac{1}{n}\|X\theta - y\|_2^2 + \lambda\|\theta\|_2^2$.
Tasks (by hand, then verify with NumPy):
- Compute $X^\top X$ and $X^\top y$.
- Write the Ridge closed form $\hat\theta = (X^\top X + n\lambda I)^{-1} X^\top y$.
- Evaluate $\hat\theta$ for $\lambda = 0$, $\lambda = 1$, $\lambda = 10$.
- Compute the effective degrees of freedom $\text{df}(\lambda) = \sum_j \sigma_j^2 / (\sigma_j^2 + n\lambda)$ for each $\lambda$.
Expected results:
| $\lambda$ | $\hat\theta$ | $\text{df}(\lambda)$ |
|---|---|---|
| 0 | 1.4286 | 1.0000 |
| 1 | 1.2353 | 0.8235 |
| 10 | 0.5946 | 0.3108 |
# Verify your hand calculations here
X_ex1 = np.array([[1], [2], [3]], dtype=float)
y_ex1 = np.array([2, 5, 4], dtype=float)
n_ex1 = len(y_ex1)
# TODO: compute X^T X, X^T y
# XtX = ...
# Xty = ...
# TODO: solve for theta_hat at each lambda
# for lam in [0, 1, 10]:
# theta_hat = ...
# df_lam = ... # hint: use SVD of X
# print(f'lambda={lam:>2} theta={theta_hat[0]:.4f} df={df_lam:.4f}')
Exercise 2 — Conceptual: L1 vs L2 Geometry¶
Consider the 2D constraint sets for Ridge ($\|\theta\|_2^2 \le t$) and Lasso ($\|\theta\|_1 \le t$).
Questions:
Why does the L1 diamond produce exact zeros while the L2 circle does not? (Hint: think about where elliptical contours can be tangent to each shape.)
In dimensions $p > 2$, how many corners does the L1 constraint polytope have? What fraction of those corners lie on a coordinate axis (i.e., have exactly one non-zero coordinate)?
Elastic Net uses $\lambda_1\|\theta\|_1 + \lambda_2\|\theta\|_2^2$. Sketch (or describe) the shape of the Elastic Net constraint set in 2D. Why does it combine sparsity with grouping of correlated features?
Exercise 3 — Coding: Soft-Thresholding Operator¶
The soft-thresholding operator is the proximal map for the L1 norm:
$$ S_\lambda(z) = \text{sign}(z) \max(|z| - \lambda, 0) $$
Tasks:
- Implement
soft_threshold(z, lam)that works on both scalars and arrays. - Verify against these test cases:
| Input $z$ | $\lambda$ | Expected $S_\lambda(z)$ |
|---|---|---|
| 3.0 | 1.0 | 2.0 |
| -3.0 | 1.0 | -2.0 |
| 0.5 | 1.0 | 0.0 |
| -0.5 | 1.0 | 0.0 |
| 0.0 | 1.0 | 0.0 |
| 1.0 | 1.0 | 0.0 |
- Plot $S_\lambda(z)$ for $z \in [-5, 5]$ with $\lambda \in \{0.5, 1.0, 2.0\}$ on the same axes. Include the identity line $y = z$ for reference.
def soft_threshold(z, lam):
"""Soft-thresholding: sign(z) * max(|z| - lam, 0)."""
# TODO: implement (1 line)
pass
Solution 3¶
def soft_threshold(z, lam):
"""Soft-thresholding: sign(z) * max(|z| - lam, 0)."""
return np.sign(z) * np.maximum(np.abs(z) - lam, 0.0)
# Deterministic check
test_cases = [
(3.0, 1.0, 2.0),
(-3.0, 1.0, -2.0),
(0.5, 1.0, 0.0),
(-0.5, 1.0, 0.0),
(0.0, 1.0, 0.0),
(1.0, 1.0, 0.0),
]
for z, lam, expected in test_cases:
result = soft_threshold(z, lam)
assert np.isclose(result, expected, atol=1e-12), (
f'FAIL: S_{lam}({z}) = {result}, expected {expected}'
)
print('All soft-thresholding tests passed.')
All soft-thresholding tests passed.
# TODO: plot S_lambda(z) for lambda in {0.5, 1.0, 2.0}
# z_range = np.linspace(-5, 5, 300)
# ...
# Include identity line y = z as dashed reference
Exercise 4 — Coding: Lasso Recovery on a Known Sparse Signal¶
Generate a dataset with $n = 150$ samples and $p = 20$ features, where only the first 3 features have non-zero true coefficients:
$$ \theta^\ast = (4, -3, 2, 0, 0, \dots, 0)^\top $$
Tasks:
- Generate $X \sim \mathcal{N}(0, 1)$, $y = X\theta^\ast + \epsilon$ with $\epsilon \sim \mathcal{N}(0, 0.5)$.
- Fit Lasso with your from-scratch implementation at $\lambda = 0.05$.
- Verify that the active set matches $\{0, 1, 2\}$ and that the recovered coefficients satisfy $\|\hat\theta_{\text{active}} - \theta^\ast_{\text{active}}\|_\infty < 0.3$.
- Compare with Ridge at the same $\lambda$ — does Ridge zero out any features?
Deterministic check (with SEED = 42):
- Lasso active set =
[0, 1, 2] - Ridge non-zero features = 20 (all features)
# TODO: generate data, fit Lasso and Ridge, compare active sets
# Use the RidgeRegressionScratch and LassoRegressionScratch from
# first_principles.ipynb, or import from src/ml_first_principles/linear_models.py
Exercise 5 — Conceptual: Bayesian Interpretation¶
Questions:
Ridge regression is equivalent to MAP estimation under Gaussian likelihood and a Gaussian prior $\theta \sim \mathcal{N}(0, \tau^2 I)$. Derive the relationship $\lambda = \sigma^2 / (n\tau^2)$, starting from the log-posterior.
If you double the prior variance $\tau^2$ (= weaker prior belief that coefficients are small), what happens to $\lambda$? Does the model become more or less regularized? Explain intuitively.
Lasso corresponds to a Laplace prior. The Laplace density has a cusp at zero. Why does the MAP estimate produce exact zeros, even though the Laplace distribution itself is continuous (i.e., $P(\theta_j = 0) = 0$)?
Exercise 6 — Failure Analysis: When Lasso Misleads¶
Consider two highly correlated features $x_1$ and $x_2$ where $x_2 = x_1 + \epsilon$ with $\epsilon \sim \mathcal{N}(0, 0.01)$, and the true model is $y = x_1 + x_2 + \text{noise}$.
Questions:
Run Lasso on this dataset 20 times with slight noise perturbations. What do you observe about which feature Lasso selects?
What does Elastic Net ($\lambda_1 > 0$, $\lambda_2 > 0$) do differently in this scenario? Why does the L2 component help with grouped selection?
Name one practical scenario where Lasso's instability with correlated features could lead to incorrect scientific conclusions.
# TODO: generate correlated features, run Lasso 20 times,
# observe selection instability
# n = 100
# x1 = rng.normal(size=n)
# x2 = x1 + rng.normal(0, 0.01, size=n)
# X_corr = np.column_stack([x1, x2])
# y_corr = x1 + x2 + rng.normal(0, 0.3, size=n)
# ...