1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED=7
  6rng=np.random.default_rng(SEED)
  7phi=(1+math.sqrt(5))/2
  8
  9# Composite convex objective: h is logistic loss, f is elastic-net, g is box indicator on Kx=x.
 10def make_problem(n=320,d=12):
 11    X=rng.normal(size=(n,d)); X[:,0]*=12.0; X[:,1]*=4.0  # sharp and flatter directions
 12    w=rng.normal(size=d); y=(X@w+rng.normal(scale=1.0,size=n)>0).astype(float)
 13    return X,y
 14
 15def sigmoid(t): return 1/(1+np.exp(-np.clip(t,-40,40)))
 16def loss_grad(x,X,y,alpha,rho):
 17    p=sigmoid(X@x)
 18    h=np.mean(np.logaddexp(0,X@x)-y*(X@x))
 19    # h excludes f and g; derivative is minibatch/full gradient
 20    return h, X.T@(p-y)/len(y)
 21
 22def prox_f(v,lam,alpha=2e-3,rho=2e-3):
 23    return np.sign(v)*np.maximum(np.abs(v)-lam*alpha,0)/(1+lam*rho)
 24
 25def gstar_prox(y,sigma,c):
 26    # prox_{sigma g*}(y)=y-sigma clip(y/sigma,-c,c), K=I
 27    return y-sigma*np.clip(y/sigma,-c,c)
 28
 29def full_obj(x,X,y,alpha=2e-3,rho=2e-3,c=.8):
 30    h,_=loss_grad(x,X,y,alpha,rho)
 31    return h+alpha*np.abs(x).sum()+rho*.5*(x@x), max(0.,np.abs(x).max()-c)
 32
 33def grpd(X,y,steps,lam0,eta=.95,sigma=.5,c=.8,alpha=2e-3,rho=2e-3):
 34    d=X.shape[1]; x=np.zeros(d); z=x.copy(); dual=np.zeros(d); oldx=x.copy(); oldq=np.zeros(d); Lhat=1.0
 35    rec=[]
 36    for k in range(steps):
 37        h,q=loss_grad(x,X,y,alpha,rho)
 38        if k>0:
 39            Lhat=np.linalg.norm(q-oldq)/(np.linalg.norm(x-oldx)+1e-8)
 40            # robust EMA as explicitly allowed by idea
 41            Lhat=.8*Lhat+.2*Lhat_prev
 42        Lhat_prev=Lhat
 43        lam=eta/(Lhat+sigma+1e-8)
 44        # retain the requested initial-step stress: first iteration starts from lam0
 45        if k==0: lam=lam0
 46        xn=prox_f(z-lam*(q+dual),lam,alpha,rho)
 47        zn=((phi-1)/phi)*xn+(1/phi)*z
 48        dualn=gstar_prox(dual+sigma*zn,sigma,c)
 49        oldx,oldq=x.copy(),q.copy(); x,z,dual=xn,zn,dualn
 50        obj,v=full_obj(x,X,y,alpha,rho,c)
 51        rec.append((obj,v,lam,Lhat,np.linalg.norm(q)))
 52        if not np.isfinite(obj) or np.linalg.norm(x)>1e8: break
 53    return rec,x
 54
 55def proxgrad(X,y,steps,lam,alpha=2e-3,rho=2e-3,c=.8):
 56    x=np.zeros(X.shape[1]); rec=[]
 57    for k in range(steps):
 58        _,q=loss_grad(x,X,y,alpha,rho)
 59        x=prox_f(x-lam*q,lam,alpha,rho)
 60        obj,v=full_obj(x,X,y,alpha,rho,c); rec.append((obj,v,lam,np.nan,np.linalg.norm(q)))
 61        if not np.isfinite(obj) or np.linalg.norm(x)>1e8: break
 62    return rec,x
 63
 64def projected_pg(X,y,steps,lam,alpha=2e-3,rho=2e-3,c=.8):
 65    x=np.zeros(X.shape[1]); rec=[]
 66    for k in range(steps):
 67        _,q=loss_grad(x,X,y,alpha,rho)
 68        x=np.clip(prox_f(x-lam*q,lam,alpha,rho),-c,c)
 69        obj,v=full_obj(x,X,y,alpha,rho,c); rec.append((obj,v,lam,np.nan,np.linalg.norm(q)))
 70        if not np.isfinite(obj) or np.linalg.norm(x)>1e8: break
 71    return rec,x
 72
 73def adamw(X,y,steps,lr,alpha=2e-3,rho=2e-3,c=.8):
 74    x=np.zeros(X.shape[1]); m=np.zeros_like(x); v=np.zeros_like(x); rec=[]
 75    for k in range(1,steps+1):
 76        _,q=loss_grad(x,X,y,alpha,rho); m=.9*m+.1*q; v=.999*v+.001*q*q
 77        x=x-lr*(m/(1-.9**k))/(np.sqrt(v/(1-.999**k))+1e-8)-lr*rho*x
 78        obj,viol=full_obj(x,X,y,alpha,rho,c); rec.append((obj,viol,lr,np.nan,np.linalg.norm(q)))
 79        if not np.isfinite(obj) or np.linalg.norm(x)>1e8: break
 80    return rec,x
 81
 82def run():
 83    X,y=make_problem(); steps=150
 84    # Initial step sweep is the falsifiable stability test. PG has no constraint dual.
 85    # Use identical full-gradient evaluations; constants chosen around sharp curvature.
 86    spectral=np.linalg.eigvalsh(X.T@X/len(y)).max()/4
 87    candidates=[.02,.05,.1,.2,.5,1.,2.]
 88    rows=[]
 89    for lr in candidates:
 90        for name,fn in [('PG',lambda lr:proxgrad(X,y,steps,lr)),('ProjectedPG',lambda lr:projected_pg(X,y,steps,lr)),('GRPD',lambda lr:grpd(X,y,steps,lr)) ,('AdamW',lambda lr:adamw(X,y,steps,lr))]:
 91            rec,_=fn(lr); finite=bool(rec) and np.isfinite(rec[-1][0]) and rec[-1][0]<1e6
 92            rows.append({'method':name,'lr':lr,'final_loss':float(rec[-1][0]) if rec else None,'final_violation':float(rec[-1][1]) if rec else None,'finite':bool(finite),'min_loss':float(min(r[0] for r in rec)) if rec else None,'steps':len(rec)})
 93    # Representative trajectories at a stress level where behavior differs, plus math checks.
 94    checks={'phi':phi,'extrapolation_weights':[(phi-1)/phi,1/phi], 'weights_sum':((phi-1)/phi+1/phi)}
 95    # finite-difference curvature estimate agrees with local Hessian directional curvature
 96    x=rng.normal(size=X.shape[1])*0.1; dx=rng.normal(size=x.shape); dx/=np.linalg.norm(dx)
 97    _,q1=loss_grad(x,X,y,0,0); _,q2=loss_grad(x+1e-5*dx,X,y,0,0)
 98    est=np.linalg.norm(q2-q1)/(1e-5+1e-8)
 99    p=sigmoid(X@x); Hdx=X.T@((p*(1-p))*(X@dx))/len(y)
100    true=float(np.linalg.norm(Hdx))
101    checks.update({'curvature_estimate':float(est),'analytic_hessian_action':true,'curvature_relative_error':float(abs(est-true)/(true+1e-12))})
102    out={'seed':SEED,'spectral_quadratic_curvature_bound':float(spectral),'checks':checks,'results':rows}
103    Path('results.json').write_text(json.dumps(out,indent=2))
104    print(json.dumps(out,indent=2))
105if __name__=='__main__': run()