Bias-Variance & Regularization

The goal of supervised learning is not to memorize training labels — any lookup table can do that. The goal is generalization: performing well on new, unseen examples drawn from the same distribution as the training data.

Formally, suppose we draw a training set $S = {(x_i, y_i)}_{i=1}^n$ i.i.d. from some unknown distribution $P(X, Y)$. We learn a hypothesis $h_S$ from $S$ and care about its generalization error (also called true risk or expected loss):

$$R(h_S) = \mathbb{E}_{(x,y) \sim P}[\ell(h_S(x), y)]$$

We can only estimate this from a held-out test set — data not used during training. The train/test gap $R(h_S) - \hat{R}_{\text{train}}(h_S)$ measures how much worse we do on new data versus training data. When this gap is large, the model has overfit.

The key insight: training error is not what we care about. A model that achieves zero training error by memorizing every label (e.g., a lookup table, a degree-1000 polynomial through 10 noisy points) is useless. This is why evaluation methodology — always testing on held-out data — is a sacred rule in ML.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

np.random.seed(42)
X = np.sort(np.random.uniform(0, 1, 20)).reshape(-1, 1)
y = np.sin(2 * np.pi * X).ravel() + np.random.normal(0, 0.2, 20)

X_test = np.linspace(0, 1, 200).reshape(-1, 1)
y_test_true = np.sin(2 * np.pi * X_test).ravel()

for degree in [1, 4, 15]:
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(X, y)
    train_err = np.mean((model.predict(X) - y) ** 2)
    test_err  = np.mean((model.predict(X_test) - y_test_true) ** 2)
    print(f"Degree {degree:2d}: train MSE={train_err:.4f}, test MSE={test_err:.4f}")
# Degree  1: train MSE=0.3226, test MSE=0.3217  <- underfits
# Degree  4: train MSE=0.0388, test MSE=0.0412  <- good balance
# Degree 15: train MSE=0.0024, test MSE=3.1700  <- memorizes noise

These two failure modes sit at opposite ends of the model complexity axis.

Underfitting (high bias): The model is too simple to capture the true structure of the data. A linear model fit to data with a quadratic trend will consistently predict the wrong values regardless of how much data you collect — the error is systematic. Adding more training data will not help; the model family is fundamentally wrong. You will see: high training error AND high test error.

Overfitting (high variance): The model is so flexible that it fits not just the true signal but also the random noise in the training set. A degree-15 polynomial through 20 noisy points will wiggle violently to pass through every sample, failing catastrophically on new points. The model has essentially memorized a particular quirky training set. Adding more training data helps a lot. You will see: low training error BUT high test error.

The diagnostic tool — learning curves: Plot training error and validation error as a function of training set size $n$.

• For an underfitting model: both curves plateau at a high error and meet early. More data does not help.
• For an overfitting model: the curves have a large gap that closes slowly as $n$ grows. The model would do well with more data.

The right level of complexity is where the model captures the true underlying function without fitting noise. In practice, you select complexity via model selection (cross-validation).

# Learning curve diagnostic (pseudocode)
for n in [10, 50, 100, 500, 1000, 5000]:
    X_train, y_train = sample(n points from distribution)
    model.fit(X_train, y_train)
    train_errors.append(mse(model, X_train, y_train))
    val_errors.append(mse(model, X_val, y_val))  # fixed held-out set

plot(n, train_errors, label='Train')
plot(n, val_errors,   label='Validation')
# Underfit: both curves high, close together, flat as n grows
# Overfit:  train low and falling, val high, large gap
# Good fit: train slightly below val, both converging

Here is the central theorem of this module. For squared loss and a fixed test point $x$, the expected test error over all possible training sets $S$ of size $n$ decomposes as:

$$\mathbb{E}_S\left[(y - \hat{f}_S(x))^2\right] = \underbrace{\left(\mathbb{E}_S[\hat{f}S(x)] - f(x)\right)^2}{\text{Bias}^2} + \underbrace{\mathbb{E}_S\left[\left(\hat{f}_S(x) - \mathbb{E}S[\hat{f}S(x)]\right)^2\right]}{\text{Variance}} + \underbrace{\sigma^2}{\text{Noise}}$$

where $f(x) = \mathbb{E}[y|x]$ is the true function and $\sigma^2 = \text{Var}(y|x)$ is irreducible noise.

