III-12 · Regularization

Expose the price of a coefficient

Compare Ridge and Lasso on the same observations, inspect their real penalty terms, and choose whether the intercept is penalized.

Fit and optimization laboratory

Fit error + penalty → shrinkage

Hand-authored toy observations, editable in this sheet. No seed and no measured population.

● Training observation · ◇ Prediction · dashed residual. Axis domain adjusts to the data; x/y units share one scale.

Change the active experimentPenaltyPenalty strength λ0 to 10 · Enter to applyPenalize the intercept

Fit satisfies the KKT tolerance ≤ 1e−8.

Selected observation

Observation x-5 to 5 · Enter to applyObservation target-5 to 5 · Enter to apply

X · augmented designSamples in rows; first column is 1.: 5 × 2
W · coefficientsIntercept first, then feature weights.: 2 × 1
y, predictions, residuals: 5 × 1
NormalizationNo half factor; evaluation excludes penalties.: MSE / n

W = [≈ 0.8, ≈ 1.44]ᵀ

Live state

Data MSE: ≈ 0.3392
Penalty term: ≈ 1.0368
Total objective: ≈ 1.376
KKT residual: ≈ 6.6613e−16

\hat y=XW,\quad r=\hat y-y

L_{\mathrm{data}}=\frac1{5}\sum_i r_i^2\approx 0.3392

J=L_{\mathrm{data}}+\lambda\sum_{j\in P}w_j^2\approx 1.376

Exact inputs; approximate computed results from the active state
ID	x	Design row	Target	Prediction	Residual
s1	-2	[1, -2]	-3	≈ −2.08	≈ 0.92
s2	-1	[1, -1]	-1	≈ −0.64	≈ 0.36
s3	0	[1, 0]	1	≈ 0.8	≈ −0.2
s4	1	[1, 1]	3	≈ 2.24	≈ −0.76
s5	2	[1, 2]	4	≈ 3.68	≈ −0.32

Same observations, lambda and intercept policy for both penalties
Penalty	b	w	Data MSE	Penalty	Objective
L2	≈ 0.8	≈ 1.44	≈ 0.3392	≈ 1.0368	≈ 1.376
L1	≈ 0.8	≈ 1.675	≈ 0.1113	≈ 0.8375	≈ 0.9487

Change the controls first, then inspect the active equation, numeric rows and shapes. Exact inputs and approximate computed values are distinguished; Reset restores the hand-authored experiment.

Watch shrinkage change the actual objective

Switch between L2/Ridge and L1/Lasso and change $\lambda$ . The predictions, residuals, coefficients and penalty value update together. The comparison table fits both penalties to the same observations with the same intercept policy. Check Penalize the intercept to make that decision explicit; its default is off.

At $\lambda=0$ , the penalty disappears and the ordinary least-squares solution returns. Raising $\lambda$ trades fit to these training targets for smaller coefficients. A small training objective does not by itself establish lower evaluation error; the displayed data MSE deliberately excludes the penalty.

Define the loss, penalty and policy separately

Use augmented design $X=[\mathbf1,x]$ , coefficient column $W=[b,w]^\mathsf T$ and data loss $L_{\mathrm{data}}=\|XW-y\|^2/n$ . Let $P$ contain the parameter indices to penalize. By default $P$ excludes the intercept.

J_{L2}(W)=L_{\mathrm{data}}+\lambda\sum_{j\in P}w_j^2, \qquad J_{L1}(W)=L_{\mathrm{data}}+\lambda\sum_{j\in P}|w_j|.

Neither definition has a half factor. Consequently the Ridge gradient contribution is $2\lambda w_j$ , whereas the L1 nonzero-coordinate subgradient is $\lambda\operatorname{sign}(w_j)$ . At zero, the L1 subgradient is the interval $[-\lambda,\lambda]$ . Changing normalization changes the numerical interpretation of $\lambda$ .

The intercept shifts the prediction origin. Penalizing it pulls that origin toward zero, so whether it belongs to $P$ is a modeling choice. Excluding it means a constant offset in all targets can move the intercept without incurring a direct penalty.

