1import json, math
 2import numpy as np
 3
 4SEED = 2807
 5
 6def boundary(lam, r, sig, B):
 7    return 2*lam/(r*(lam*lam + sig*sig/B))
 8
 9def toy():
10    rng=np.random.default_rng(SEED); lam,r,sig=1.,.8,1.5
11    contraction=[]
12    for alpha in (.2,.5,1.0,1.4):
13        a=1-alpha*r*lam; x=1.; ratios=[]
14        for _ in range(30):
15            y=a*x; ratios.append(abs(y/x)); x=y
16        contraction.append({'alpha':alpha,'predicted':abs(a),'observed':float(np.mean(ratios))})
17    # Directly estimate E[(x_{k+1}/x_k)^2], avoiding rare-event trajectory bias.
18    rows=[]; Bs=[1,2,4,8,16,32]
19    for B in Bs:
20        pred=boundary(lam,r,sig,B)
21        grid=np.linspace(pred*.72,pred*1.28,41)
22        vals=[]
23        for alpha in grid:
24            xi=rng.normal(0,sig/math.sqrt(B),size=500000)
25            factor=1-alpha*r*(lam+xi)
26            m2=float(np.mean(factor*factor))
27            vals.append((alpha,m2-1))
28        cross=min(vals,key=lambda z:abs(z[1]))[0]
29        for (a,s),(b,t) in zip(vals[:-1],vals[1:]):
30            if s*t<=0:
31                cross=a+(0-s)*(b-a)/(t-s); break
32        rows.append({'B':B,'predicted_alpha_boundary':pred,'observed_alpha_boundary':float(cross),'relative_error':float(abs(cross-pred)/pred)})
33    # Additive response-noise floor: stationary E[x^2] = alpha^2*sigma^2/B/(1-a^2).
34    alpha=.25; add_sig=1.2; scale=[]; a=1-alpha*r*lam
35    predicted_per_B=(alpha*alpha*add_sig*add_sig)/(1-a*a)
36    for B in (1,4,16):
37        x=np.zeros(50000); vals=[]
38        for k in range(2200):
39            x=a*x + alpha*rng.normal(0,add_sig/math.sqrt(B),size=x.size)
40            if k>=1000: vals.append(np.mean(x*x))
41        obs=float(np.mean(vals))
42        scale.append({'B':B,'predicted_second_moment':predicted_per_B/B,'observed_second_moment':obs,'B_times_moment':B*obs})
43    return {'contraction':contraction,'boundary':rows,'variance_scaling':scale}
44
45def optimizer_demo():
46    rng=np.random.default_rng(SEED+1); n=1600; d=12
47    X=rng.normal(size=(n,d)); true=rng.normal(size=d); y=X@true+rng.normal(0,.25,n)
48    def run(kind, seed):
49        rr=np.random.default_rng(seed); w=np.zeros(d); m=np.zeros(d); v=np.zeros(d); losses=[]; alpha=.35; B=32
50        for k in range(300):
51            ix=rr.integers(0,n,B); xb=X[ix]; yb=y[ix]; g=xb.T@(xb@w-yb)/B
52            if kind=='adam':
53                m=.9*m+.1*g; v=.999*v+.001*g*g; step=.06*m/(np.sqrt(v)+1e-8)
54            elif kind=='prox': step=.35*g
55            else: step=alpha*.35*g
56            w-=step; losses.append(float(np.mean((X@w-y)**2)))
57        return losses
58    out={k:run(k,SEED+2) for k in ('adam','prox','relaxed_prox')}
59    return {k:{'final_loss':v[-1],'best_loss':min(v)} for k,v in out.items()}
60
61if __name__=='__main__':
62    result={'toy':toy(),'optimizer_demo':optimizer_demo()}
63    with open('results.json','w') as f: json.dump(result,f,indent=2)
64    print(json.dumps(result,indent=2))