Bias measures how far the average prediction (averaged over all possible training sets) is from the true answer. A linear model fitting sinusoidal data has high bias — no matter which 20 points you give it, its average prediction will be a line, not a sine curve.

Variance measures how much the model's prediction varies as we change the training set. A degree-15 polynomial has high variance — give it 20 different points and you get a completely different wild curve each time.

Noise $\sigma^2$ is irreducible — it is the inherent randomness in the data-generating process. No model can do better than this on average.

The tradeoff: Simple models have high bias, low variance. Complex models have low bias, high variance. Total error = Bias² + Variance + Noise. We cannot reduce both simultaneously; we minimize their sum.

Proof sketch: Write $y = f(x) + \epsilon$ and let $\bar{f} = \mathbb{E}_S[\hat{f}_S(x)]$. Then:

$$\mathbb{E}[(y - \hat{f})^2] = \mathbb{E}[(f + \epsilon - \hat{f})^2]$$

Add and subtract $\bar{f}$, expand $(a+b+c)^2$, and use independence of $\epsilon$ from the training set to show all cross terms involving $\epsilon$ vanish.

import numpy as np

def empirical_bias_variance(degree, n_trials=500, n_train=20):
    """Empirically estimate bias^2 and variance at test point x=0.5."""
    from sklearn.pipeline import make_pipeline
    from sklearn.preprocessing import PolynomialFeatures
    from sklearn.linear_model import LinearRegression

    x_test = np.array([[0.5]])
    f_true = np.sin(2 * np.pi * 0.5)  # true function value at 0.5
    preds = []
    for _ in range(n_trials):
        X = np.sort(np.random.uniform(0,1, n_train)).reshape(-1,1)
        y = np.sin(2*np.pi*X).ravel() + np.random.normal(0, 0.3, n_train)
        m = make_pipeline(PolynomialFeatures(degree),
                          LinearRegression()).fit(X, y)
        preds.append(m.predict(x_test)[0])
    preds = np.array(preds)
    bias2    = (np.mean(preds) - f_true) ** 2
    variance = np.var(preds)
    print(f"Degree {degree:2d}: Bias^2={bias2:.4f} "
          f"Variance={variance:.4f} Sum={bias2+variance:.4f}")

for d in [1, 3, 9, 15]:
    empirical_bias_variance(d)

Regularization is any technique that reduces model complexity — and thus variance — by adding a constraint or penalty to the learning objective beyond just fitting the training data.

The general framework: instead of minimizing just the empirical risk (training loss), we minimize a regularized objective:

$$\min_\theta \underbrace{\frac{1}{n}\sum_{i=1}^n \ell(h_\theta(x_i), y_i)}{\text{fit the data}} + \underbrace{\lambda \cdot \Omega(\theta)}{\text{complexity penalty}}$$

where $\lambda \geq 0$ is the regularization strength (a hyperparameter) and $\Omega(\theta)$ is a complexity measure on the parameters. Larger $\lambda$ means stronger regularization, pushing toward simpler models.

Why does this work? The penalty prevents the optimizer from going to extreme parameter values that would be needed to fit noise. It introduces a deliberate bias (toward small weights) in exchange for reduced variance. The bias-variance tradeoff is now directly controlled by $\lambda$.

Intuition: Think of the training loss as a landscape with many valleys, some of which are artifact canyons created by noise in the finite training set. The regularization penalty adds a bowl (L2) or a tent (L1) centered at the origin, tilting the landscape so that solutions far from zero are penalized. The optimizer naturally settles in flatter, more central, more generalizable regions.

Hyperparameter selection: $\lambda$ is not learned from training data — it must be selected via cross-validation on held-out data. If you tune $\lambda$ on your test set, you leak test information and your results are invalid.

# Generic regularized training loop (pseudocode)
def train_regularized(X, y, lambda_reg, penalty='L2'):
    theta = initialize_weights_to_zero()
    for step in range(max_steps):
        grad_loss = compute_gradient_of_loss(X, y, theta)
        if penalty == 'L2':
            grad_reg = 2 * lambda_reg * theta      # d/dtheta ||theta||^2
        elif penalty == 'L1':
            grad_reg = lambda_reg * sign(theta)    # subgradient of ||theta||_1
        theta -= learning_rate * (grad_loss + grad_reg)
    return theta

