1import json
  2import math
  3import numpy as np
  4
  5SEED = 7
  6np.random.seed(SEED)
  7
  8# Bounded diagonal quadratic manifold. The barrier is smooth in (-5,5)^d.
  9LO, HI, LAM, GAMMA = -5.0, 5.0, 0.015, 0.03
 10H = np.linspace(0.7, 1.8, 8)
 11
 12
 13def potential(x):
 14    if np.any(x <= LO) or np.any(x >= HI):
 15        return np.inf
 16    return 0.5 * np.sum(H*x*x) - LAM*np.log(x-LO).sum() - LAM*np.log(HI-x).sum()
 17
 18
 19def grad_potential(x):
 20    return H*x + LAM*(-1.0/(x-LO) + 1.0/(HI-x))
 21
 22
 23def metric(x):
 24    return 1.0 + GAMMA*x*x
 25
 26
 27def state(x):
 28    G = metric(x)
 29    g = grad_potential(x)/G
 30    r = math.sqrt(float(np.sum(G*g*g)))
 31    return g, r
 32
 33
 34def step_fixed(x, eta=0.035, a=0.75, b=0.10, p=0.5, q=2.0, eps=1e-10):
 35    g, r0 = state(x)
 36    r = math.sqrt(r0*r0 + eps*eps)
 37    return x - eta*(a*r**(p-1.0) + b*r**(q-1.0))*g
 38
 39
 40def step_gd(x, eta=0.035):
 41    return x - eta*grad_potential(x)/metric(x)
 42
 43
 44def run(step, x0, n=1200, tol=1e-5):
 45    x = x0.copy(); vstar = potential(np.zeros_like(x)); vals=[]
 46    feasible = True; settle = None
 47    for t in range(n+1):
 48        v = potential(x)-vstar
 49        vals.append(v)
 50        feasible = feasible and np.all(np.isfinite(x)) and np.all(x > LO) and np.all(x < HI)
 51        if settle is None and v <= tol: settle=t
 52        if t < n: x=step(x)
 53    return {"settle": settle if settle is not None else n+1,
 54            "residual": float(vals[-1]), "min_residual": float(np.min(vals)),
 55            "feasible": bool(feasible)}
 56
 57
 58def math_check():
 59    # Check the exact claimed inequality using the metric gradient. Estimate
 60    # the largest valid mu on the same central neighborhood, then use 90% of it.
 61    p,q,a,b=.5,2.,.75,.10
 62    alpha,beta=(p+1)/2,(q+1)/2
 63    ratios=[]; xs=np.linspace(-3.8,3.8,2001)
 64    v0=potential(np.zeros(1))
 65    for z in xs:
 66        x=np.array([z]); V=potential(x)-v0
 67        if V > 1e-12:
 68            _,r=state(x); ratios.append(r*r/(2*V))
 69    mu=.9*min(ratios)
 70    A=a*(2*mu)**alpha; BB=b*(2*mu)**beta
 71    worst=-np.inf; violations=0
 72    for z in xs:
 73        x=np.array([z]); V=potential(x)-v0
 74        if V <= 1e-12: continue
 75        _,r=state(x)
 76        dV=-a*r**(p+1)-b*r**(q+1)
 77        rhs=-A*V**alpha-BB*V**beta
 78        worst=max(worst,dV-rhs)
 79        violations += int(dV > rhs+1e-10)
 80    bound=1/(A*(1-alpha))+1/(BB*(beta-1))
 81    return {"metric_mu_estimate":float(mu),"alpha":alpha,"beta":beta,
 82            "settling_bound":float(bound),"max_dV_minus_rhs":float(worst),
 83            "violations":violations,"samples":len(xs)-1}
 84
 85
 86def epsilon_check():
 87    # For the smoothed norm, measure residual potential at a fixed long horizon.
 88    # This is a small numerical check of the stated O(epsilon) neighborhood claim.
 89    out=[]; x0=np.full(8,3.0); vstar=potential(np.zeros(8))
 90    for eps in [1e-4, 3e-4, 1e-3, 3e-3, 1e-2]:
 91        x=x0.copy()
 92        for _ in range(3000): x=step_fixed(x,eps=eps)
 93        out.append([eps,float(potential(x)-vstar)])
 94    logs=np.log(np.array(out)); slope=float(np.polyfit(logs[:,0],logs[:,1],1)[0])
 95    return {"epsilon_residual_pairs":out,"loglog_slope":slope}
 96
 97
 98def main():
 99    results={"math_check":math_check(),"epsilon_check":epsilon_check(),"runs":[]}
100    for scale in [0.1,0.5,1.0,2.0,3.5,4.5]:
101        x0=np.full(8,scale)
102        results["runs"].append({"scale":scale,
103          "fixed_time":run(step_fixed,x0),"gradient_descent":run(step_gd,x0)})
104    results["summary"]={
105      "fixed_time_settles":[r["fixed_time"]["settle"] for r in results["runs"]],
106      "gd_settles":[r["gradient_descent"]["settle"] for r in results["runs"]],
107      "fixed_time_residuals":[r["fixed_time"]["residual"] for r in results["runs"]],
108      "gd_residuals":[r["gradient_descent"]["residual"] for r in results["runs"]]}
109    with open("results.json","w") as f: json.dump(results,f,indent=2)
110    print(json.dumps(results,indent=2))
111
112if __name__ == "__main__": main()