Derive both fitting mechanisms

For L2, append synthetic rows to the design. A penalized coefficient gets a row with $\sqrt{n\lambda}$ in its position and target zero. Then

\frac1n\left\|\begin{bmatrix}X\\\sqrt{n\lambda}D\end{bmatrix}W- \begin{bmatrix}y\\0\end{bmatrix}\right\|^2 =L_{\mathrm{data}}+\lambda\|DW\|^2,

where $D$ selects penalized coordinates. The rank-aware SVD solves this augmented problem; it never inverts normal equations. Original observations still determine $n$ in the objective, not the number of augmented rows.

For L1, hold all other coefficients fixed and let the partial residual exclude coordinate $j$ . Define $\rho=(1/n)\sum_i x_{ij}(y_i-\sum_{k\ne j}x_{ik}w_k)$ and $a=(1/n)\sum_i x_{ij}^2$ . Completing the one-coordinate quadratic plus absolute-value penalty yields

w_j^+=\frac{\operatorname{sign}(\rho)\max(|\rho|-\lambda/2,0)}a

for a penalized nonzero-design coordinate. An unpenalized coordinate uses threshold zero; an all-zero feature takes representative coefficient zero. Coordinate descent repeats these exact coordinate minimizations. The sheet reports the KKT residual, making an approximate iteration-limited solution distinguishable from a converged one.

Work an orthogonal two-row example

For inputs $[-1,1]$ , targets $[3,7]$ and unpenalized intercept, the data loss is

L_{\mathrm{data}}=(b-5)^2+(w-2)^2.

At $\lambda=2$ , Ridge gives $b=5,w=2/(1+2)=2/3$ . Data MSE is $16/9$ , penalty $2(2/3)^2=8/9$ , and objective $8/3$ . Lasso gives $b=5,w=\operatorname{soft}(2,1)=1$ . Its data MSE is $1$ , penalty $2$ , and objective $3$ . L1 can make an exact zero: at $\lambda=4$ , $w=0$ .

If the intercept is also penalized at $\lambda=2$ , Ridge gives $b=5/3$ and Lasso gives $b=4$ . These are not the predictions of the default policy. The checkbox changes the objective and recomputes the fit.

The opening five-row data are the hand-authored III-07 observations. Their unpenalized optimum is $b=0.8,w=1.8$ . At default $\lambda=0.5$ , Ridge gives $w=1.44$ and Lasso $w=1.675$ , both with $b=0.8$ . Try these numbers in the comparison table. There is no seed or held-out set here.

KKT, shapes and interpretation

For L1 at a zero penalized coefficient, $|\partial_j L_{\mathrm{data}}|\le\lambda$ is the optimality condition; at nonzero coefficients use $\partial_jL_{\mathrm{data}}+\lambda\operatorname{sign}(w_j)=0$ . Unpenalized coordinates require a zero data gradient. L2 uses its smooth total gradient. The maximum violation is the shown KKT residual, with acceptance tolerance $10^{-8}$ .

Samples remain rows, $X$ is $n\times2$ , $W$ is $2\times1$ and predictions are $n\times1$ . Synthetic penalty rows are a solving device, not additional observations or an evaluation dataset. The pure implementation supports up to 64 original rows, including Ridge’s internal augmentation.

Mistakes, limits and connections

Do not compare numerical L1 and L2 objectives as if their penalties were interchangeable units. Do not add a penalty to evaluation MSE. Feature scaling changes what a coefficient magnitude means; standardization learned from training data matters when comparing different feature units.

The sheet accepts $0\le\lambda\le10$ . L1 coordinate descent has a 3000-sweep budget and parameter-change threshold $10^{-12}$ , but the final KKT check independently decides convergence. Correlated features can slow it; the UI labels a budget-limited approximate solution. Uniqueness is not promised for rank-deficient L1 or zero-penalty designs. Finite arithmetic is validated; overflow preserves the last supported state. NumPy’s canonical regularized uses these mechanisms; separate analytic shrinkage fixtures, Ridge lstsq checks and gradient/KKT tests verify them.