# Effect of lambda on parameter norm (qualitative)
# lambda=0      -> theta = OLS solution, can be arbitrarily large
# lambda=small  -> theta close to OLS but slightly shrunk
# lambda=large  -> theta near zero, high bias low variance
# lambda=inf    -> theta = 0 (predicts mean of y for all inputs)

Ridge regression adds an L2 penalty on the weight vector to ordinary least squares:

$$\hat{\theta}^{\text{ridge}} = \arg\min_\theta |y - X\theta|_2^2 + \lambda|\theta|_2^2$$

Closed-form solution: Ridge has a unique closed-form:

$$\hat{\theta}^{\text{ridge}} = (X^TX + \lambda I)^{-1} X^T y$$

The $\lambda I$ term adds a constant to the diagonal of $X^TX$, which has two beautiful effects: (1) it makes the matrix invertible even when $X^TX$ is singular (more features than samples), and (2) it shrinks the effective singular values, reducing the amplification of noise in low-data directions.

SVD interpretation: Write $X = U\Sigma V^T$. Then:

$$\hat{\theta}^{\text{ridge}} = V,\text{diag}!\left(\frac{\sigma_i^2}{\sigma_i^2 + \lambda}\right) U^T y$$

Each singular direction is shrunk by $\frac{\sigma_i^2}{\sigma_i^2 + \lambda}$. Directions with large singular values (well-supported by data) are shrunk little; directions with small singular values (noisy, data-poor) are shrunk heavily. This is spectral regularization — elegant and interpretable.

Bayesian interpretation: Ridge regression is equivalent to MAP estimation with a Gaussian prior $\theta \sim \mathcal{N}(0, \frac{1}{\lambda}I)$ on the weights. The L2 penalty IS the negative log-prior. This gives a principled probabilistic justification for why L2 makes sense.

Key property: Ridge shrinks all coefficients toward zero but never exactly to zero. It performs shrinkage, not selection. Use this when you believe all features are somewhat relevant.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.datasets import make_regression

X, y, true_coef = make_regression(n_samples=100, n_features=20,
                                   n_informative=5, noise=10,
                                   coef=True, random_state=0)

for lam in [0.001, 1.0, 100.0, 10000.0]:
    ridge = Ridge(alpha=lam).fit(X, y)
    coef_norm = np.linalg.norm(ridge.coef_)
    print(f"lambda={lam:8.3f}: ||theta||_2 = {coef_norm:.4f}")
# lambda=   0.001: ||theta||_2 = 22.09  (near OLS)
# lambda=   1.000: ||theta||_2 = 18.34
# lambda= 100.000: ||theta||_2 =  6.21
# lambda=10000.000: ||theta||_2 =  0.31  (very small, high bias)

def ridge_closed_form(X, y, lam):
    """Numerically stable closed-form ridge regression."""
    n, p = X.shape
    # Use solve instead of explicit matrix inverse
    return np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

def ridge_via_svd(X, y, lam):
    """Ridge via SVD — exposes spectral shrinkage."""
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    d = s**2 / (s**2 + lam)   # shrinkage factors
    return Vt.T @ (d * (U.T @ y))

Lasso (Least Absolute Shrinkage and Selection Operator) replaces the L2 penalty with an L1 penalty:

$$\hat{\theta}^{\text{lasso}} = \arg\min_\theta |y - X\theta|_2^2 + \lambda|\theta|_1$$

where $|\theta|_1 = \sum_j |\theta_j|$.

Lasso has one remarkable property: it drives many coefficients to exactly zero, performing automatic feature selection. This seems like a small difference from ridge, but the implications are enormous.

Why does L1 produce sparsity? The geometric argument: Consider the constrained form (equivalent by Lagrangian duality):

$$\min_\theta |y - X\theta|_2^2 \quad \text{s.t.} \quad |\theta|_1 \leq t$$

The L1 ball ${\theta : |\theta|_1 \leq t}$ in 2D is a diamond (rotated square). Its corners lie on the axes. The least-squares contours (ellipses) expanding outward from their unconstrained minimum will typically first hit the diamond at a corner, where one coordinate is zero. The L2 ball is a sphere — smooth, with no corners — so the first-contact point is generically not on an axis.

