1import json
  2import numpy as np
  3
  4
  5def simulate(lambdas, eta, kappa, x0, y0, steps):
  6    H = np.diag(np.asarray(lambdas, dtype=float))
  7    x, y = np.array(x0, dtype=float).copy(), np.array(y0, dtype=float).copy()
  8    xs, ys, losses = [x.copy()], [y.copy()], []
  9    for _ in range(steps):
 10        # Exact discrete Euler version of the proposed coupled optimizer.
 11        # Simultaneous Euler update: both right-hand sides use the old state.
 12        xo, yo = x.copy(), y.copy()
 13        gx, gy = H @ xo, H @ yo
 14        x = xo - eta * (gx + kappa * (xo - yo))
 15        y = yo - eta * (gy + kappa * (yo - xo))
 16        xs.append(x.copy()); ys.append(y.copy())
 17        losses.append(0.5 * (x @ H @ x + y @ H @ y))
 18    return np.asarray(xs), np.asarray(ys), np.asarray(losses)
 19
 20
 21def modal_rate(series, eta, start=0, end=None):
 22    z = np.abs(np.asarray(series)[start:end])
 23    z = z[z > 1e-14]
 24    if len(z) < 2:
 25        return float("nan")
 26    # Geometric per-step rate, converted to continuous-time equivalent.
 27    return -np.polyfit(np.arange(len(z)), np.log(z), 1)[0] / eta
 28
 29
 30def run():
 31    # Diagonal quadratic is a local Hessian model, with an intentionally separated
 32    # slow mode so the predicted transition is visible.
 33    lam = np.array([0.10, 0.70, 2.00])
 34    eta = 0.05
 35    x0 = np.array([1.0, 0.8, 0.4])
 36    results = {"lambdas": lam.tolist(), "eta": eta}
 37
 38    # Prediction 1: prepared packet cancels c_1 exactly, while ordinary packet
 39    # retains c_1=x_1. Also test a tenfold cancellation perturbation.
 40    y_unprepared = x0.copy()
 41    y_prepared = x0.copy(); y_prepared[0] = -x0[0]
 42    y_perturbed = y_prepared.copy(); y_perturbed[0] += 0.10 * x0[0]
 43    rows = []
 44    for name, y0 in [("unprepared", y_unprepared), ("prepared", y_prepared),
 45                     ("10pct_residual", y_perturbed)]:
 46        xs, ys, _ = simulate(lam, eta, 0.35, x0, y0, 30)
 47        c1 = 0.5 * (xs[:, 0] + ys[:, 0])
 48        rows.append({"case": name, "initial_abs_c1_over_abs_a1": float(abs(c1[0]) / abs(x0[0])),
 49                     "abs_c1_after_30": float(abs(c1[-1]))})
 50    results["cancellation"] = rows
 51
 52    # Prediction 2: after cancellation, rate is min(lambda_2, lambda_1+2*kappa).
 53    # Rates are measured from the exact modal trajectories and compared to theory.
 54    rate_rows = []
 55    for k in [0.0, 0.05, 0.15, 0.30, 0.60]:
 56        xs, ys, _ = simulate(lam, eta, k, x0, y_prepared, 180)
 57        c2 = 0.5 * (xs[:, 1] + ys[:, 1])
 58        d1 = 0.5 * (xs[:, 0] - ys[:, 0])
 59        # The slowest surviving term is whichever of c2 and d1 decays slower.
 60        r2 = -np.log(abs(1 - eta * lam[1])) / eta
 61        r1d = -np.log(abs(1 - eta * (lam[0] + 2*k))) / eta
 62        predicted = min(r2, r1d)
 63        def fit_rate(z):
 64            zz = np.abs(z[:80])
 65            idx = zz > 1e-12
 66            return float(-np.polyfit(np.arange(len(zz))[idx], np.log(zz[idx]), 1)[0] / eta)
 67        observed = min(fit_rate(c2), fit_rate(d1))
 68        rate_rows.append({"kappa": k, "predicted_rate": float(predicted),
 69                          "observed_modal_rate": float(observed),
 70                          "lambda1_plus_2k": float(lam[0] + 2*k),
 71                          "lambda2": float(lam[1])})
 72    results["rate_sweep"] = rate_rows
 73
 74    # Prediction 3: Euler stability boundary eta*(lambda_max+2*kappa)=2.
 75    # Sweep around predicted kappa critical value.
 76    kcrit = (2.0 / eta - lam.max()) / 2.0
 77    stability_rows = []
 78    for k in [kcrit - 0.5, kcrit - 0.05, kcrit, kcrit + 0.05, kcrit + 0.5]:
 79        xs, ys, loss = simulate(lam, eta, k, x0, -x0, 200)
 80        rho = max(abs(1 - eta * lam.max()), abs(1 - eta * (lam.max() + 2*k)))
 81        stability_rows.append({"kappa": float(k), "predicted_spectral_radius": float(rho),
 82                               "predicted_stable": bool(rho <= 1 + 1e-12),
 83                               "observed_unstable_200_steps": bool((not np.all(np.isfinite(loss))) or loss[-1] > loss[0] * 1.01),
 84                               "loss_final": float(loss[-1])})
 85    results["stability_sweep"] = {"predicted_kappa_critical": float(kcrit), "rows": stability_rows}
 86
 87    # Baseline comparison at equal coupled update steps: prepared packet has a
 88    # substantially lower late loss once the slow mode is removed.
 89    xb, yb, lb = simulate(lam, eta, 0.35, x0, y_unprepared, 100)
 90    xp, yp, lp = simulate(lam, eta, 0.35, x0, y_prepared, 100)
 91    results["baseline_vs_idea"] = {"unprepared_final_loss": float(lb[-1]),
 92                                    "prepared_final_loss": float(lp[-1]),
 93                                    "loss_ratio_idea_over_baseline": float(lp[-1] / lb[-1])}
 94    with open("results.json", "w") as f:
 95        json.dump(results, f, indent=2)
 96    print(json.dumps(results, indent=2))
 97
 98
 99if __name__ == "__main__":
100    run()