Regularization changes the hypothesis preference, not the observation truth. It connects to Book II’s constrained optimization and MAP interpretation; choosing $\lambda$ using held-out evaluation belongs to later chapters. You should be able to write both penalties, account for the intercept policy, derive soft thresholding and distinguish data loss from total objective.

Execute the actual Python and NumPy

The imported source below is the runnable file used in canonical comparisons. Read the named mechanism from this chapter; the other functions coordinate the same Part experiment. From the book root:

.venv/bin/python -B code/python/p02/fit_reference.py

code/python/p02/fit_reference.py

"""III-08–12: centered fits, frozen regression updates, small SVD and explicit penalties.

Run from the book root: .venv/bin/python -B code/python/p02/fit_reference.py
Rows are samples; W and y are columns. Intercept is first in augmented X/W.
"""
import json
import math
import numpy as np
from linear_regression import CANONICAL_REGRESSION, derive_regression, regression_trace

RCOND = 1e-10
POLYNOMIAL_DATA = {
    "trainX": [-1, -.6, -.2, .2, .6, 1], "trainY": [1.1, .46, 0, .16, .26, 1.2],
    "heldX": [-.8, -.4, 0, .4, .8], "heldY": [.64, .16, 0, .16, .64],
}


def finite(values, name):
    if not np.isfinite(values).all():
        raise ValueError(f"{name} must remain finite")
    return values


def design(X, y):
    # Validate before coercion: booleans, strings and ragged rows are not numeric data.
    if not isinstance(X, (list, tuple)) or not 1 <= len(X) <= 64:
        raise ValueError("X requires 1–64 rows")
    if not isinstance(X[0], (list, tuple)) or not 1 <= len(X[0]) <= 7:
        raise ValueError("X requires 1–7 columns")
    d = len(X[0])
    for row in X:
        if not isinstance(row, (list, tuple)) or len(row) != d:
            raise ValueError("X must be rectangular")
        for v in row:
            if isinstance(v, bool) or not isinstance(v, (int, float)) or not math.isfinite(v):
                raise ValueError("X must contain finite numbers")
    column(y, len(X), "Targets")
    return finite(np.array(X, dtype=float), "X"), np.array(y, dtype=float).reshape(-1, 1)


def column(values, size, name):
    if not isinstance(values, (list, tuple)) or len(values) != size:
        raise ValueError(f"{name} requires {size} entries (a column)")
    if any(isinstance(v, bool) or not isinstance(v, (int, float)) or not math.isfinite(v) for v in values):
        raise ValueError(f"{name} must contain finite numbers")
    return np.array(values, dtype=float).reshape(-1, 1)


def linear(X, y, weights):
    A, target = design(X, y)
    W = column(weights, A.shape[1], "Weights")
    with np.errstate(over="ignore", invalid="ignore", under="ignore"):
        prediction = finite(A @ W, "Predictions")
        residual = finite(prediction - target, "Residuals")
        mse = float(finite(np.sum(residual * residual) / len(A), "MSE"))
        gradient = finite(2 / len(A) * A.T @ residual, "Gradient")
    n, d = A.shape
    return {"predictions": prediction[:, 0].tolist(), "residuals": residual[:, 0].tolist(),
            "mse": mse, "gradient": gradient[:, 0].tolist(),
            "shapes": {"X": [n, d], "W": [d, 1], "y": [n, 1], "prediction": [n, 1], "gradient": [d, 1]}}