The subdifferential argument: At $\theta_j = 0$, the L1 norm is not differentiable; its subdifferential is $[-1, +1]$. The optimality condition says a coefficient stays at zero when the gradient of the loss at that point is small enough to be absorbed by the subgradient: $|\nabla_j \text{loss}|_{\theta_j=0}| \leq \lambda$. For large $\lambda$, many features satisfy this, and their coefficients stay exactly zero.

Bayesian interpretation: Lasso corresponds to MAP estimation with a Laplace prior $\theta_j \sim \text{Laplace}(0, 1/\lambda)$. The Laplace distribution is sharply peaked at zero with heavy tails — it strongly encourages sparsity while allowing nonzero features to have large values.

No closed form: Unlike ridge, lasso requires iterative algorithms like coordinate descent. The update for each coordinate uses the soft-thresholding operator:

$$\hat{\theta}_j \leftarrow \text{sign}(\rho_j) \max(|\rho_j| - \lambda, 0)$$

where $\rho_j = X_j^T(y - X_{-j}\hat{\theta}_{-j}) / n$ is the partial residual correlation for feature $j$.

import numpy as np
from sklearn.linear_model import Lasso, Ridge
from sklearn.datasets import make_regression

X, y, true_coef = make_regression(n_samples=100, n_features=20,
                                   n_informative=5, noise=10,
                                   coef=True, random_state=0)

lasso = Lasso(alpha=10.0).fit(X, y)
ridge = Ridge(alpha=10.0).fit(X, y)

print("Lasso nonzero coefs:", np.sum(lasso.coef_ != 0))  # ~5 (sparse!)
print("Ridge nonzero coefs:", np.sum(ridge.coef_ != 0))  # 20 (all nonzero)

# Soft-thresholding: the core lasso coordinate descent update
def soft_threshold(rho, lam):
    return np.sign(rho) * max(abs(rho) - lam, 0.0)

# Lasso coordinate descent (one full sweep over all features)
def lasso_cd_sweep(X, y, theta, lam):
    n, p = X.shape
    residual = y - X @ theta
    for j in range(p):
        residual  += X[:, j] * theta[j]         # add back feature j
        rho_j      = X[:, j] @ residual / n     # partial residual corr
        theta[j]   = soft_threshold(rho_j, lam) # update with soft threshold
        residual  -= X[:, j] * theta[j]         # remove updated feature j
    return theta

# Full lasso via coordinate descent
def lasso_cd(X, y, lam, tol=1e-6, max_iter=1000):
    n, p = X.shape
    theta = np.zeros(p)
    for _ in range(max_iter):
        theta_old = theta.copy()
        theta = lasso_cd_sweep(X, y, theta, lam)
        if np.max(np.abs(theta - theta_old)) < tol:
            break
    return theta

Elastic Net interpolates between ridge and lasso by including both penalties:

$$\hat{\theta}^{\text{EN}} = \arg\min_\theta |y - X\theta|_2^2 + \lambda_1|\theta|_1 + \lambda_2|\theta|_2^2$$

In scikit-learn's parameterization, this is controlled by alpha (overall strength) and l1_ratio $\rho \in [0,1]$: the L1 penalty is $\rho \cdot \alpha$ and L2 is $(1-\rho) \cdot \alpha$.

When do you need elastic net?

1. Correlated features (grouping effect): Lasso tends to arbitrarily select one feature from a group of highly correlated features and zero out the rest. Elastic net keeps them all (at roughly equal coefficients) because the L2 term prevents any one from dominating. This is critical in genomics where many SNPs are in linkage disequilibrium.

2. The $p \gg n$ setting: When you have far more features than samples, lasso can select at most $n$ features (a limitation from the geometry of the problem). Elastic net can select more because the L2 penalty stabilizes the optimization in this regime.

Geometric intuition: The elastic net constraint region is a rounded diamond — it has the corners of the L1 ball (promoting sparsity) but softened by the L2 term (promoting grouping and stability). It inherits the best practical properties of both.

Applications: Gene expression analysis, financial factor models with correlated predictors, text regression with synonym features.

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

