import sys, json, math, time from pathlib import Path import numpy as np import torch from torch import nn sys.path.insert(0, '/home/maxwelhelp/all/math2nn') from bench import get_dataset, train_model, sweep_baseline, make_report SEED0=3104 NFEAT=48 BATCH=128 EPOCHS=18 class RFFMLP(nn.Module): def __init__(self, xtr, mode='gaussian', kappa=0.75, whiten=True, seed=0): super().__init__() x = xtr.detach().float() self.register_buffer('mu', x.mean(0)) self.register_buffer('sd', x.std(0).clamp_min(1e-3)) d=x.shape[1] rng=np.random.default_rng(seed) if mode=='gaussian': w=rng.normal(0.0, 1.0, (NFEAT,d)).astype('float32') elif mode=='matched': # exp(-2*kappa*r), including the d-dimensional Jacobian: r is Gamma(d, rate=2*kappa) r=rng.gamma(shape=d, scale=1.0/(2.0*kappa), size=NFEAT) u=rng.normal(size=(NFEAT,d)); u/=np.linalg.norm(u,axis=1,keepdims=True) w=(r[:,None]*u).astype('float32') else: raise ValueError(mode) b=rng.uniform(0,2*np.pi,NFEAT).astype('float32') self.register_buffer('w', torch.from_numpy(w)); self.register_buffer('b',torch.from_numpy(b)) # Fit the whitening transform on training inputs only, as a frozen input projection. z=np.sqrt(2.)*np.cos(((x.numpy()-self.mu.numpy())/self.sd.numpy())@w.T+b[None,:]) g=(z.T@z)/len(z); eps=1e-5*float(np.trace(g))/NFEAT vals, vecs=np.linalg.eigh(g+eps*np.eye(NFEAT)); vals=np.maximum(vals,eps) T=(vecs*(1/np.sqrt(vals)))@vecs.T if whiten else np.eye(NFEAT) self.register_buffer('T',torch.from_numpy(T.astype('float32'))) self.head=nn.Sequential(nn.Linear(NFEAT,64),nn.Tanh(),nn.Linear(64,64),nn.Tanh(),nn.Linear(64,1)) self.raw_cond=float(np.linalg.eigvalsh(g)[-1]/max(np.linalg.eigvalsh(g)[0],1e-12)) self.radius=float(np.linalg.norm(w,axis=1).mean()) def forward(self,x): q=(x-self.mu)/self.sd z=math.sqrt(2.)*torch.cos(q@self.w.T+self.b) return self.head(z@self.T) def make_runner(cfg, mode): def run(seed, capture=False): d=get_dataset('tabular',seed,n_train=400,n_test=200) torch.manual_seed(SEED0+seed); np.random.seed(SEED0+seed) net=RFFMLP(d['xtr'],mode=mode,kappa=cfg.get('kappa',.75),whiten=True,seed=SEED0+100*seed+int(cfg.get('kappa',.75)*100)) net, metric, hist=train_model(net,d,epochs=EPOCHS,lr=cfg['lr'],batch=BATCH,weight_decay=cfg.get('weight_decay',0.0),log=lambda *a,**k:None) if not capture: return float(metric) with torch.no_grad(): dev=next(net.parameters()).device xte=d['xte'].to(dev); yte=d['yte'].to(dev) pred=net(xte).cpu().numpy().ravel(); y=yte.cpu().numpy().ravel() # NN-scale mechanism check: finite-difference output sensitivity versus sampled spectral radius. x=xte[:128]; delta=torch.zeros_like(x); delta[:,0]=1e-3 sens=float(((net(x+delta)-net(x))/1e-3).abs().mean().item()) return {'metric':float(metric),'radius':net.radius,'raw_cond':net.raw_cond,'sensitivity':sens,'pred':pred,'y':y} return run def main(): t=time.time(); seeds=tuple(range(8)) # Same lr union on both sides; whitening is fixed as the explicitly proposed stabilizer. lrs=[0.0015,0.003,0.006] grid=[{'lr':v} for v in lrs] base=sweep_baseline(lambda c: make_runner(c,'gaussian'),grid,seeds=seeds) # Idea sweep has same size and exact lr union, with a priori regularity values. idea_trials=[] # Three idea settings: baseline-best lr plus two nearby lrs; kappa is fixed # a priori to the smooth Gevrey choice used by the implementation. for lr in lrs: cfg={'lr':lr,'kappa':.5} vals=[make_runner(cfg,'matched')(s) for s in seeds] idea_trials.append({'cfg':cfg,'mean':float(np.mean(vals)),'per_seed':vals}) best=min(idea_trials,key=lambda z:z['mean']); idea_cfg=best['cfg'] idea_vals=[make_runner(idea_cfg,'matched')(s) for s in seeds] # Full baseline best config, explicitly paired with the same seeds. base_vals=[make_runner(base['best_cfg'],'gaussian')(s) for s in seeds] diffs=[float(a-b) for a,b in zip(idea_vals,base_vals)] # Reuse harness report; it computes paired delta and sign permutation p-value. base_block={'best_cfg':base['best_cfg'],'sweep':base['sweep'],'full':{'per_seed':base_vals,'mean':float(np.mean(base_vals)),'std':float(np.std(base_vals))}} idea_res={'cfg':idea_cfg,'sweep':idea_trials,'per_seed':idea_vals,'mean':float(np.mean(idea_vals)),'std':float(np.std(idea_vals))} # trained behavior signature, predicted inequality is lower radius/sensitivity for matched. bs=[make_runner(base['best_cfg'],'gaussian')(s,True) for s in seeds] ins=[make_runner(idea_cfg,'matched')(s,True) for s in seeds] observed_radius_gap=float(np.mean([z['radius'] for z in ins])-np.mean([z['radius'] for z in bs])) observed_sensitivity_gap=float(np.mean([z['sensitivity'] for z in ins])-np.mean([z['sensitivity'] for z in bs])) # Quantitative prediction made before inspecting trained models: in d=10, # exp(-2*kappa*r), kappa=.5 gives E[r]=d/(2*kappa)=10, versus about sqrt(d) # for N(0,I); larger frequencies should therefore produce larger input slopes. predicted_radius_gap=10.0-math.sqrt(10.0) predicted_sensitivity_gap=predicted_radius_gap sig={'prediction':'in standardized d=10 inputs, the matched exponential radial law (kappa=.5) predicts larger mean frequency radius than Gaussian and consequently larger input sensitivity','predicted_radius_gap':predicted_radius_gap,'observed_radius_gap':observed_radius_gap,'predicted_sensitivity_gap':predicted_sensitivity_gap,'observed_sensitivity_gap':observed_sensitivity_gap,'baseline_radius_mean':float(np.mean([z['radius'] for z in bs])),'idea_radius_mean':float(np.mean([z['radius'] for z in ins])),'baseline_sensitivity_mean':float(np.mean([z['sensitivity'] for z in bs])),'idea_sensitivity_mean':float(np.mean([z['sensitivity'] for z in ins])),'confirmed':bool(observed_radius_gap>0 and observed_sensitivity_gap>0)} rep=make_report('tabular','mlp_tiny',base_block,idea_res,{'mechanism_signature':sig,'protocol':'8 paired seeds; baseline and idea share frozen-RFF+MLP architecture; baseline lr union equals idea lr union','runtime_sec':time.time()-t}) Path('bench_report.json').write_text(json.dumps(rep,indent=2)); print(json.dumps(rep,indent=2)) if __name__=='__main__': main()