import json, math, random, time from pathlib import Path import numpy as np SEED = 1552 np.random.seed(SEED) random.seed(SEED) def eig_sensitivity(A, B): """d lambda/dp for A(p)=A+pB, using the stated left/right formula.""" vals, vr = np.linalg.eig(A) vals_l, vl = np.linalg.eig(A.T.conj()) out = [] for i, lam in enumerate(vals): # left eigenvector w satisfies w^H A=lambda w^H; column eigvec of A^H j = int(np.argmin(np.abs(vals_l - np.conj(lam)))) v = vr[:, i] w = vl[:, j] d = (np.conj(w) @ B @ v) / (np.conj(w) @ v) out.append(d) return vals, np.asarray(out) def radial_slope(lam, dlam): if abs(lam) < 1e-12: return abs(dlam) return float(np.real(np.exp(-1j*np.angle(lam))*dlam)) def toy_verification(): # A(p)=A0+pB has a known dominant mode and a nontrivial off-diagonal coupling. # It makes the first-order crossing prediction directly falsifiable. rows = [] for margin in [0.05, 0.10, 0.20, 0.35]: for gain in [0.25, 0.50, 1.0, 1.75]: a = 1.0 - margin A = np.array([[a, 0.65], [0.0, 0.42]], dtype=float) B = np.array([[gain, -0.3], [0.15, 0.05]], dtype=float) vals, ds = eig_sensitivity(A, B) i = int(np.argmax(np.abs(vals))) lam, dlam = vals[i], ds[i] slope = radial_slope(lam, dlam) pred = (1-abs(lam))/slope if slope > 1e-10 else np.inf # exact observed crossing by a dense one-dimensional sweep ps = np.linspace(0, min(2.0, max(.02, 2.5*pred)), 20001) radii = np.array([max(abs(np.linalg.eigvals(A+p*B))) for p in ps]) ix = np.where(radii >= 1.0)[0] obs = float(ps[ix[0]]) if len(ix) else float('nan') # independently check derivative with centered finite differences eps = 1e-6 fd = (np.linalg.eigvals(A+eps*B)[i]-np.linalg.eigvals(A-eps*B)[i])/(2*eps) relsens = abs(dlam-fd)/max(1e-12, abs(fd)) rows.append(dict(margin=margin, gain=gain, predicted_radius=float(pred), observed_radius=obs, relative_prediction_error=(abs(obs-pred)/pred if np.isfinite(obs) else None), sensitivity=float(abs(dlam)), finite_difference_relative_error=float(relsens))) # Three explicit predictions: formula correctness, linear gain scaling, and margin scaling. formula_err = max(r['finite_difference_relative_error'] for r in rows) valid = [r for r in rows if r['observed_radius'] == r['observed_radius']] cross_err = float(np.median([r['relative_prediction_error'] for r in valid])) # Since derivative is gain-proportional, delta_hat * gain should be constant at fixed margin. scaled = [] for m in [0.05, .10, .20, .35]: x = [r['predicted_radius']*r['gain'] for r in rows if r['margin']==m] scaled.append(float(np.std(x)/np.mean(x))) gain_scaling = float(np.median(scaled)) # delta_hat is proportional to margin for this affine family. margin_ratios = [] for g in [.25, .5, 1., 1.75]: x = [r['predicted_radius']/r['margin'] for r in rows if r['gain']==g] margin_ratios.append(float(np.std(x)/np.mean(x))) margin_scaling = float(np.median(margin_ratios)) return rows, dict(max_formula_relative_error=formula_err, median_crossing_relative_error=cross_err, gain_scaling_cv=gain_scaling, margin_scaling_cv=margin_scaling, predictions={'eigen_derivative_matches_finite_difference': formula_err < 1e-5, 'crossing_radius_is_margin_over_slope': cross_err < .03, 'radius_scales_inverse_linearly_with_gain': gain_scaling < .02, 'radius_scales_linearly_with_margin': margin_scaling < .02}) def train_rnn(regularized, device='cpu', steps=350): import torch torch.manual_seed(SEED + int(regularized)) D, T, B = 12, 25, 64 Wh = torch.nn.Parameter(torch.randn(D,D,device=device)*0.28) Wx = torch.nn.Parameter(torch.randn(D,1,device=device)*0.35) Wo = torch.nn.Parameter(torch.randn(1,D,device=device)*0.15) b = torch.nn.Parameter(torch.zeros(D,device=device)) params=[Wh,Wx,Wo,b] opt=torch.optim.Adam(params, lr=.012) history=[] def jacobian(h, u, gain=1.0): # Jacobian of one state update at one representative operating point. def fn(x): return torch.tanh(gain*Wh@x + Wx@u + b) return torch.autograd.functional.jacobian(fn, h, create_graph=True) for step in range(steps): u=torch.randn(B,T,1,device=device) # Stable teacher: exponentially filtered input, trained as sequence predictor. y=torch.zeros(B,T,1,device=device) for t in range(1,T): y[:,t]=.82*y[:,t-1]+.18*u[:,t] h=torch.zeros(B,D,device=device); loss=0. for t in range(T): h=torch.tanh(h@Wh.T + u[:,t]@Wx.T + b) loss=loss+torch.mean((h@Wo.T-y[:,t])**2) loss=loss/T reg=torch.tensor(0.,device=device) if regularized and step % 8 == 0: # The perturbation is recurrent gain p. Penalize predicted/actual local # spectral radius at p=0 and p=+0.08, a differentiable proxy for margin. h0=h[0].detach().requires_grad_(True); u0=u[0,-1].detach() A=jacobian(h0,u0,1.0) Ap=jacobian(h0,u0,1.08) rho=torch.max(torch.abs(torch.linalg.eigvals(A))) rhop=torch.max(torch.abs(torch.linalg.eigvals(Ap))) # margin target and sensitivity-aware forward perturbation penalty reg=0.8*torch.relu(rho-.88)**2 + 1.5*torch.relu(rhop-.93)**2 total=loss+reg opt.zero_grad(); total.backward(); torch.nn.utils.clip_grad_norm_(params, 1.0); opt.step() if step % 50 == 0: history.append(float(loss.detach().cpu())) return [x.detach().cpu() for x in params], history def evaluate(params, device='cpu'): import torch Wh,Wx,Wo,b=params Wh,Wx,Wo,b=[x.to(device) for x in params] D=Wh.shape[0] def run(gain, noise, T=180): h=torch.zeros(1,D,device=device); norms=[]; outs=[] u=torch.randn(1,T,1,device=device)*noise for t in range(T): h=torch.tanh(gain*(h@Wh.T)+u[:,t]@Wx.T+b) norms.append(float(torch.linalg.vector_norm(h).item())); outs.append(float((h@Wo.T).item())) return max(norms), float(np.mean(np.square(outs[-40:]))) # local eigenvalue and first-order prediction for gain perturbation, via numpy fd h=torch.zeros(D,device=device,requires_grad=True); u=torch.zeros(1,device=device) def fn(x,g): return torch.tanh(g*(Wh@x)+Wx@u+b) A=torch.autograd.functional.jacobian(lambda x:fn(x,1.),h).detach().cpu().numpy() Ap=torch.autograd.functional.jacobian(lambda x:fn(x,1.001),h).detach().cpu().numpy() vals=np.linalg.eigvals(A); rho=max(abs(vals)); dr=(max(abs(np.linalg.eigvals(Ap)))-rho)/.001 pred=(1-rho)/max(1e-9,dr) return dict(rho=float(rho), predicted_radius=float(pred), gain1_noise1=run(1.,1.), gain1p_noise1=run(1.08,1.), gain1_noise2=run(1.,2.)) def main(): rows,toy=toy_verification() result={'seed':SEED,'toy_summary':toy,'toy_rows':rows} try: import torch device='cuda' if torch.cuda.is_available() else 'cpu' try: base,hb=train_rnn(False,device); idea,hi=train_rnn(True,device) except Exception as e: device='cpu'; base,hb=train_rnn(False,device); idea,hi=train_rnn(True,device) result.update({'device':device,'baseline':evaluate(base,device),'idea':evaluate(idea,device), 'training_loss_samples':{'baseline':hb,'idea':hi}}) except Exception as e: result['rnn_error']=repr(e) Path('results.json').write_text(json.dumps(result,indent=2)) print(json.dumps(result,indent=2)) if __name__=='__main__': main()