np.random.seed(0)
n = 100
base = np.random.randn(n)         # shared factor
# Features 0,1,2 are nearly identical (correlated triplet)
X = np.column_stack([
    base + 0.01 * np.random.randn(n),
    base + 0.01 * np.random.randn(n),
    base + 0.01 * np.random.randn(n),
    np.random.randn(n),            # independent feature 3
    np.random.randn(n),            # independent feature 4
])
# True signal uses ALL three correlated features equally
y = 2*base + 3*X[:,3] + 0.5*np.random.randn(n)

lasso = Lasso(alpha=0.5).fit(X, y)
en    = ElasticNet(alpha=0.5, l1_ratio=0.5).fit(X, y)

print("Lasso coefs:     ", np.round(lasso.coef_, 2))
# -> e.g. [6.0, 0.0, 0.0, 3.1, 0.0]  picks ONE of triplet
print("ElasticNet coefs:", np.round(en.coef_, 2))
# -> e.g. [2.1, 1.9, 2.0, 2.9, 0.0]  distributes weight evenly

How do you choose $\lambda$? Or the polynomial degree? Or the number of trees? The answer is cross-validation — the gold-standard methodology for model selection.

K-Fold Cross-Validation: Split your training data into $k$ equal folds. For each fold $i$:

1. Train on the remaining $k-1$ folds
2. Evaluate on fold $i$ (the validation fold)
3. Average the $k$ validation errors

This estimates the generalization error of a given hyperparameter choice using only training data.

Why $k = 5$ or $k = 10$? Leave-one-out CV ($k = n$) has lowest bias but highest variance and is computationally expensive. $k = 5$ or $k = 10$ strikes the right empirical balance, as validated in Kohavi (1995) across many datasets.

The right experimental protocol — often violated in practice:

1. Split data into train and test. Lock away the test set.
2. Use cross-validation within the train set to select hyperparameters.
3. Retrain the final model on ALL training data with the selected hyperparameters.
4. Evaluate once on the test set to report generalization performance.

Validation curve: Plot train and CV error vs. $\lambda$. As $\lambda$ increases (stronger regularization), train error rises (constraining the fit) and CV error first falls (reducing overfitting) then rises again (underfitting). The minimum of the CV curve is optimal $\lambda$.

The one-SE rule: Prefer the simplest model (largest $\lambda$) whose CV error is within one standard error of the minimum. This favors interpretability when performance is similar.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score, KFold
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

np.random.seed(0)
X = np.sort(np.random.uniform(0,1, 50)).reshape(-1,1)
y = np.sin(2*np.pi*X).ravel() + np.random.normal(0, 0.2, 50)

# Split: 80% train, 20% test -- TEST SET LOCKED AWAY
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=0)

kf = KFold(n_splits=5, shuffle=True, random_state=42)

best_lam, best_score = None, -np.inf
for lam in np.logspace(-4, 4, 30):
    model = make_pipeline(
        PolynomialFeatures(degree=9),
        StandardScaler(),
        Ridge(alpha=lam)
    )
    # negative MSE (sklearn: higher=better)
    scores = cross_val_score(model, X_train, y_train, cv=kf,
                             scoring='neg_mean_squared_error')
    if scores.mean() > best_score:
        best_score = scores.mean()
        best_lam   = lam

print(f"Best lambda: {best_lam:.4f}")
# Retrain on ALL training data
final_model = make_pipeline(
    PolynomialFeatures(9), StandardScaler(), Ridge(best_lam)
).fit(X_train, y_train)
# ONE final evaluation on locked test set
test_mse = np.mean((final_model.predict(X_test) - y_test)**2)
print(f"Final test MSE: {test_mse:.4f}")

The regularization path traces how the coefficients $\hat{\theta}(\lambda)$ change as $\lambda$ varies from 0 (no regularization, full OLS) to $\infty$ (all coefficients zero). This gives a rich picture of your data.

Ridge path: All coefficients shrink smoothly and continuously toward zero as $\lambda$ increases. No coefficient ever reaches exactly zero. The path is smooth and differentiable — a direct consequence of the smooth L2 ball.

Lasso path: Coefficients take piecewise-linear paths as a function of $\lambda$. Some features enter the model (leave zero) at specific $\lambda$ values, and others may return to zero at different $\lambda$ values. This is exploited by the LARS (Least Angle Regression) algorithm, which traces the entire lasso path efficiently in $O(np \min(n,p))$ time — the same cost as a single OLS fit.

