1import json
  2import numpy as np
  3
  4SEED = 1053
  5rng = np.random.default_rng(SEED)
  6
  7# Two positive-definite phase Hessians. They do not commute.
  8H1 = np.array([[4.0, 1.0], [1.0, 1.0]])
  9H2 = np.array([[1.5, -0.7], [-0.7, 3.0]])
 10I = np.eye(2)
 11
 12
 13def monodromy(scale, order=(1, 2)):
 14    # eta1=1.6s and eta2=.4s: same two-step average LR s as baseline.
 15    e1, e2 = 1.6 * scale, 0.4 * scale
 16    J1, J2 = I - e1 * H1, I - e2 * H2
 17    return (J2 @ J1) if order == (1, 2) else (J1 @ J2)
 18
 19
 20def rho(M):
 21    return float(np.max(np.abs(np.linalg.eigvals(M))))
 22
 23
 24def boundary():
 25    # Largest stable scale, found from the exact rho(M)=1 crossing.
 26    lo, hi = 0.0, 2.0
 27    while rho(monodromy(hi)) < 1.0:
 28        hi *= 2
 29    for _ in range(70):
 30        mid = (lo + hi) / 2
 31        if rho(monodromy(mid)) < 1:
 32            lo = mid
 33        else:
 34            hi = mid
 35    return (lo + hi) / 2
 36
 37
 38def growth_check(s, periods=80):
 39    M = monodromy(s)
 40    vals, vecs = np.linalg.eig(M)
 41    j = int(np.argmax(np.abs(vals)))
 42    v = np.real(vecs[:, j])
 43    v /= np.linalg.norm(v)
 44    x = v.copy()
 45    norms = []
 46    for _ in range(periods):
 47        x = M @ x
 48        norms.append(np.linalg.norm(x))
 49    # Geometric mean removes initial normalization and reports empirical Floquet growth.
 50    empirical = float((norms[-1] / norms[0]) ** (1.0 / (periods - 1)))
 51    return rho(M), empirical, float(norms[-1])
 52
 53
 54def commutator_scaling():
 55    C = H2 @ H1 - H1 @ H2
 56    c_norm = np.linalg.norm(C, 2)
 57    rows = []
 58    for s in [0.01, 0.02, 0.04, 0.08, 0.16]:
 59        # Exact antisymmetric part is eta1*eta2*(H2 H1-H1 H2).
 60        M = monodromy(s)
 61        observed = np.linalg.norm(M - M.T, 2)
 62        predicted = (1.6 * s) * (0.4 * s) * c_norm
 63        rows.append({"scale": s, "observed": observed, "predicted": predicted,
 64                     "ratio": observed / predicted})
 65    return {"commutator_norm": c_norm, "rows": rows}
 66
 67
 68def quadratic_optimizer(scale=0.20, periods=60):
 69    # Both methods see the same alternating phase gradients. Baseline uses constant
 70    # LR s; Floquet uses the two phase LRs while preserving the same average LR.
 71    x0 = np.array([2.0, -1.5])
 72    x_const, x_periodic = x0.copy(), x0.copy()
 73    e1, e2 = 1.6 * scale, .4 * scale
 74    trace = []
 75    Hbar = .5 * (H1 + H2)
 76    def objective(x):
 77        return float(.5 * x @ Hbar @ x)
 78    for p in range(periods):
 79        # phase 1
 80        x_const -= scale * (H1 @ x_const)
 81        x_periodic -= e1 * (H1 @ x_periodic)
 82        # phase 2
 83        x_const -= scale * (H2 @ x_const)
 84        x_periodic -= e2 * (H2 @ x_periodic)
 85        if p in (0, 4, 9, 19, 39, periods - 1):
 86            trace.append({"period": p + 1, "constant_loss": objective(x_const),
 87                          "periodic_loss": objective(x_periodic)})
 88    return {"scale": scale, "periods": periods, "trace": trace,
 89            "final_constant_loss": objective(x_const),
 90            "final_periodic_loss": objective(x_periodic)}
 91
 92
 93def main():
 94    s_star = boundary()
 95    # Boundary prediction is the formula rho[(I-e2 H2)(I-e1 H1)]=1;
 96    # observed transition is independently classified by long-run norm behavior.
 97    scales = [0.90 * s_star, 0.99 * s_star, 1.01 * s_star, 1.10 * s_star]
 98    growth = []
 99    for s in scales:
100        r, empirical, final_norm = growth_check(s)
101        growth.append({"scale": s, "rho_predicted": r,
102                       "empirical_per_period_growth": empirical,
103                       "final_perturbation_norm": final_norm,
104                       "classification": "stable" if empirical < 1 else "unstable"})
105    # A direct sweep gives an independently observable transition bracket.
106    sweep = []
107    for s in np.linspace(.8 * s_star, 1.2 * s_star, 17):
108        r, empirical, _ = growth_check(float(s), periods=50)
109        sweep.append([float(s), r, empirical])
110    order_difference = np.linalg.norm(monodromy(.25, (1, 2)) - monodromy(.25, (2, 1)))
111    result = {
112        "seed": SEED,
113        "hessians": {"H1": H1.tolist(), "H2": H2.tolist()},
114        "predictions": {
115            "critical_scale_exact": s_star,
116            "stability_rule": "rho(M)<1 predicts decay and rho(M)>1 predicts growth",
117            "commutator_rule": "||M-M.T|| = eta1*eta2*||[H2,H1]|| for symmetric H phases"
118        },
119        "growth_near_boundary": growth,
120        "boundary_sweep": sweep,
121        "commutator_scaling": commutator_scaling(),
122        "reversed_order_matrix_difference_at_scale_.25": order_difference,
123        "quadratic_optimizer": quadratic_optimizer()
124    }
125    print(json.dumps(result, indent=2))
126
127if __name__ == "__main__":
128    main()