Weakly Normally Hyperbolic Cyclic Optimizer / cyclic_optimizer_experiment.py
Failed on benchmark
1import json, math, random
2from pathlib import Path
3import numpy as np
4
5SEED = 1645
6np.random.seed(SEED)
7random.seed(SEED)
8
9# Weakly normally hyperbolic cyclic optimizer on a quadratic.
10# State x=(theta,m), phase advances by one fixed step per update.
11def step(x, phi, eta0, amp, mu0, mu_amp, curvature=1.0, forcing=0.0):
12 theta, mom = float(x[0]), float(x[1])
13 eta = eta0 * (1.0 + amp * math.sin(phi))
14 mu = mu0 + mu_amp * math.cos(phi)
15 grad = curvature * theta + forcing
16 mom2 = mu * mom + grad
17 return np.array([theta - eta * mom2, mom2], dtype=float)
18
19def jacobian(phi, eta0, amp, mu0, mu_amp, curvature=1.0):
20 eta = eta0 * (1.0 + amp * math.sin(phi))
21 mu = mu0 + mu_amp * math.cos(phi)
22 # derivative of [theta-eta*(mu*m+curvature*theta), mu*m+curvature*theta]
23 return np.array([[1.0 - eta * curvature, -eta * mu],
24 [curvature, mu]], dtype=float)
25
26def cycle_map(x, pars, return_path=False):
27 P = pars['P']; dphi = 2.0 * math.pi / P
28 path = []
29 y = np.asarray(x, dtype=float).copy()
30 for k in range(P):
31 if return_path: path.append(y.copy())
32 y = step(y, k*dphi, **{q: pars[q] for q in ('eta0','amp','mu0','mu_amp','curvature','forcing')})
33 return (y, path) if return_path else y
34
35def monodromy(pars):
36 P = pars['P']; dphi = 2.0*math.pi/P
37 M = np.eye(2)
38 # perturbations propagate J_k, hence M=J_{P-1}...J_0
39 for k in range(P):
40 M = jacobian(k*dphi, pars['eta0'], pars['amp'], pars['mu0'], pars['mu_amp'], pars['curvature']) @ M
41 return M
42
43def periodic_orbit(pars):
44 # affine map is solved robustly by fixed-point iteration; fallback linear solve
45 x = np.zeros(2)
46 for _ in range(10000):
47 y = cycle_map(x, pars)
48 if np.linalg.norm(y-x) < 1e-13: return x
49 x = y
50 # Numerical finite difference around zero gives affine offset.
51 c = cycle_map(np.zeros(2), pars)
52 return np.linalg.solve(np.eye(2)-monodromy(pars), c)
53
54def empirical_decay(pars, perturb=1e-6, cycles=12):
55 xstar = periodic_orbit(pars)
56 y0 = cycle_map(xstar, pars)
57 # Start on the same phase with a transverse state perturbation.
58 y = xstar + np.array([perturb, -0.7*perturb])
59 norms = []
60 for _ in range(cycles+1):
61 norms.append(float(np.linalg.norm(y-xstar)))
62 y = cycle_map(y, pars)
63 # fit log norm, excluding any underflow tail
64 vals = np.maximum(np.asarray(norms), 1e-300)
65 slope = float(np.polyfit(np.arange(len(vals)), np.log(vals), 1)[0])
66 return norms, slope
67
68def noise_gain(pars, sigma=1e-5, trials=200, cycles=80):
69 # IID gradient perturbation added at each update; compare stationary RMS.
70 xstar = periodic_orbit(pars)
71 rng = np.random.default_rng(SEED+7)
72 vals=[]
73 for t in range(trials):
74 y=xstar.copy()
75 for k in range(cycles*pars['P']):
76 phi=2*math.pi*(k % pars['P'])/pars['P']
77 # forcing is zero here; inject noise into gradient via equivalent forcing
78 q=dict(pars); q['forcing']=float(rng.normal(0,sigma))
79 y=step(y,phi, **{z:q[z] for z in ('eta0','amp','mu0','mu_amp','curvature','forcing')})
80 vals.append(np.linalg.norm(y-xstar))
81 return float(np.sqrt(np.mean(np.square(vals))))
82
83def stability_sweep():
84 base={'P':12,'eta0':0.8,'amp':0.0,'mu0':0.75,'mu_amp':0.0,'curvature':1.0,'forcing':0.0}
85 # Prediction 1: analytic one-cycle Jacobian and finite differences agree.
86 rows=[]
87 for amp in [0.0,0.2,0.5,0.8]:
88 p=dict(base,amp=amp,mu_amp=0.08)
89 M=monodromy(p); rho=max(abs(np.linalg.eigvals(M)))
90 eps=1e-7; x=periodic_orbit(p); c0=cycle_map(x,p)
91 fd=np.column_stack([(cycle_map(x+eps*np.eye(2)[j],p)-c0)/eps for j in range(2)])
92 rho_fd=max(abs(np.linalg.eigvals(fd)))
93 rows.append({'amp':amp,'rho_analytic':float(rho),'rho_finite_difference':float(rho_fd),
94 'relative_error':float(abs(rho-rho_fd)/max(rho,1e-15))})
95
96 # Prediction 2: rho<1 gives decay and rho>1 gives growth. Find predicted
97 # boundary (rho=1) by bisection, then test both sides empirically.
98 def rho_at(eta):
99 p=dict(base,eta0=float(eta),amp=0.65,mu_amp=0.10)
100 return float(max(abs(np.linalg.eigvals(monodromy(p)))))
101 lo,hi=2.3,2.4
102 for _ in range(50):
103 mid=(lo+hi)/2
104 if rho_at(mid)<1: lo=mid
105 else: hi=mid
106 predicted_boundary=(lo+hi)/2
107 crossing=[]
108 for eta in [predicted_boundary-0.05,predicted_boundary-0.01,
109 predicted_boundary+0.01,predicted_boundary+0.05]:
110 p=dict(base,eta0=float(eta),amp=0.65,mu_amp=0.10)
111 rho=rho_at(eta); norms,slope=empirical_decay(p,cycles=12)
112 crossing.append({'eta0':float(eta),'rho':rho,
113 'predicted_log_multiplier':float(math.log(rho)),
114 'observed_log_decay_per_cycle':slope,
115 'decays':bool(norms[-1] < norms[0])})
116
117 # Prediction 3: near the stable cycle, noise amplification grows with the
118 # resolvent scale 1/(1-rho). Vary momentum to span distinct rho values.
119 ng=[]
120 for mu in [0.10,0.25,0.40,0.55,0.70,0.82]:
121 p=dict(base,eta0=0.35,amp=0.35,mu0=mu,mu_amp=0.0)
122 rho=float(max(abs(np.linalg.eigvals(monodromy(p)))))
123 if rho < .98:
124 gain=noise_gain(p,sigma=2e-5,trials=150,cycles=80)
125 ng.append({'mu0':mu,'rho':rho,'noise_rms':gain,
126 'predicted_resolvent':float(1/(1-rho))})
127 # rank correlation is a scale-free check of monotonicity.
128 order=np.argsort([x['rho'] for x in ng])
129 noise_monotone=all(ng[order[i]]['noise_rms'] <= ng[order[i+1]]['noise_rms']
130 for i in range(len(order)-1))
131 return rows,crossing,ng,{'predicted_eta_boundary_rho1':predicted_boundary,
132 'noise_monotone_with_rho':bool(noise_monotone)}
133
134# Small nonlinear regression: same data, update count, and average LR.
135def mlp_compare():
136 rng=np.random.default_rng(SEED)
137 X=rng.normal(size=(512,2)).astype(np.float64)
138 y=(np.sin(X[:,0])+0.35*X[:,1]**2).astype(np.float64)
139 split=384; Xtr,ytr=X[:split],y[:split]; Xte,yte=X[split:],y[split:]
140 def run(cyclic):
141 W1=rng.normal(0,.35,(2,16)); b1=np.zeros(16); W2=rng.normal(0,.2,16); b2=0.
142 m=[np.zeros_like(W1),np.zeros_like(b1),np.zeros_like(W2),0.]
143 eta0=.025; mu0=.85; P=20; amp=.55; muamp=.06
144 for k in range(1800):
145 ix=rng.choice(split,64,replace=False); a=Xtr[ix]; target=ytr[ix]
146 h=np.tanh(a@W1+b1); pred=h@W2+b2; d=(2*(pred-target)/len(ix))
147 dW2=h.T@d; db2=d.sum(); dh=d[:,None]*W2[None,:]; dz=dh*(1-h*h)
148 dW1=a.T@dz; db1=dz.sum(0)
149 if cyclic:
150 ph=2*math.pi*(k%P)/P; eta=eta0*(1+amp*math.sin(ph)); mu=mu0+muamp*math.cos(ph)
151 else: eta=eta0; mu=mu0
152 for j,g in enumerate([dW1,db1,dW2,db2]):
153 m[j]=mu*m[j]+g
154 W1-=eta*m[0]; b1-=eta*m[1]; W2-=eta*m[2]; b2-=eta*m[3]
155 pred=np.tanh(Xte@W1+b1)@W2+b2
156 return float(np.mean((pred-yte)**2))
157 # reset is intentionally deterministic via local seeds for fair pair.
158 np.random.seed(SEED); random.seed(SEED); cyc=run(True)
159 np.random.seed(SEED); random.seed(SEED); const=run(False)
160 return {'constant_average_momentum_test_mse':const,'cyclic_test_mse':cyc}
161
162def main():
163 rows,crossing,ng,summary=stability_sweep()
164 result={'seed':SEED,'predictions':[
165 'Floquet multiplier from the analytic one-cycle Jacobian equals finite-difference cycle response.',
166 'The decay/growth transition occurs at spectral radius rho=1.',
167 'Stable noisy response increases with the predicted resolvent scale 1/(1-rho).'],
168 'floquet_check':rows,'stability_boundary_sweep':crossing,'noise_scaling':ng,'summary':summary,
169 'mlp_comparison':mlp_compare()}
170 Path('results.json').write_text(json.dumps(result,indent=2))
171 print(json.dumps(result,indent=2))
172if __name__=='__main__': main()