def _least_norm_validated(A, target):
    # Trusted Ridge augmentation may add ≤7 rows beyond the 64 public observation rows.
    d = A.shape[1]
    scale = float(np.max(np.abs(A)))
    if scale == 0:
        return {"weights": [0.] * d, "rank": 0, "singularValues": [0.] * d, "cutoff": 0., "condition": None}
    B = A / scale
    V = np.eye(d)
    for sweep in range(100):
        rotated = False
        for p in range(d):
            for q in range(p + 1, d):
                a, b, c = float(B[:, p] @ B[:, p]), float(B[:, q] @ B[:, q]), float(B[:, p] @ B[:, q])
                if a == 0 or b == 0 or abs(c) <= 1e-14 * math.sqrt(a) * math.sqrt(b):
                    continue
                tau = (b - a) / (2 * c)
                t = 1. if tau == 0 else math.copysign(1., tau) / (abs(tau) + math.hypot(1., tau))
                cs, sn = 1 / math.hypot(1., t), t / math.hypot(1., t)
                for matrix in (B, V):
                    left, right = matrix[:, p].copy(), matrix[:, q].copy()
                    matrix[:, p], matrix[:, q] = cs * left - sn * right, sn * left + cs * right
                rotated = True
        if not rotated:
            break
    else:
        raise ValueError("SVD did not converge in 100 sweeps; no fit returned")
    norms = np.array([math.hypot(*B[:, j]) for j in range(d)])
    cutoff = RCOND * float(np.max(norms))
    weights = np.zeros(d)
    for j in range(d):
        if norms[j] > cutoff:
            coefficient = float((B[:, j] / norms[j]) @ target[:, 0]) / norms[j] / scale
            weights += V[:, j] * coefficient
    finite(weights, "Fitted weights")
    residual = finite(A @ weights.reshape(-1, 1) - target, "Fitted residuals")
    finite(np.sum(residual * residual) / len(A), "Fitted MSE")
    values = sorted(finite(norms * scale, "Singular values").tolist(), reverse=True)
    rank = int(np.sum(norms > cutoff))
    condition = values[0] / values[-1] if rank == d else None
    return {"weights": weights.tolist(), "rank": rank, "singularValues": values,
            "cutoff": cutoff * scale, "condition": condition}


def least_norm(X, y):
    A, target = design(X, y)
    return _least_norm_validated(A, target)


def closed_form(samples):
    derive_regression({"samples": samples, "weight": 0., "bias": 0., "learningRate": .1})
    x = np.array([r["x"] for r in samples], dtype=float)
    y = np.array([r["y"] for r in samples], dtype=float)
    with np.errstate(over="ignore", invalid="ignore", under="ignore"):
        mean_x, mean_y = float(finite(np.mean(x), "Mean x")), float(finite(np.mean(y), "Mean y"))
        dx, dy = finite(x - mean_x, "Centered x"), finite(y - mean_y, "Centered y")
        sxx, sxy = float(finite(dx @ dx, "Sxx")), float(finite(dx @ dy, "Sxy"))
    if sxx == 0 and any(v != x[0] for v in x):
        raise ValueError("Input variance underflowed; slope is numerically unresolved")
    weight = 0. if sxx == 0 else sxy / sxx
    bias = mean_y - weight * mean_x
    mse = derive_regression({"samples": samples, "weight": weight, "bias": bias, "learningRate": .1})["mse"]
    return {"weight": weight, "bias": bias, "meanX": mean_x, "meanY": mean_y, "sxx": sxx, "sxy": sxy, "mse": mse, "unique": sxx != 0}


def descent(state):
    result = derive_regression(state)  # Frozen III-07 is the predictor, gradient and update authority.
    x = np.array([r["x"] for r in state["samples"]])
    a, b, c = 2 * float(np.mean(x * x)), 2 * float(np.mean(x)), 2.
    maximum = (a + c + math.hypot(a - c, 2 * b)) / 2
    factor = 1 - state["learningRate"] * maximum
    behavior = ("undamped oscillation" if factor == -1 else "divergence") if factor <= -1 else ("convergent oscillation" if factor < 0 else "convergence without fast-mode sign flips")
    return {**result, "lambdaMax": maximum, "criticalRate": 2 / maximum, "factor": factor, "behavior": behavior}


