Tikhonov Regularisation — L2 Penalty and Bias-Variance Trade-off

Companion notebook to the Tikhonov Regularisation article (Regularisation & Stability, Module 2).

§1 — Filter factor computation and selective suppression of small singular values
§2 — Bias-Variance decomposition: Monte Carlo verification of Theorem 2.2
§3 — L-curve generation and corner detection for λ selection
§4 — GCV minimisation vs. oracle optimal λ
§5 — Validation summary

💡 Run cells in order — each cell uses variables from the previous ones.

§0 — Imports

import numpy as np

RNG = np.random.default_rng(42)
print("✓ NumPy", np.__version__, "ready.")

§1 — Filter Factors: Selective Suppression

The Tikhonov filter factor for singular value σᵢ is fᵢ(λ) = σᵢ²/(σᵢ² + λ²). We verify that: (a) for σᵢ ≫ λ, fᵢ ≈ 1 (well-identified, unfiltered); (b) for σᵢ ≪ λ, fᵢ → 0 (ill-conditioned, suppressed); (c) the crossover σᵢ = λ gives fᵢ = 0.5 exactly.

def filter_factor(sigma: float, lam: float) -> float:
    """Tikhonov filter factor f_i(lambda) = sigma^2 / (sigma^2 + lambda^2)."""
    return sigma**2 / (sigma**2 + lam**2)

# Singular value spectrum: from large (well-identified) to small (ill-conditioned)
sigmas = [10.0, 1.0, 0.1, 0.01, 0.001]
lambdas = [0.001, 0.01, 0.1, 1.0, 10.0]

print(f"{'sigma':>8} | " + " | ".join(f"λ={lam:>5}" for lam in lambdas))
print("-" * 60)
for sigma in sigmas:
    row = f"{sigma:8.3f} | "
    row += " | ".join(f"{filter_factor(sigma, lam):8.4f}" for lam in lambdas)
    print(row)

print()
# Verify crossover: f(sigma, lambda=sigma) = 0.5
for sigma in sigmas:
    f = filter_factor(sigma, sigma)
    print(f"f(σ={sigma}, λ=σ) = {f:.4f}  (expected 0.5000) ✓" if abs(f - 0.5) < 1e-10 else f"FAIL: f={f}")

print()
print("→ For λ = 0.1: singular values σ >> 0.1 are unfiltered; σ << 0.1 are suppressed.")
print("→ Tikhonov is a smooth threshold at σ = λ.")

§2 — Bias-Variance Decomposition (Theorem 2.2)

We construct a synthetic ill-posed problem with known singular values and true parameters, then run Monte Carlo to estimate Bias² + Variance and compare to the theoretical formula.

def make_jacobian_prescribed(n_obs: int, sigmas: list, seed: int = 0) -> np.ndarray:
    """Build n_obs x p Jacobian with prescribed singular values."""
    p = len(sigmas)
    rng = np.random.default_rng(seed)
    U_raw = rng.standard_normal((n_obs, p))
    U, _ = np.linalg.qr(U_raw)
    U = U[:, :p]
    V, _ = np.linalg.qr(rng.standard_normal((p, p)))
    Sigma = np.diag(sigmas)
    return U @ Sigma @ V.T, np.diag(sigmas), V

def tikhonov_solution(J: np.ndarray, y: np.ndarray, lam: float) -> np.ndarray:
    """Tikhonov regularised least-squares: (J'J + lam^2 I)^{-1} J' y."""
    p = J.shape[1]
    A = J.T @ J + lam**2 * np.eye(p)
    b = J.T @ y
    return np.linalg.solve(A, b)

N_OBS = 20
SIGMAS = [3.0, 1.0, 0.1, 0.01]   # 4 parameters, decreasing identifiability
DELTA = 0.05                       # noise standard deviation
N_MC = 2000

J, _, V = make_jacobian_prescribed(N_OBS, SIGMAS)
theta_star = np.array([1.0, 0.5, 0.3, 0.1])   # true parameters

lambdas_test = [0.001, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0]

print(f"{'lambda':>8} | {'MC MSE':>12} | {'Theory MSE':>12} | {'Bias^2':>10} | {'Variance':>10}")
print("-" * 60)

for lam in lambdas_test:
    # Monte Carlo
    theta_estimates = np.zeros((N_MC, 4))
    for trial in range(N_MC):
        eps = RNG.standard_normal(N_OBS) * DELTA
        y_noisy = J @ theta_star + eps
        theta_estimates[trial] = tikhonov_solution(J, y_noisy, lam)

    mc_mse = np.mean(np.sum((theta_estimates - theta_star)**2, axis=1))

    # Theoretical Bias^2 + Variance (SVD formula, Theorem 2.2)
    _, s, Vt = np.linalg.svd(J, full_matrices=False)
    theta_v = Vt @ theta_star     # coordinates in V-basis
    bias_sq = np.sum((lam**2 / (s**2 + lam**2))**2 * theta_v**2)
    variance = DELTA**2 * np.sum(s**2 / (s**2 + lam**2)**2)
    theory_mse = bias_sq + variance

    print(f"{lam:8.3f} | {mc_mse:12.4e} | {theory_mse:12.4e} | {bias_sq:10.4e} | {variance:10.4e}")

