1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5from scipy.optimize import brentq
  6
  7
  8def g_filter(s, nu, t):
  9    d = float(s - nu)
 10    if d == 0.0:
 11        return float(t)
 12    return -math.expm1(-d * t) / d
 13
 14
 15def f_filter(s, nu, t):
 16    return float(s) * g_filter(s, nu, t)
 17
 18
 19def discrete_g(s, nu, eta, K):
 20    """Coefficient filter of K exact Euler updates from beta_0=0."""
 21    q = 1.0 - eta * (s - nu)
 22    if abs(q - 1.0) < 1e-12:
 23        return eta * K
 24    return eta * (1.0 - q ** K) / (1.0 - q)
 25
 26
 27def root_filter(nu, t, lo=1e-8, hi=5.0):
 28    grid = np.linspace(lo, hi, 4001)
 29    vals = [f_filter(s, nu, t) - 1.0 for s in grid]
 30    for a, b, fa, fb in zip(grid[:-1], grid[1:], vals[:-1], vals[1:]):
 31        if fa * fb < 0:
 32            return float(brentq(lambda s: f_filter(s, nu, t) - 1.0, a, b))
 33    return None
 34
 35
 36def core_checks():
 37    nu, t = 1.0, 1.5
 38    # Prediction 1: the pole is removable, with f(nu)=nu*t, and error is O(|s-nu|).
 39    eps = 10.0 ** (-np.arange(2, 9))
 40    errors = np.array([abs(f_filter(nu + e, nu, t) - nu * t) for e in eps])
 41    slopes = np.diff(np.log(errors)) / np.diff(np.log(eps))
 42    pole = {
 43        "predicted_limit": nu * t,
 44        "observed_at_pole": f_filter(nu, nu, t),
 45        "exact_error": abs(f_filter(nu, nu, t) - nu * t),
 46        "errors_eps_1e-2_to_1e-8": errors.tolist(),
 47        "median_loglog_slope": float(np.median(slopes[-4:])),
 48        "prediction": "limit is nu*t and off-pole error scales linearly in |s-nu|",
 49    }
 50
 51    # Prediction 2: for fixed nu,t, f(s)-1 changes sign once; compare analytic
 52    # root with a grid observation, and sweep t to verify the predicted transition.
 53    crossover_rows = []
 54    for tt in [0.5, 1.0, 1.5, 2.0, 3.0]:
 55        pred = root_filter(nu, tt)
 56        grid = np.linspace(1e-5, 5.0, 100001)
 57        obs = float(grid[np.argmin(np.abs(np.array([f_filter(s, nu, tt) for s in grid]) - 1))])
 58        crossover_rows.append({"t": tt, "predicted_root": pred, "observed_grid_root": obs,
 59                               "absolute_error": abs(pred - obs) if pred is not None else None,
 60                               "f_at_nu": f_filter(nu, nu, tt)})
 61    mixed = {"rows": crossover_rows,
 62             "example_low_f": f_filter(0.2, nu, t),
 63             "example_high_f": f_filter(2.0, nu, t),
 64             "mixed_sign_observed": f_filter(0.2, nu, t) < 1 < f_filter(2.0, nu, t),
 65             "prediction": "crossover moves with t and f(nu)=nu*t crosses one at t=1/nu"}
 66
 67    # Prediction 3: for a positive-curvature mode s>nu, Euler stability ends at
 68    # eta*=2/(s-nu). Verify by sweeping eta and checking the exact multiplier.
 69    s = 2.0
 70    predicted_eta = 2.0 / (s - nu)
 71    etas = np.linspace(0.01, 3.0, 10000)
 72    factors = np.abs(1.0 - etas * (s - nu))
 73    stable = etas[factors <= 1.0 + 1e-12]
 74    observed_eta = float(stable.max())
 75    stability = {"s": s, "nu": nu, "predicted_boundary": predicted_eta,
 76                 "observed_sweep_boundary": observed_eta,
 77                 "absolute_error": abs(predicted_eta - observed_eta),
 78                 "factor_below_boundary": float(1 - (observed_eta - .001) * (s - nu)),
 79                 "factor_above_boundary": float(1 - (observed_eta + .001) * (s - nu)),
 80                 "prediction": "|1-eta*(s-nu)|<=1; instability beyond eta=2/(s-nu)"}
 81
 82    # Small-step Euler-vs-flow check: discrete filter converges to continuous filter.
 83    flow_rows = []
 84    for eta in [0.1, 0.05, 0.02, 0.01]:
 85        K = int(round(t / eta))
 86        ss = 0.7
 87        flow_rows.append({"eta": eta, "K": K,
 88                          "flow_g": g_filter(ss, nu, t),
 89                          "euler_g": discrete_g(ss, nu, eta, K),
 90                          "abs_error": abs(discrete_g(ss, nu, eta, K) - g_filter(ss, nu, t))})
 91    return {"removable_pole": pole, "mixed_sign_crossover": mixed,
 92            "euler_stability_boundary": stability, "euler_to_flow": flow_rows}
 93
 94
 95def regression_experiment():
 96    rng = np.random.default_rng(1285)
 97    d, n_train, n_test = 30, 120, 3000
 98    Q, _ = np.linalg.qr(rng.normal(size=(d, d)))
 99    eig = np.geomspace(0.05, 3.0, d)
100    A = Q @ np.diag(np.sqrt(eig))
101    X = rng.normal(size=(n_train, d)) @ A.T
102    Xt = rng.normal(size=(n_test, d)) @ A.T
103    beta_eig = np.zeros(d); beta_eig[:5] = [2.0, 1.5, 1.0, .8, .6]
104    beta_true = Q @ beta_eig
105    y = X @ beta_true + rng.normal(scale=.35, size=n_train)
106    yt = Xt @ beta_true + rng.normal(scale=.35, size=n_test)
107    S, b = X.T @ X / n_train, X.T @ y / n_train
108    K, eta, nu = 30, .05, 1.0
109    def run(kind):
110        beta = np.zeros(d)
111        for _ in range(K):
112            grad = S @ beta - b
113            if kind == 'gd': beta -= eta * grad
114            elif kind == 'ns-gd': beta = (1 + eta * nu) * beta - eta * grad
115            else: beta -= eta * (grad + .1 * beta)
116        return beta
117    out = {}
118    for kind in ('gd', 'positive-ridge', 'ns-gd'):
119        beta = run(kind)
120        out[kind] = {'test_mse': float(np.mean((Xt @ beta - yt)**2)),
121                     'train_mse': float(np.mean((X @ beta - y)**2)),
122                     'parameter_norm': float(np.linalg.norm(beta))}
123    out['settings'] = {'K': K, 'eta': eta, 'nu': nu, 'n_train': n_train, 'dimension': d}
124    return out
125
126
127def main():
128    out = {'core_checks': core_checks(), 'regression': regression_experiment()}
129    Path('results.json').write_text(json.dumps(out, indent=2))
130    print(json.dumps(out, indent=2))
131
132if __name__ == '__main__': main()