def polynomial_design(x, degree):
    if type(degree) is not int or not 1 <= degree <= 5:
        raise ValueError("Degree must be an integer from 1 to 5")
    column(x, len(x), "Polynomial input")
    return finite(np.array([[v ** j for j in range(degree + 1)] for v in x]), "Polynomial features").tolist()


def polynomial_fit(samples, degree):
    closed_form(samples)
    X, y = polynomial_design([r["x"] for r in samples], degree), [r["y"] for r in samples]
    solve = least_norm(X, y)
    return {**solve, "train": linear(X, y, solve["weights"]), "held": linear(polynomial_design(POLYNOMIAL_DATA["heldX"], degree), POLYNOMIAL_DATA["heldY"], solve["weights"]), "X": X}


def regularized(X, y, lam, kind, penalize_intercept):
    A, target = design(X, y)
    if isinstance(lam, bool) or not isinstance(lam, (int, float)) or not math.isfinite(lam) or not 0 <= lam <= 10 or kind not in ("l1", "l2") or type(penalize_intercept) is not bool:
        raise ValueError("Penalty requires L1/L2, lambda in [0,10] and an explicit intercept policy")
    if not np.all(A[:, 0] == 1):
        raise ValueError("First design column must be the intercept ones")
    n, d = A.shape
    mask = np.array([j != 0 or penalize_intercept for j in range(d)])
    iterations = 0
    if kind == "l2" or lam == 0:
        extra = np.diag(np.full(d, math.sqrt(n * lam)))[mask] if lam > 0 else np.empty((0, d))
        solve = _least_norm_validated(np.vstack([A, extra]), np.vstack([target, np.zeros((len(extra), 1))]))
        weights = np.array(solve["weights"])
    else:
        weights = np.zeros(d)
        for iterations in range(1, 3001):
            old = weights.copy()
            for j in range(d):
                partial = target[:, 0] - A @ weights + A[:, j] * weights[j]
                rho, curvature = float(A[:, j] @ partial / n), float(A[:, j] @ A[:, j] / n)
                threshold = lam / 2 if mask[j] else 0
                weights[j] = 0 if curvature == 0 else math.copysign(max(0, abs(rho) - threshold), rho) / curvature
            finite(weights, "L1 weights")
            if np.max(np.abs(weights - old)) < 1e-12:
                break
    result = linear(X, y, weights.tolist())
    penalty = float(finite(lam * np.sum((weights ** 2 if kind == "l2" else np.abs(weights))[mask]), "Penalty"))
    gradient = np.array(result["gradient"])
    kkt = []
    for j, w in enumerate(weights):
        g = gradient[j]
        kkt.append(abs(g) if not mask[j] else abs(g + 2 * lam * w) if kind == "l2" else max(0, abs(g) - lam) if w == 0 else abs(g + lam * np.sign(w)))
    residual = float(max(kkt))
    return {**result, "weights": weights.tolist(), "penalty": penalty, "objective": float(finite(result["mse"] + penalty, "Objective")), "kkt": residual, "iterations": iterations, "converged": residual <= 1e-8}


def demonstration():
    samples = CANONICAL_REGRESSION["samples"]
    poly = [dict(id=f"s{i+1}", x=x, y=POLYNOMIAL_DATA["trainY"][i]) for i, x in enumerate(POLYNOMIAL_DATA["trainX"])]
    return {"closedForm": closed_form(samples), "descent": descent(CANONICAL_REGRESSION),
            "descentTrace": regression_trace(CANONICAL_REGRESSION, 4),
            "multiple": least_norm([[1, -1, 0], [1, 0, 1], [1, 1, 0], [1, 2, 1]], [-1, 0, 3, 4]),
            "polynomial": {str(d): polynomial_fit(poly, d) for d in (1, 2, 5)},
            "regularization": {kind: regularized([[1, r["x"]] for r in samples], [r["y"] for r in samples], .5, kind, False) for kind in ("l1", "l2")}}


if __name__ == "__main__":
    print(json.dumps(demonstration(), indent=2, allow_nan=False))