print()
print("→ MC MSE closely tracks theoretical MSE confirming Theorem 2.2.")

§3 — L-Curve Generation and Corner Detection

The L-curve plots ‖r(θ_λ)‖ vs ‖θ_λ‖ on a log-log scale for a range of λ. The optimal λ is at the corner — where curvature is maximised. We approximate curvature using finite differences of the parametric curve.

✦ The corner of the L-curve sits near the optimal λ that minimises the true MSE. This works without knowing the noise level δ — a key practical advantage.

def lcurve_point(J: np.ndarray, y: np.ndarray, lam: float):
    """Returns (log‖r‖, log‖θ‖) for the Tikhonov solution at regularisation λ."""
    theta_lam = tikhonov_solution(J, y, lam)
    r_norm = np.linalg.norm(y - J @ theta_lam)
    theta_norm = np.linalg.norm(theta_lam)
    return np.log10(r_norm + 1e-15), np.log10(theta_norm + 1e-15)

def lcurve_curvature(log_r_arr, log_t_arr):
    """
    Approximate curvature of the L-curve at each interior point.
    Uses the formula κ = |x'y'' - y'x''| / (x'^2 + y'^2)^(3/2)
    where x = log‖r‖, y = log‖θ‖ parameterised by index.
    """
    dx = np.gradient(log_r_arr)
    dy = np.gradient(log_t_arr)
    d2x = np.gradient(dx)
    d2y = np.gradient(dy)
    num = np.abs(dx * d2y - dy * d2x)
    denom = (dx**2 + dy**2)**1.5 + 1e-30
    return num / denom

# One noisy realisation
y_obs = J @ theta_star + RNG.standard_normal(N_OBS) * DELTA

lambdas_lc = np.logspace(-4, 2, 80)
log_r = np.zeros(len(lambdas_lc))
log_t = np.zeros(len(lambdas_lc))

for k, lam in enumerate(lambdas_lc):
    log_r[k], log_t[k] = lcurve_point(J, y_obs, lam)

# Find the corner (max curvature)
curvature = lcurve_curvature(log_r, log_t)
corner_idx = np.argmax(curvature)
lambda_corner = lambdas_lc[corner_idx]

# Also find the oracle optimal lambda (minimises true MSE)
oracle_mses = []
for lam in lambdas_lc:
    theta_lam = tikhonov_solution(J, y_obs, lam)
    oracle_mses.append(np.sum((theta_lam - theta_star)**2))
lambda_oracle = lambdas_lc[np.argmin(oracle_mses)]

print(f"L-curve results (one noise realisation, δ = {DELTA}):")
print(f"  L-curve corner    λ* = {lambda_corner:.4f}  (log₁₀ = {np.log10(lambda_corner):.2f})")
print(f"  Oracle optimal    λ* = {lambda_oracle:.4f}  (log₁₀ = {np.log10(lambda_oracle):.2f})")
print(f"  Ratio oracle/corner: {lambda_oracle/lambda_corner:.2f}  (should be O(1))")
print()
# Print the L-curve around the corner
print(f"{'lambda':>10} | {'log10‖r‖':>10} | {'log10‖θ‖':>10} | {'curvature':>12} | corner?")
print("-" * 60)
for k in range(max(0, corner_idx-3), min(len(lambdas_lc), corner_idx+4)):
    marker = " ← CORNER" if k == corner_idx else ""
    print(f"{lambdas_lc[k]:10.4f} | {log_r[k]:10.3f} | {log_t[k]:10.3f} | {curvature[k]:12.4f} |{marker}")

§4 — Generalised Cross-Validation (GCV) vs. Oracle λ

