""" Training is numerical integration — experiments for the article "Training Is Numerical Integration: What Sixty Years of Numerical Analysis Say About Learning Rates" (numericalmethod.com, 2026). Requirements: Python >= 3.9, numpy >= 1.24. No other dependencies. Run: python3 experiments.py (writes results.json next to this file) """ import json import math import os import time import warnings import numpy as np # numpy 2.x + Apple Accelerate can emit spurious 'divide by zero/overflow in matmul' warnings. # The results were checked against exactly rounded sums (math.fsum); the warnings are silenced here. warnings.filterwarnings('ignore', message='.*encountered in matmul') # --------------------------------------------------------------------------- # A small non-convex model: one-hidden-layer tanh network on a 1-D regression task # --------------------------------------------------------------------------- class TanhNet: def __init__(self, width=16, n=64, seed=0): rng = np.random.default_rng(seed) self.x = np.linspace(-2, 2, n) self.y = np.sin(2 * self.x) + 0.3 * self.x self.width = width self.theta0 = np.concatenate([ rng.normal(0, 1.0, width), # input weights a rng.normal(0, 0.5, width), # biases b rng.normal(0, 1.0 / math.sqrt(width), width), # output weights c ]) def unpack(self, th): w = self.width return th[:w], th[w:2 * w], th[2 * w:] def loss(self, th): a, b, c = self.unpack(th) h = np.tanh(np.outer(self.x, a) + b) r = h @ c - self.y return 0.5 * float(np.mean(r ** 2)) def grad(self, th): a, b, c = self.unpack(th) z = np.outer(self.x, a) + b h = np.tanh(z) r = (h @ c - self.y) / len(self.x) gc = h.T @ r dz = np.outer(r, c) * (1 - h ** 2) ga = dz.T @ self.x gb = dz.sum(axis=0) return np.concatenate([ga, gb, gc]) def hvp(self, th, v, eps=1e-6): """Hessian-vector product by central differences of the analytic gradient.""" return (self.grad(th + eps * v) - self.grad(th - eps * v)) / (2 * eps) def rk4(f, x0, T, steps): x = x0.copy() dt = T / steps for _ in range(steps): k1 = f(x) k2 = f(x + 0.5 * dt * k1) k3 = f(x + 0.5 * dt * k2) k4 = f(x + dt * k3) x = x + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4) return x # --------------------------------------------------------------------------- # Experiment 1: backward error analysis — GD follows the modified flow to O(h^2) # gradient flow : x' = -grad L # modified flow : x' = -grad L - (h/2) H grad L ( = -grad[L + (h/4)|grad L|^2] ) # --------------------------------------------------------------------------- def experiment_1(): net = TanhNet() T = 2.0 out = [] for h in (0.2, 0.1, 0.05, 0.025, 0.0125): k = int(round(T / h)) x = net.theta0.copy() for _ in range(k): x = x - h * net.grad(x) flow = rk4(lambda th: -net.grad(th), net.theta0, T, 4000) mod = rk4(lambda th: -net.grad(th) - 0.5 * h * net.hvp(th, net.grad(th)), net.theta0, T, 4000) d_flow = float(np.linalg.norm(x - flow)) d_mod = float(np.linalg.norm(x - mod)) out.append({"h": h, "steps": k, "dist_to_gradient_flow": d_flow, "dist_to_modified_flow": d_mod}) print(f"E1 h={h:<7} |GD - flow|={d_flow:.3e} |GD - modified flow|={d_mod:.3e}") return out # --------------------------------------------------------------------------- # Experiment 2: early-stopped gradient descent is ridge regression in disguise # Synthetic return prediction: 200 weak signals, 300 training days, 5000 test days. # --------------------------------------------------------------------------- def make_signals(p=200, n_train=300, n_test=5000, seed=7): rng = np.random.default_rng(seed) # correlated signals: a few common drivers plus idiosyncratic noise F = rng.normal(size=(p, 5)) def draw(n): Z = rng.normal(size=(n, 5)) @ F.T * 0.6 + rng.normal(size=(n, p)) return Z / math.sqrt(p) beta = rng.normal(0, 1, p) * 0.35 # many weak true coefficients Xtr, Xte = draw(n_train), draw(n_test) noise = 1.0 ytr = Xtr @ beta + rng.normal(0, noise, n_train) yte = Xte @ beta + rng.normal(0, noise, n_test) return Xtr, ytr, Xte, yte, beta def r2(y, yhat): return 1 - float(np.sum((y - yhat) ** 2) / np.sum((y - y.mean()) ** 2)) def experiment_2(): Xtr, ytr, Xte, yte, beta = make_signals() n = Xtr.shape[0] A = Xtr.T @ Xtr / n b = Xtr.T @ ytr / n lam_max = float(np.linalg.eigvalsh(A).max()) out = {"lambda_max": lam_max, "gd": {}, "ridge": []} # GD paths for three step sizes, recorded against "time" t = k h for frac in (0.1, 0.5, 0.95): h = frac * 2 / lam_max w = np.zeros(A.shape[0]) path = [] t, k = 0.0, 0 checkpoints = set(np.unique(np.round(np.logspace(0, math.log10(60000), 70)).astype(int))) while k <= 60000: if k in checkpoints: path.append({"k": k, "t": k * h, "r2_test": r2(yte, Xte @ w), "r2_train": r2(ytr, Xtr @ w)}) w = w - h * (A @ w - b) k += 1 out["gd"][f"{frac:.2f}"] = {"h": h, "path": path} best = max(path, key=lambda r: r["r2_test"]) print(f"E2 GD h={h:.4f} (h*lmax={h * lam_max:.2f}): best test R2={best['r2_test']:.4f} at t=kh={best['t']:.2f}") # ridge path, parameterised by alpha; compare at alpha = 1/t evals, evecs = np.linalg.eigh(A) bt = evecs.T @ b for alpha in np.logspace(-4, 2, 80): w = evecs @ (bt / (evals + alpha)) out["ridge"].append({"alpha": float(alpha), "r2_test": r2(yte, Xte @ w), "r2_train": r2(ytr, Xtr @ w)}) best = max(out["ridge"], key=lambda r: r["r2_test"]) print(f"E2 ridge: best test R2={best['r2_test']:.4f} at alpha={best['alpha']:.4f} (1/alpha={1 / best['alpha']:.2f})") # effective degrees of freedom: df(t) = sum_i (1 - exp(-t s_i)) for gradient flow, sum_i s_i/(s_i+alpha) for ridge best_gd = max(out["gd"]["0.10"]["path"], key=lambda r: r["r2_test"]) out["eigenvalues"] = {"min": float(evals.min()), "max": float(evals.max()), "count": int(len(evals))} out["df_flow_at_best_t"] = {"t": best_gd["t"], "df": float(np.sum(1 - np.exp(-best_gd["t"] * evals)))} out["df_ridge_at_best_alpha"] = {"alpha": best["alpha"], "df": float(np.sum(evals / (evals + best["alpha"])))} print(f"E2 effective df: gradient flow at t={best_gd['t']:.1f}: {out['df_flow_at_best_t']['df']:.1f}; " f"ridge at alpha={best['alpha']:.4f}: {out['df_ridge_at_best_alpha']['df']:.1f} (of {len(evals)} parameters)") w_ols = np.linalg.lstsq(Xtr, ytr, rcond=None)[0] out["ols_r2_test"] = r2(yte, Xte @ w_ols) out["oracle_r2_test"] = r2(yte, Xte @ beta) print(f"E2 OLS test R2={out['ols_r2_test']:.4f} oracle (true beta) test R2={out['oracle_r2_test']:.4f}") return out # --------------------------------------------------------------------------- # Experiment 3: the spectral filters themselves (exact, no data needed) # GD after k steps: phi_k(s) = 1 - (1 - h s)^k ridge: s / (s + alpha) # --------------------------------------------------------------------------- def experiment_3(): out = [] t = 10.0 for frac in (0.1, 0.5, 0.95): lam_max = 1.0 h = frac * 2 / lam_max k = int(round(t / h)) ss = np.logspace(-3, 0, 200) phi_gd = 1 - (1 - h * ss) ** k phi_flow = 1 - np.exp(-t * ss) phi_ridge = ss / (ss + 1 / t) out.append({"h": h, "k": k, "s": ss.tolist(), "gd": phi_gd.tolist(), "flow": phi_flow.tolist(), "ridge": phi_ridge.tolist(), "max_gap_gd_flow": float(np.max(np.abs(phi_gd - phi_flow))), "max_gap_flow_ridge": float(np.max(np.abs(phi_flow - phi_ridge)))}) print(f"E3 h={h:.2f}: max|GD filter - flow filter|={out[-1]['max_gap_gd_flow']:.3f} max|flow - ridge|={out[-1]['max_gap_flow_ridge']:.3f}") return out if __name__ == "__main__": t0 = time.time() res = {"experiment_1": experiment_1(), "experiment_2": experiment_2(), "experiment_3": experiment_3(), "numpy": np.__version__} path = os.path.join(os.path.dirname(os.path.abspath(__file__)), "results.json") with open(path, "w") as f: json.dump(res, f, indent=1) print(f"done in {time.time() - t0:.0f}s -> {path}")