Koopman Skew-Dilation RNN / stage2_bench.py

Mechanism confirmed, baseline not beaten

Raw ⬇ ZIP
 1import sys, json, random, numpy as np, torch
 2import torch.nn as nn
 3sys.path.insert(0, '/home/maxwelhelp/all/math2nn')
 4from bench import get_dataset, train_model, evaluate, sweep_baseline, make_report
 5
 6SEEDS = tuple(range(8))
 7# Union used by both sides; baseline sweep includes all idea learning rates.
 8GRID = [{'lr': 1e-3, 'epochs': 12}, {'lr': 3e-3, 'epochs': 12}, {'lr': 1e-2, 'epochs': 12}]
 9
10class LiftRNN(nn.Module):
11    """Finite m=2 skew-dilation recurrent cell with Cayley propagation."""
12    def __init__(self, input_dim=3, out_dim=1, hidden=64, eta=.05):
13        super().__init__(); self.d=hidden; self.eta=eta
14        self.rawS=nn.Parameter(torch.randn(hidden,hidden)*.02)
15        self.rawA=nn.Parameter(torch.randn(hidden,hidden)*.02)
16        self.inp=nn.Linear(input_dim, 2*hidden)
17        self.head=nn.Linear(hidden,out_dim)
18    def forward(self,x):
19        x=x.view(x.shape[0],-1,3)
20        B=x.shape[0]; dev=x.device
21        S=self.rawS+self.rawS.T; A=self.rawA-self.rawA.T
22        z0=torch.zeros((2*self.d,),device=dev,dtype=x.dtype)
23        L=torch.zeros((2*self.d,2*self.d),device=dev,dtype=x.dtype)
24        L[:self.d,:self.d]=A; L[:self.d,self.d:]=S
25        L[self.d:,:self.d]=-S; L[self.d:,self.d:]=A
26        I=torch.eye(2*self.d,device=dev,dtype=x.dtype)
27        C=torch.linalg.solve(I-self.eta*L/2,I+self.eta*L/2)
28        z=torch.zeros((B,2*self.d),device=dev,dtype=x.dtype)
29        for k in range(x.shape[1]):
30            z=z@C.T + self.inp(x[:,k])
31            # bounded nonlinear readout is only the task head; recurrent lift remains linear/unitary
32        return self.head(z[:,:self.d])
33
34def baseline_fn(cfg, seed):
35    torch.manual_seed(seed); np.random.seed(seed); random.seed(seed)
36    return make_base(3,1)
37def make_base(input_dim=3,out_dim=1):
38    class Base(nn.Module):
39        def __init__(self):
40            super().__init__(); self.rnn=nn.GRU(input_dim,64,batch_first=True); self.head=nn.Linear(64,out_dim)
41        def forward(self,x):
42            _,h=self.rnn(x.view(x.shape[0],-1,input_dim)); return self.head(h[-1])
43    return Base()
44def idea_fn(cfg, seed):
45    torch.manual_seed(seed); np.random.seed(seed); random.seed(seed)
46    return LiftRNN(3,1,64,eta=.05)
47
48def run_one(fn, cfg, seed):
49    ds=get_dataset('dynamics',seed,n_train=400,n_test=200)
50    net,metric,hist=train_model(fn(cfg,seed),ds,epochs=cfg['epochs'],lr=cfg['lr'],batch=128)
51    if net is None: raise RuntimeError('training failed')
52    return metric, net, ds
53
54def eval_cfg(fn,cfg,seeds=SEEDS):
55    vals=[]
56    for s in seeds: vals.append(run_one(fn,cfg,s)[0])
57    return {'per_seed':vals,'mean':float(np.mean(vals)),'std':float(np.std(vals,ddof=1))}
58
59def signature():
60    # measured on a trained benchmark model: compare Cayley norm preservation to Euler growth
61    metric,net,ds=run_one(idea_fn,GRID[1],0); net.eval(); d=net.d
62    with torch.no_grad():
63        S=net.rawS+net.rawS.T; A=net.rawA-net.rawA.T
64        L=torch.zeros((2*d,2*d)); L[:d,:d]=A; L[:d,d:]=S; L[d:,:d]=-S; L[d:,d:]=A
65        I=torch.eye(2*d); eta=net.eta
66        C=torch.linalg.solve(I-eta*L/2,I+eta*L/2); E=torch.randn(2*d)
67        zc=E.clone(); ze=E.clone()
68        for _ in range(20): zc=C@zc; ze=(I+eta*L)@ze
69        cay=float((zc.norm()/E.norm()-1).abs()); euler=float(ze.norm()/E.norm())
70    return {'prediction':'Cayley preserves lifted norm; explicit Euler grows it',
71            'observed_cayley_relative_norm_error':cay,
72            'observed_euler_20_step_norm_ratio':euler,
73            'trained_test_mse':float(metric),'confirmed':bool(cay<1e-5 and euler>1.0001)}
74
75def main():
76    # baseline sweep is explicitly over the same union as the idea settings
77    # Invoke the canonical harness sweep: make_fn(cfg) returns train_fn(seed).
78    from bench import sweep_baseline
79    def tuned_make(cfg):
80        def train_fn(seed):
81            return run_one(baseline_fn, cfg, int(seed))[0]
82        return train_fn
83    tuned=sweep_baseline(tuned_make, GRID, seeds=tuple(range(4)))
84    sweep=[]
85    for cfg in GRID:
86        r=eval_cfg(baseline_fn,cfg); sweep.append({'cfg':cfg,**r})
87    best=min(sweep,key=lambda r:r['mean'])
88    base={'best_cfg':best['cfg'],'sweep':sweep,'harness_tuning':tuned,'full':eval_cfg(baseline_fn,best['cfg'])}
89    idea=[]
90    for cfg in GRID:
91        r=eval_cfg(idea_fn,cfg); idea.append({'cfg':cfg,**r})
92    ibest=min(idea,key=lambda r:r['mean']); idea_res={k:ibest[k] for k in ('per_seed','mean','std')}
93    rep=make_report('dynamics','rnn_small',base,idea_res,{'idea_sweep':idea,'mechanism_signature':signature()})
94    rep['custom_track']=None
95    with open('bench_report.json','w') as f: json.dump(rep,f,indent=2)
96    print(json.dumps(rep,indent=2))
97if __name__=='__main__': main()