GCV minimises N·‖r(θ_λ)‖² / [tr(I − H_λ)]² where H_λ = J(J'J + λ²I)⁻¹J' is the hat matrix. We compute the GCV score in the SVD basis (efficient) and compare the GCV-selected λ to the oracle.

def gcv_score(J: np.ndarray, y: np.ndarray, lam: float) -> float:
    """
    GCV score for linear Tikhonov: numerator = ‖(I - H_λ)y‖², denominator = [tr(I - H_λ)]².
    Computed efficiently in the SVD basis.
    """
    U, s, Vt = np.linalg.svd(J, full_matrices=True)
    p = len(s)
    n = J.shape[0]

    # Residuals in U basis
    Uty = U.T @ y
    # Filter factors
    f = s**2 / (s**2 + lam**2)     # shape (p,)

    # Numerator: ‖(I - H_λ)y‖² = Σᵢ (1 - fᵢ)² (uᵢ'y)² + Σᵢ>p (uᵢ'y)²
    resid_sq = np.sum(((1 - f) * Uty[:p])**2) + np.sum(Uty[p:]**2)
    # Denominator: tr(I - H_λ) = Σᵢ (1 - fᵢ) + (n - p)
    dof = np.sum(1 - f) + (n - p)
    return resid_sq / (dof**2 + 1e-30)

lambdas_gcv = np.logspace(-4, 2, 120)
gcv_scores = np.array([gcv_score(J, y_obs, lam) for lam in lambdas_gcv])
lambda_gcv = lambdas_gcv[np.argmin(gcv_scores)]

print(f"GCV results (same noise realisation):")
print(f"  GCV-selected  λ* = {lambda_gcv:.4f}  (log₁₀ = {np.log10(lambda_gcv):.2f})")
print(f"  Oracle optimal λ* = {lambda_oracle:.4f}  (log₁₀ = {np.log10(lambda_oracle):.2f})")
print(f"  Ratio oracle/GCV: {lambda_oracle/lambda_gcv:.2f}  (should be O(1))")
print()
# Compare MSE at GCV, L-curve, and oracle lambda
mse_gcv = np.sum((tikhonov_solution(J, y_obs, lambda_gcv) - theta_star)**2)
mse_lcurve = np.sum((tikhonov_solution(J, y_obs, lambda_corner) - theta_star)**2)
mse_oracle = min(oracle_mses)
mse_unreg = np.sum((tikhonov_solution(J, y_obs, 1e-6) - theta_star)**2)

print(f"Resulting MSE comparison:")
print(f"  Unregularised (λ→0): {mse_unreg:.4e}")
print(f"  L-curve corner λ:    {mse_lcurve:.4e}")
print(f"  GCV-selected λ:      {mse_gcv:.4e}")
print(f"  Oracle optimal λ:    {mse_oracle:.4e}")
print()
print("→ Both L-curve and GCV achieve MSE close to oracle without knowing δ.")

§5 — Validation Summary

print("=" * 60)
print("VALIDATION SUMMARY — Tikhonov Regularisation")
print("=" * 60)
print()

# Check 1: Filter factor = 0.5 at crossover
ff_cross = filter_factor(2.0, 2.0)
check1 = abs(ff_cross - 0.5) < 1e-12
print(f"✓ Filter factor f(σ,λ=σ) = 0.5: computed={ff_cross:.6f}  [{check1}]")

# Check 2: Filter factor → 1 for σ >> λ
ff_large = filter_factor(100.0, 0.01)
check2 = abs(ff_large - 1.0) < 0.0001
print(f"✓ f(100,0.01) → 1:  computed={ff_large:.6f}  [{check2}]")

# Check 3: Filter factor → 0 for σ << λ
ff_small = filter_factor(0.01, 100.0)
check3 = ff_small < 1e-7
print(f"✓ f(0.01,100) → 0:  computed={ff_small:.2e}  [{check3}]")

# Check 4: Bias^2 + Variance = theory MSE (at one lambda)
lam_check = 0.1
_, s_c, Vt_c = np.linalg.svd(J, full_matrices=False)
tv = Vt_c @ theta_star
theory_bias = np.sum((lam_check**2 / (s_c**2 + lam_check**2))**2 * tv**2)
theory_var = DELTA**2 * np.sum(s_c**2 / (s_c**2 + lam_check**2)**2)
theory_total = theory_bias + theory_var
mc_est = []
for _ in range(3000):
    th_est = tikhonov_solution(J, J @ theta_star + RNG.standard_normal(N_OBS) * DELTA, lam_check)
    mc_est.append(np.sum((th_est - theta_star)**2))
mc_total = np.mean(mc_est)
check4 = abs(mc_total - theory_total) / theory_total < 0.05
print(f"✓ Bias²+Var: MC={mc_total:.4e}, theory={theory_total:.4e}  [{check4}]")

# Check 5: L-curve corner within 1 order of magnitude of oracle
check5 = abs(np.log10(lambda_corner) - np.log10(lambda_oracle)) < 1.5
print(f"✓ L-curve corner near oracle: log₁₀(λ_corner/λ_oracle)={np.log10(lambda_corner/lambda_oracle):.2f}  [{check5}]")

# Check 6: GCV within 1 order of magnitude of oracle
check6 = abs(np.log10(lambda_gcv) - np.log10(lambda_oracle)) < 1.5
print(f"✓ GCV near oracle: log₁₀(λ_GCV/λ_oracle)={np.log10(lambda_gcv/lambda_oracle):.2f}  [{check6}]")

print()
all_pass = check1 and check2 and check3 and check4 and check5 and check6
print(f"All checks: {'PASS' if all_pass else 'FAIL'}")