1import json
  2import math
  3import numpy as np
  4
  5# Continuation Maps toy MVP.
  6# Quadratic objective: L(w)=1/2 w^T A w, with A=diag(1,4).
  7# Decoupled weight decay changes the update multiplier to
  8# q_i = 1 - eta*(eigenvalue_i + lambda).  Thus stability requires
  9# max_i |q_i| < 1, and the first boundary is eta=2/(4+lambda).
 10
 11SEED = 2208
 12A = np.array([1.0, 4.0])
 13W0 = np.array([1.0, 1.0])
 14STEPS = 160
 15WINDOW = 30
 16
 17def run(eta, lam, steps=STEPS):
 18    w = W0.copy()
 19    norms = [np.linalg.norm(w)]
 20    losses = [0.5 * np.sum(A * w * w)]
 21    for _ in range(steps):
 22        g = A * w
 23        w = w - eta * (g + lam * w)
 24        norms.append(np.linalg.norm(w))
 25        losses.append(0.5 * np.sum(A * w * w))
 26    return np.asarray(losses), np.asarray(norms)
 27
 28def feature(eta, lam):
 29    # Late-time normalized temporal growth, analogous to a drift/oscillation
 30    # feature. For this linear system F=0 exactly at unit spectral radius.
 31    _, n = run(eta, lam)
 32    r = n[-WINDOW:] / np.maximum(n[-WINDOW-1:-1], 1e-300)
 33    return float(np.mean(np.log(np.maximum(r, 1e-300))))
 34
 35def exact_boundary(lam):
 36    return 2.0 / (4.0 + lam)
 37
 38def bisect_eta(lam, lo, hi, tol=2e-7, maxit=60):
 39    flo, fhi = feature(lo, lam), feature(hi, lam)
 40    # Expand only if a caller supplied a non-bracketing interval.
 41    for _ in range(20):
 42        if flo <= 0 <= fhi: break
 43        if flo > 0:
 44            lo *= 0.5; flo = feature(lo, lam)
 45        if fhi < 0:
 46            hi *= 1.25; fhi = feature(hi, lam)
 47    for _ in range(maxit):
 48        mid = (lo + hi) / 2
 49        fm = feature(mid, lam)
 50        if abs(fm) < 1e-9 or hi - lo < tol: return mid
 51        if fm > 0: hi = mid
 52        else: lo = mid
 53    return (lo + hi) / 2
 54
 55def main():
 56    np.random.seed(SEED)
 57    # Prediction 1: divergence starts at eta*(4+lambda)=2.
 58    lam_sweep = np.linspace(0, 4, 9)
 59    boundary_rows = []
 60    for lam in lam_sweep:
 61        pred = exact_boundary(lam)
 62        # Direct sign sweep, with a fine grid, independently of continuation.
 63        etas = np.linspace(0.05, 0.55, 1001)
 64        fs = np.array([feature(e, lam) for e in etas])
 65        ix = np.flatnonzero(fs >= 0)
 66        measured = float(etas[ix[0]]) if len(ix) else float('nan')
 67        boundary_rows.append({"lambda": float(lam), "predicted_eta": pred,
 68                              "measured_eta": measured,
 69                              "relative_error": abs(measured-pred)/pred})
 70
 71    # Prediction 2: eta_boundary*(4+lambda) is constant (=2).
 72    scaled = [r["measured_eta"] * (4 + r["lambda"]) for r in boundary_rows]
 73    scaling_error = abs(np.mean(scaled) - 2.0) / 2.0
 74
 75    # Prediction 3: secant predictor + one-dimensional correction follows the
 76    # same curve. Start with two measured boundary points and advance lambda.
 77    cont_lams = np.linspace(0, 4, 9)
 78    p0 = bisect_eta(cont_lams[0], 0.2, 0.7)
 79    p1 = bisect_eta(cont_lams[1], 0.2, 0.7)
 80    points = [p0, p1]
 81    for lam in cont_lams[2:]:
 82        pred = points[-1] + (points[-1] - points[-2])
 83        # A bracket around secant prediction; the correction is a scalar sweep.
 84        lo, hi = max(0.03, pred * 0.65), pred * 1.35
 85        corrected = bisect_eta(float(lam), lo, hi)
 86        points.append(corrected)
 87    continuation_err = [abs(e-exact_boundary(l))/exact_boundary(l)
 88                        for e,l in zip(points, cont_lams)]
 89
 90    # Local normal-direction sign switch at a regular crossing.
 91    lam0 = 2.0
 92    e0 = exact_boundary(lam0)
 93    eps = 0.01 * e0
 94    local = {"eta_below": e0-eps, "F_below": feature(e0-eps, lam0),
 95             "eta_above": e0+eps, "F_above": feature(e0+eps, lam0),
 96             "sign_switch": feature(e0-eps,lam0) < 0 < feature(e0+eps,lam0)}
 97
 98    # Standard baseline: evaluate a uniform 20x20 parameter grid and infer
 99    # each lambda-column boundary from the first unstable eta.
100    grid_etas = np.linspace(0.05, 0.55, 20)
101    grid_lams = np.linspace(0.0, 4.0, 20)
102    grid_F = np.array([[feature(e, l) for e in grid_etas] for l in grid_lams])
103    grid_boundary = []
104    for j, l in enumerate(grid_lams):
105        unstable = np.flatnonzero(grid_F[j] >= 0)
106        grid_boundary.append(float(grid_etas[unstable[0]]) if len(unstable) else float('nan'))
107    grid_boundary = np.asarray(grid_boundary)
108    grid_valid = np.isfinite(grid_boundary)
109    exact_grid = np.array([exact_boundary(l) for l in grid_lams])
110    grid_rel_error = (np.abs(grid_boundary[grid_valid] - exact_grid[grid_valid]) /
111                      exact_grid[grid_valid])
112
113    result = {
114      "seed": SEED,
115      "predictions": {
116        "stability": "divergence when eta*(4+lambda)>2",
117        "scaling": "eta_boundary*(4+lambda)=2",
118        "continuation": "secant plus scalar correction tracks boundary"
119      },
120      "direct_sweep": boundary_rows,
121      "mean_relative_boundary_error": float(np.mean([r["relative_error"] for r in boundary_rows])),
122      "scaling_mean": float(np.mean(scaled)),
123      "scaling_relative_error": float(scaling_error),
124      "continuation": {"lambdas": cont_lams.tolist(), "etas": points,
125                        "relative_errors": continuation_err,
126                        "mean_relative_error": float(np.mean(continuation_err)),
127                        "runs": 2 + 7 * 2, "dense_grid_20x20_runs": 400,
128                        "grid_mean_relative_error": float(np.mean(grid_rel_error)),
129                        "grid_boundary_points": int(np.sum(grid_valid))},
130      "local_normal_check": local
131    }
132    with open("results.json", "w") as f: json.dump(result, f, indent=2)
133    print(json.dumps(result, indent=2))
134
135if __name__ == "__main__": main()