Reading the lasso path as data storytelling: The features that remain nonzero for the largest range of $\lambda$ are the most robustly predictive. The first feature to enter as $\lambda$ decreases from $\infty$ is the single most correlated feature with the response.

The one-SE rule on the path: After computing the path and running CV, many practitioners choose not the minimum CV-error $\lambda$ but the largest $\lambda$ within one standard error of the minimum. This prefers sparser models when performance is similar, improving interpretability without significantly hurting accuracy.

import numpy as np
from sklearn.linear_model import lasso_path
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=100, n_features=10,
                        n_informative=4, noise=10, random_state=0)
# Always standardize features before lasso (L1 is not scale-invariant)
X = (X - X.mean(0)) / X.std(0)

# Compute the entire regularization path efficiently
alphas, coefs, _ = lasso_path(X, y, n_alphas=100)
# coefs shape: (n_features=10, n_alphas=100)
# coefs[:, 0]  -> all zeros (max lambda, full regularization)
# coefs[:, -1] -> near-OLS solution (min lambda, no regularization)

print(f"Alpha range: [{alphas[-1]:.3f}, {alphas[0]:.3f}]")
print(f"Nonzero at max lambda: {np.sum(coefs[:,0] != 0)}")
print(f"Nonzero at min lambda: {np.sum(coefs[:,-1] != 0)}")

# Order in which features enter the model (first nonzero)
for j in range(coefs.shape[0]):
    first_nonzero = np.argmax(coefs[j,:] != 0)   # index in path
    if coefs[j,:].any():
        print(f"Feature {j} enters at alpha={alphas[first_nonzero]:.3f}")

Everything we have learned about L1/L2 regularization applies directly to deep learning, though it appears in different forms.

Weight decay is L2 regularization on all network weights. In PyTorch, passing weight_decay=1e-4 to an optimizer adds $\lambda |\theta|_2^2$ to the effective loss. This is the default regularization in most deep learning recipes.

Dropout (Srivastava et al. 2014) is a very different form of regularization: randomly zero out activations during training with probability $p$. At test time, scale activations by $(1-p)$ (inverted dropout). Dropout can be interpreted as training an exponential ensemble of thinned networks sharing weights. It is empirically one of the most effective regularizers for deep networks.

Early stopping is perhaps the most widely used regularization: stop training when the validation error stops decreasing, even if training error is still falling. For gradient descent on convex objectives, early stopping is mathematically equivalent to L2 regularization — the number of steps controls the effective $\lambda$. This connection extends approximately to deep networks.

Data augmentation (random crops, flips, color jitter) increases the effective training set size, reducing variance. It is implicit regularization with no parameter to tune.

The modern picture — double descent: The classical U-shaped bias-variance curve holds for fixed model families. For overparameterized neural networks (parameters $\gg$ training samples), test error can continue decreasing past the interpolation threshold as model size grows. This "benign overfitting" occurs because large overparameterized models fit data via minimum-norm interpolation, which generalizes well when signal is aligned with large-singular-value directions. Regularization (weight decay, dropout) remains crucial to achieving this in practice.

import torch
import torch.nn as nn

class RegularizedNet(nn.Module):
    def __init__(self, dropout_p=0.3):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(20, 64),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),   # zeros activations w.p. p at train time
            nn.Linear(64, 64),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(64, 1)
        )
    def forward(self, x):
        return self.net(x)

model = RegularizedNet(dropout_p=0.3)

# L2 weight decay via optimizer (applied to all weight parameters)
optimizer = torch.optim.Adam(
    model.parameters(),
    lr=1e-3,
    weight_decay=1e-4   # lambda for L2 penalty
)

# Early stopping loop
best_val_loss = float('inf')
patience, wait = 10, 0
for epoch in range(1000):
    model.train()   # dropout ACTIVE during training
    # ... forward pass, loss, optimizer.step() ...

    model.eval()    # dropout DISABLED at evaluation
    with torch.no_grad():
        val_loss = evaluate_val_set(model)  # placeholder
    if val_loss < best_val_loss - 1e-5:
        best_val_loss = val_loss
        torch.save(model.state_dict(), 'checkpoint.pt')
        wait = 0
    else:
        wait += 1
        if wait >= patience:
            print(f'Early stopping at epoch {epoch}')
            model.load_state_dict(torch.load('checkpoint.pt'))
            break