Record the matching battery: fifteen solvers, one amplifier

Adds MATCHING_RESULTS.md. Entropic GW annealed and PATH both reach 0.961
on the reference field, above the 0.958 the project's own pipeline
reached after months, and neither had been run. GW reaches 0.881 on the
rank-8 field where GRAMPA reaches 0.076 -- which retires the morning's
rank-ladder conclusion, since that ladder was run entirely with GRAMPA
and GRAMPA degrades on clustered eigenvalues.

The hard instance yields to amplification rather than a better solver:
descent multiplies a partial answer by about four, so the job is to feed
it a start that is 10-20% correct rather than to replace it.

Co-Authored-By: Claude <noreply@anthropic.com>

14 files changed, 460 insertions, 0 deletions

diff --git a/artifacts/spectral_frontier_probe/band.py b/artifacts/spectral_frontier_probe/band.py
new file mode 100644
index 0000000..cca368a
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/band.py
@@ -0,0 +1,26 @@
+import numpy as np, torch
+from scipy.optimize import linear_sum_assignment
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V); truth=np.arange(N)
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1]
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+acc=lambda p: float((p==truth).mean())
+CONST=float((T*T).sum()+(V*V).sum())
+en=lambda p:(CONST-2*(T[np.ix_(p,p)]*V).sum())/(N*(N-1))
+for r in (12,16,24):
+    XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r]))
+    u,s,vt=np.linalg.svd(XT.T@XV); O=u@vt
+    idx=np.arange(r); band_mass=[]
+    for b in (1,2,3,4,6,r):
+        Mk=(np.abs(idx[:,None]-idx[None,:])<=b).astype(float)
+        frac=float((O**2*Mk).sum()/ (O**2).sum())
+        Ob=O*Mk; uu,ss,vv=np.linalg.svd(Ob); Ob=uu@vv
+        p=hung(XV,XT@Ob)
+        band_mass.append((b,round(frac,3),round(acc(p),3),round(en(p),4)))
+    print(f"r={r}: (bandwidth, mass of oracle-O inside band, acc after re-orthogonalising, E) -> {band_mass}")
+    print(f"     free params: full {r*(r-1)//2}, band3 {sum(min(3,r-1-i) for i in range(r))}")
diff --git a/artifacts/spectral_frontier_probe/confirm.log b/artifacts/spectral_frontier_probe/confirm.log
new file mode 100644
index 0000000..ff6ac9d
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/confirm.log
@@ -0,0 +1,7 @@
+  r=  4: pre-descent best E 0.7185 (acc 0.004) -> post-descent E 0.6493 acc 0.016
+  r=  6: pre-descent best E 0.5638 (acc 0.117) -> post-descent E 0.5333 acc 0.148
+  r=  8: pre-descent best E 0.4786 (acc 0.395) -> post-descent E 0.3597 acc 0.750
+  r= 10: pre-descent best E 0.6224 (acc 0.180) -> post-descent E 0.5519 acc 0.172
+  r= 12: pre-descent best E 0.3977 (acc 0.660) -> post-descent E 0.3396 acc 0.879
+  r= 14: pre-descent best E 0.3671 (acc 0.707) -> post-descent E 0.3396 acc 0.836
+  r= 16: pre-descent best E 0.6410 (acc 0.145) -> post-descent E 0.5474 acc 0.277
diff --git a/artifacts/spectral_frontier_probe/confirm.py b/artifacts/spectral_frontier_probe/confirm.py
new file mode 100644
index 0000000..d0fc003
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/confirm.py
@@ -0,0 +1,49 @@
+import sys, time, torch, numpy as np
+sys.path.insert(0,'/home/yurenh2/emm')
+from scipy.optimize import linear_sum_assignment
+from scipy.stats import ortho_group
+from worldalign.synth_fast_gate import fast_pair_descent
+dev='cuda:3'
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(T,dtype=torch.float32,device=dev)
+CONST=float((Tt*Tt).sum()+(Vt*Vt).sum())
+def energy(p):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1))
+def descend(p,steps=4000):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy()
+truth=np.arange(N); acc=lambda p: float((p==truth).mean())
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1]
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+rng=np.random.default_rng(20260801)
+def icp(XV,XT,O,iters=30):
+    for _ in range(iters):
+        p=hung(XV,XT@O)
+        u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt
+        if np.allclose(On,O,atol=1e-10): O=On; break
+        O=On
+    return hung(XV,XT@O)
+t0=time.time(); pool=[]
+for r in (4,6,8,10,12,14,16,20,24):
+    XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r]))
+    cand=[]
+    for t in range(120):
+        p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30))))
+        cand.append((energy(p),p))
+    cand.sort(key=lambda z:z[0])
+    best=None
+    for e,p in cand[:5]:
+        pd=descend(p); ed=energy(pd)
+        if best is None or ed<best[0]: best=(ed,pd)
+    pool.append((best[0],best[1],r))
+    print(f"  r={r:3d}: pre-descent best E {cand[0][0]:.4f} (acc {acc(cand[0][1]):.3f}) -> post-descent E {best[0]:.4f} acc {acc(best[1]):.3f}",flush=True)
+pool.sort(key=lambda z:z[0])
+print(f"\nBLIND PICK (lowest E over all r): r={pool[0][2]}  E={pool[0][0]:.4f}  ACC={acc(pool[0][1]):.3f}")
+print(f"E(truth)={energy(truth):.4f}   total wall {time.time()-t0:.0f}s")
diff --git a/artifacts/spectral_frontier_probe/diag10.py b/artifacts/spectral_frontier_probe/diag10.py
new file mode 100644
index 0000000..752612c
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag10.py
@@ -0,0 +1,41 @@
+import numpy as np, torch
+from collections import defaultdict
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy()
+T=standardise(T0); V=standardise(V0); N=len(T)
+# exact automorphism-by-identical-rows of T (blind: uses T only)
+key={}; cls=defaultdict(list)
+Toff=T.copy(); np.fill_diagonal(Toff,0.0)
+# group scenes whose T rows agree after removing the two swapped coords
+groups=[]; used=np.zeros(N,bool)
+for i in range(N):
+    if used[i]: continue
+    g=[i]; used[i]=True
+    for j in range(i+1,N):
+        if used[j]: continue
+        a=np.delete(Toff[i],[i,j]); b=np.delete(Toff[j],[i,j])
+        if np.abs(a-b).max()<1e-9 and abs(Toff[i,j]-max(Toff[i,i],0))<1e6:
+            g.append(j); used[j]=True
+    groups.append(g)
+sizes=np.array([len(g) for g in groups])
+print("T exact-twin classes: total",len(groups),"; size histogram",np.bincount(sizes)[1:])
+excess=int((sizes-1).sum())
+print("scenes in non-trivial classes:",int(sizes[sizes>1].sum()),"; log|Aut| classes:",int((sizes>1).sum()))
+ceil=(N-int(sizes[sizes>1].sum())+int((sizes>1).sum()))/N
+print(f"blind accuracy ceiling if class is identified but member picked at random: {ceil:.3f}")
+# check the objective really is invariant: swap two members of a class
+import itertools
+p=np.arange(N)
+E0=(( (T*T).sum()+(V*V).sum() )-2*(T[np.ix_(p,p)]*V).sum())/(N*(N-1))
+bad=0; tested=0
+for g in groups:
+    if len(g)>1:
+        i,j=g[0],g[1]; q=p.copy(); q[[i,j]]=q[[j,i]]
+        E1=(((T*T).sum()+(V*V).sum())-2*(T[np.ix_(q,q)]*V).sum())/(N*(N-1))
+        tested+=1
+        if abs(E1-E0)>1e-9: bad+=1
+print(f"swapping twins changes E in {bad}/{tested} classes (0 means exact symmetry of the objective)")
+print(f"E(truth) = {E0:.10f}")
diff --git a/artifacts/spectral_frontier_probe/diag11.log b/artifacts/spectral_frontier_probe/diag11.log
new file mode 100644
index 0000000..dcdcefe
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag11.log
@@ -0,0 +1,8 @@
+r=24 symNMF resid V 0.043 T 0.034 [325s]
+   oracle col cos: [0.983 0.983 0.961 0.957 0.955 0.955 0.951 0.949 0.934 0.908]
+   BLIND stat col-match agrees w/ oracle 0.12 -> scene acc 0.020
+   after co-alternation (1 iters): col agree 0.12 -> scene acc 0.020
+r=42 symNMF resid V 0.023 T 0.014 [741s]
+   oracle col cos: [0.973 0.949 0.948 0.943 0.936 0.921 0.912 0.901 0.892 0.89 ]
+   BLIND stat col-match agrees w/ oracle 0.00 -> scene acc 0.027
+   after co-alternation (2 iters): col agree 0.02 -> scene acc 0.012
diff --git a/artifacts/spectral_frontier_probe/diag11.py b/artifacts/spectral_frontier_probe/diag11.py
new file mode 100644
index 0000000..4977c3c
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag11.py
@@ -0,0 +1,52 @@
+import numpy as np, torch, warnings, time
+warnings.filterwarnings('ignore')
+from scipy.optimize import linear_sum_assignment
+from sklearn.decomposition import NMF
+np.set_printoptions(precision=3,suppress=True,linewidth=200)
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy(); N=len(V0)
+rng=np.random.default_rng(777); sigma=rng.permutation(N)
+T0s=T0[np.ix_(sigma,sigma)]; truth=np.argsort(sigma)
+acc=lambda p: float((p==truth).mean())
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+def sym_nmf(M,r,iters=3000,seed=0):
+    """symmetric NMF  M ~ W W^T  by multiplicative updates (Ding et al.)"""
+    M=np.clip(M,0,None).copy(); np.fill_diagonal(M,np.clip(np.diag(M),0,None))
+    rg=np.random.default_rng(seed); W=np.abs(rg.standard_normal((len(M),r)))*np.sqrt(M.mean()/r)
+    for _ in range(iters):
+        num=M@W; den=W@(W.T@W)+1e-12
+        W=W*(0.5+0.5*num/den)
+    return W
+for r in (24,42):
+    t0=time.time()
+    WV=sym_nmf(V0,r,seed=1); WT=sym_nmf(T0s,r,seed=1)
+    print(f"r={r} symNMF resid V {np.linalg.norm(V0-WV@WV.T)/np.linalg.norm(V0):.3f} T {np.linalg.norm(T0s-WT@WT.T)/np.linalg.norm(T0s):.3f} [{time.time()-t0:.0f}s]")
+    A=WV.copy(); B=WT.copy()
+    # ORACLE column match for reference
+    An=A/np.linalg.norm(A,axis=0,keepdims=True); Bn=B/np.linalg.norm(B,axis=0,keepdims=True)
+    rr,cc=linear_sum_assignment(-(An[np.arange(N)].T@Bn[truth]))
+    print("   oracle col cos:",np.sort((An.T@Bn[truth])[rr,cc])[::-1][:10])
+    # BLIND column match by permutation-invariant column signatures
+    def sig(X):
+        Xn=X/ (np.linalg.norm(X,axis=0,keepdims=True)+1e-12)
+        q=np.quantile(Xn,np.linspace(0.5,1.0,16),axis=0).T
+        return np.hstack([q,(Xn>0.02).mean(0)[:,None],np.linalg.norm(X,axis=0)[:,None]/np.linalg.norm(X)])
+    sA=sig(A); sB=sig(B); m=hung(sA,sB)
+    agree=float((m==cc[np.argsort(rr)]).mean())
+    def scene_match(colmap):
+        AV=A; BT=B[:,colmap]
+        AV=AV/np.linalg.norm(AV,axis=1,keepdims=True).clip(1e-9); BT=BT/np.linalg.norm(BT,axis=1,keepdims=True).clip(1e-9)
+        return hung(AV,BT)
+    p=scene_match(m); print(f"   BLIND stat col-match agrees w/ oracle {agree:.2f} -> scene acc {acc(p):.3f}")
+    # alternate: scene Hungarian <-> column Hungarian
+    colmap=m.copy()
+    for it in range(25):
+        p=scene_match(colmap)
+        Ar=A/np.linalg.norm(A,axis=1,keepdims=True).clip(1e-9)
+        Br=B/np.linalg.norm(B,axis=1,keepdims=True).clip(1e-9)
+        rr2,cc2=linear_sum_assignment(-(Ar.T@Br[p])); newmap=cc2[np.argsort(rr2)]
+        if (newmap==colmap).all(): break
+        colmap=newmap
+    p=scene_match(colmap)
+    print(f"   after co-alternation ({it+1} iters): col agree {float((colmap==cc[np.argsort(rr)]).mean()):.2f} -> scene acc {acc(p):.3f}")
diff --git a/artifacts/spectral_frontier_probe/diag3.py b/artifacts/spectral_frontier_probe/diag3.py
new file mode 100644
index 0000000..45a1430
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag3.py
@@ -0,0 +1,21 @@
+import torch, numpy as np
+np.set_printoptions(precision=3, suppress=True, linewidth=200)
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1]
+K=24
+X=np.abs(UV[:,:K].T@UT[:,:K])
+print("|<u_k(V), u_l(T)>| top-16 block (rows=V index, cols=T index):")
+print(X[:16,:16])
+print("\ndiagonal |<u_k(V),u_k(T)>| k=0..23:", np.diag(X))
+print("row-max of |overlap| (best T partner for each V eigvec):", X.max(1))
+print("argmax:", X.argmax(1))
+# subspace alignment: principal angles between top-r spaces
+for r in (4,6,8,10,12,16,20,24,32,48):
+    s=np.linalg.svd(UV[:,:r].T@UT[:,:r],compute_uv=False)
+    print(f"r={r:3d}  mean cos principal angle {s.mean():.3f}  captured energy {np.sum(s**2)/r:.3f}  min {s.min():.3f}")
diff --git a/artifacts/spectral_frontier_probe/diag4.py b/artifacts/spectral_frontier_probe/diag4.py
new file mode 100644
index 0000000..0f7a75a
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag4.py
@@ -0,0 +1,29 @@
+import torch, numpy as np
+from scipy.optimize import linear_sum_assignment
+np.set_printoptions(precision=3, suppress=True, linewidth=200)
+rng=np.random.default_rng(0)
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1]
+truth=np.arange(N)
+def acc(p): return float((p==truth).mean())
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T)
+    r,c=linear_sum_assignment(C); return c
+
+print("### 1. ORACLE orthogonal mixing O between top-r spectral embeddings")
+for r in (6,8,10,12,16,20,24,32,42):
+    for scale in ('none','sqrt','lam'):
+        f=lambda w: np.ones_like(w) if scale=='none' else (np.sqrt(np.abs(w)) if scale=='sqrt' else np.abs(w))
+        XV=UV[:,:r]*f(wV[:r]); XT=UT[:,:r]*f(wT[:r])
+        # oracle Procrustes using truth
+        M=XT.T@XV; u,s,vt=np.linalg.svd(M); O=u@vt
+        p=hung(XV,XT@O)
+        # also cosine-normalised rows
+        nV=XV/np.linalg.norm(XV,axis=1,keepdims=True); nT=(XT@O); nT=nT/np.linalg.norm(nT,axis=1,keepdims=True)
+        p2=hung(nV,nT)
+        print(f"  r={r:3d} scale={scale:5s} oracle-O acc={acc(p):.3f}  rownorm acc={acc(p2):.3f}")
diff --git a/artifacts/spectral_frontier_probe/diag6.py b/artifacts/spectral_frontier_probe/diag6.py
new file mode 100644
index 0000000..0158f9d
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag6.py
@@ -0,0 +1,46 @@
+import sys, time, torch, numpy as np
+sys.path.insert(0,'/home/yurenh2/emm')
+from scipy.optimize import linear_sum_assignment
+from scipy.stats import ortho_group
+from worldalign.synth_fast_gate import fast_pair_descent
+dev='cuda:3'
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(T,dtype=torch.float32,device=dev)
+CONST=float((Tt*Tt).sum()+(Vt*Vt).sum())
+def energy(p):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1))
+def descend(p,steps=4000):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy()
+truth=np.arange(N); acc=lambda p: float((p==truth).mean())
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1]
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+rng=np.random.default_rng(1)
+
+def icp(XV,XT,O,iters=40):
+    for _ in range(iters):
+        p=hung(XV,XT@O)
+        u,s,vt=np.linalg.svd(XT[p].T@XV); Onew=u@vt
+        if np.allclose(Onew,O,atol=1e-10): O=Onew; break
+        O=Onew
+    return hung(XV,XT@O),O
+
+print("### BLIND: random-restart ICP over O(r), scored by QAP energy")
+for r in (6,8,10,12,16):
+    XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r]))
+    t0=time.time(); best=(1e9,None)
+    R=200
+    for t in range(R):
+        O=ortho_group.rvs(r,random_state=int(rng.integers(1<<30)))
+        p,_=icp(XV,XT,O)
+        e=energy(p)
+        if e<best[0]: best=(e,p)
+    p=best[1]; pd=descend(p)
+    print(f"  r={r:3d} R={R}: best-E {best[0]:.4f} acc {acc(p):.3f} | after descent E {energy(pd):.4f} acc {acc(pd):.3f}  [{time.time()-t0:.0f}s]")
diff --git a/artifacts/spectral_frontier_probe/diag7.py b/artifacts/spectral_frontier_probe/diag7.py
new file mode 100644
index 0000000..c41b3c4
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/diag7.py
@@ -0,0 +1,74 @@
+import sys, time, torch, numpy as np
+sys.path.insert(0,'/home/yurenh2/emm')
+from scipy.optimize import linear_sum_assignment
+np.set_printoptions(precision=4, suppress=True, linewidth=200)
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy()
+V=standardise(V0); T=standardise(T0); N=len(V)
+truth=np.arange(N); acc=lambda p: float((p==truth).mean())
+CONST=float((T*T).sum()+(V*V).sum())
+def energy_from_align(S): return (CONST-2.0*S)/(N*(N-1))
+def energy(p): return energy_from_align(float((T[np.ix_(p,p)]*V).sum()))
+print("E(truth) =",energy(truth))
+
+print("\n### 2. SPECTRAL CERTIFICATES (lower bounds on E over all permutations)")
+lV=np.linalg.eigvalsh(V)[::-1]; lT=np.linalg.eigvalsh(T)[::-1]
+print("  Hoffman-Wielandt / von Neumann bound  E >= ", energy_from_align(float((lV*lT).sum())))
+# projected eigenvalue bound (Hadley-Rendl-Wolkowicz): project out the all-ones direction
+Q,_=np.linalg.qr(np.hstack([np.ones((N,1))/np.sqrt(N), np.random.default_rng(0).standard_normal((N,N-1))]))
+Vp=Q[:,1:].T@V@Q[:,1:]; Tp=Q[:,1:].T@T@Q[:,1:]
+lVp=np.linalg.eigvalsh(Vp)[::-1]; lTp=np.linalg.eigvalsh(Tp)[::-1]
+sV=V.sum(); sT=T.sum(); rV=V.sum(1); rT=T.sum(1)
+# <V,PTP'> = (1/N)*?; exact decomposition: with x=ones/sqrt(N),
+# tr(V P T P^T) = tr(Vp Pp Tp Pp^T) + 2 x^T V P T P^T x*? -- use the standard HRW split
+cross = float(rV.sum()*rT.sum())/ (N*N)  # x^T V x * x^T T x term
+# rank-1 cross terms bounded by sorted-inner-product of centred row sums
+cV=np.sort(rV-rV.mean())[::-1]; cT=np.sort(rT-rT.mean())[::-1]
+proj_bound = float((lVp*lTp).sum()) + cross + 2.0*float((cV*cT).sum())/N
+print("  projected (HRW-style, upper bd on alignment) E >= ", energy_from_align(proj_bound))
+print("  (both are LOWER bounds on E; E(truth)=%.4f, solvers report E(truth)+0.44=%.4f)"%(energy(truth),energy(truth)+0.44))
+
+print("\n### 3. BLIND rotation-invariant node descriptors -> Hungarian seeds")
+def report(name,DV,DT,topk=(12,25,50)):
+    DV=DV/ (np.linalg.norm(DV,axis=1,keepdims=True)+1e-12); DT=DT/(np.linalg.norm(DT,axis=1,keepdims=True)+1e-12)
+    C=((DV**2).sum(1)[:,None]+(DT**2).sum(1)[None,:]-2*DV@DT.T)
+    r,c=linear_sum_assignment(C)
+    a=acc(c)
+    # confidence = margin between assigned cost and 2nd best in row
+    Cm=C.copy(); Cm[np.arange(N),c]=np.inf
+    margin=Cm.min(1)-C[np.arange(N),c]
+    order=np.argsort(-margin)
+    prec={k: float((c[order[:k]]==order[:k]).mean()) for k in topk}
+    print(f"  {name:34s} full-acc {a:.3f} | precision@top-margin {prec}")
+    return c
+# (a) sorted row profile
+report("sorted row profile (V vs T)", np.sort(V,1), np.sort(T,1))
+# (b) row moments
+def mom(M,K=8): return np.stack([ (M**k).mean(1) for k in range(1,K+1)],1)
+report("row power moments k=1..8", mom(V), mom(T))
+# (c) heat kernel signature on normalised Laplacian of the raw non-negative Gram
+def lap_eigs(W):
+    W=np.clip(W,0,None).copy(); np.fill_diagonal(W,0.0)
+    dg=W.sum(1); Dm=1/np.sqrt(np.maximum(dg,1e-12))
+    L=np.eye(len(W))-(Dm[:,None]*W*Dm[None,:])
+    w,U=np.linalg.eigh(L); return w,U,dg
+wV,UVl,dV=lap_eigs(V0); wT,UTl,dT=lap_eigs(T0)
+ts=np.logspace(-2,1.5,24)
+HV=np.stack([ (np.exp(-t*wV)[None,:]*UVl**2).sum(1) for t in ts],1)
+HT=np.stack([ (np.exp(-t*wT)[None,:]*UTl**2).sum(1) for t in ts],1)
+report("HKS (normalised Laplacian)", np.log(HV+1e-12), np.log(HT+1e-12))
+# (d) wave kernel signature
+def wks(w,U,M=24):
+    lw=np.log(np.maximum(w,1e-6)); e=np.linspace(lw.min(),lw.max(),M); sig=(e[1]-e[0])*2
+    return np.stack([ (np.exp(-((e_-lw)**2)/(2*sig**2))[None,:]*U**2).sum(1) for e_ in e],1)
+report("WKS (normalised Laplacian)", wks(wV,UVl), wks(wT,UTl))
+# (e) spectral graph wavelet (SGWT) coefficient energies
+def sgwt(w,U,S=10):
+    ss=np.logspace(-1.5,1.0,S)
+    return np.stack([ ((s*w*np.exp(-s*w))[None,:]*U**2).sum(1) for s in ss],1)
+report("SGWT band energies", sgwt(wV,UVl), sgwt(wT,UTl))
+# (f) degree only
+report("degree (row sum) only", dV[:,None], dT[:,None])
diff --git a/artifacts/spectral_frontier_probe/seed.log b/artifacts/spectral_frontier_probe/seed.log
new file mode 100644
index 0000000..87f0744
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/seed.log
@@ -0,0 +1,7 @@
+360 ICP solutions in 760s; corr(E,acc)=-0.520
+lowest-10 (E,acc): [(0.5186, 0.461), (0.5815, 0.27), (0.622, 0.258), (0.6385, 0.227), (0.6637, 0.051), (0.6652, 0.059), (0.6828, 0.109), (0.6903, 0.074), (0.6977, 0.016), (0.6988, 0.016)]
+E quantiles [0.519 0.759 0.835 0.91  1.031]  acc of best-E: 0.4609375
+  consensus over 5 lowest-E sols: precision@12 = 0.833  consensus over 5 lowest-E sols: precision@25 = 0.840  consensus over 5 lowest-E sols: precision@50 = 0.720  consensus over 5 lowest-E sols: precision@100 = 0.680
+  consensus over 10 lowest-E sols: precision@12 = 0.750  consensus over 10 lowest-E sols: precision@25 = 0.720  consensus over 10 lowest-E sols: precision@50 = 0.760  consensus over 10 lowest-E sols: precision@100 = 0.630
+  consensus over 20 lowest-E sols: precision@12 = 0.917  consensus over 20 lowest-E sols: precision@25 = 0.840  consensus over 20 lowest-E sols: precision@50 = 0.720  consensus over 20 lowest-E sols: precision@100 = 0.590
+  consensus over 40 lowest-E sols: precision@12 = 0.917  consensus over 40 lowest-E sols: precision@25 = 0.880  consensus over 40 lowest-E sols: precision@50 = 0.740  consensus over 40 lowest-E sols: precision@100 = 0.570
diff --git a/artifacts/spectral_frontier_probe/seed.py b/artifacts/spectral_frontier_probe/seed.py
new file mode 100644
index 0000000..5e33d7b
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/seed.py
@@ -0,0 +1,49 @@
+import sys, time, torch, numpy as np
+sys.path.insert(0,'/home/yurenh2/emm')
+from scipy.optimize import linear_sum_assignment
+from scipy.stats import ortho_group
+dev='cuda:0'
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+rng=np.random.default_rng(31337); sigma=rng.permutation(N)
+Ts=T[np.ix_(sigma,sigma)]; truth=np.argsort(sigma)
+Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(Ts,dtype=torch.float32,device=dev)
+CONST=float((Tt*Tt).sum()+(Vt*Vt).sum())
+def energy(p):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1))
+acc=lambda p: float((p==truth).mean())
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(Ts); wT=wT[::-1]; UT=UT[:,::-1]
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+def icp(XV,XT,O,iters=25):
+    for _ in range(iters):
+        p=hung(XV,XT@O); u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt
+        if np.allclose(On,O,atol=1e-10): O=On; break
+        O=On
+    return hung(XV,XT@O)
+sols=[]
+t0=time.time()
+for r in (10,12,14,16,18,20):
+    XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r]))
+    for t in range(60):
+        p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30))))
+        sols.append((energy(p),acc(p),p,r))
+sols.sort(key=lambda z:z[0])
+E=np.array([s[0] for s in sols]); A=np.array([s[1] for s in sols])
+print(f"{len(sols)} ICP solutions in {time.time()-t0:.0f}s; corr(E,acc)={np.corrcoef(E,A)[0,1]:.3f}")
+print("lowest-10 (E,acc):", [(round(e,4),round(a,3)) for e,a,_,_ in sols[:10]])
+print("E quantiles",np.quantile(E,[0,.1,.5,.9,1]).round(3)," acc of best-E:",A[0])
+for m in (5,10,20,40):
+    votes=np.zeros((N,N))
+    for e,a,p,r in sols[:m]: votes[np.arange(N),p]+=1
+    conf=votes.max(1); pick=np.argsort(-conf)
+    pm=votes.argmax(1)
+    for k in (12,25,50,100):
+        sel=pick[:k]; prec=float((pm[sel]==truth[sel]).mean())
+        print(f"  consensus over {m} lowest-E sols: precision@{k} = {prec:.3f}", end='')
+    print()
diff --git a/artifacts/spectral_frontier_probe/shufctl.log b/artifacts/spectral_frontier_probe/shufctl.log
new file mode 100644
index 0000000..8961e1c
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/shufctl.log
@@ -0,0 +1,4 @@
+E(truth)= 0.3396347943474265
+SHUFFLED r=16: pre E 0.3505 acc 0.793 -> post E 0.3396 acc 0.820 [361s]
+SHUFFLED r=6: pre E 0.5712 acc 0.078 -> post E 0.5444 acc 0.109 [195s]
+SHUFFLED r=12: pre E 0.5067 acc 0.402 -> post E 0.3540 acc 0.809 [410s]
diff --git a/artifacts/spectral_frontier_probe/shufctl.py b/artifacts/spectral_frontier_probe/shufctl.py
new file mode 100644
index 0000000..5913f53
--- /dev/null
+++ b/artifacts/spectral_frontier_probe/shufctl.py
@@ -0,0 +1,47 @@
+import sys, time, torch, numpy as np
+sys.path.insert(0,'/home/yurenh2/emm')
+from scipy.optimize import linear_sum_assignment
+from scipy.stats import ortho_group
+from worldalign.synth_fast_gate import fast_pair_descent
+dev='cuda:1'
+def standardise(M):
+    M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask]
+    out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out
+d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu')
+V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V)
+rng=np.random.default_rng(777)
+sigma=rng.permutation(N)                 # hidden shuffle applied to T
+Ts=T[np.ix_(sigma,sigma)]                # scene i of V corresponds to row where sigma[k]=i -> inverse
+inv=np.argsort(sigma)                    # truth: V index i <-> Ts index inv[i]
+truth=inv
+Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(Ts,dtype=torch.float32,device=dev)
+CONST=float((Tt*Tt).sum()+(Vt*Vt).sum())
+def energy(p):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1))
+def descend(p,steps=4000):
+    P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev)
+    return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy()
+acc=lambda p: float((p==truth).mean())
+print("E(truth)=",energy(truth),flush=True)
+wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1]
+wT,UT=np.linalg.eigh(Ts); wT=wT[::-1]; UT=UT[:,::-1]
+def hung(A,B):
+    C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c
+def icp(XV,XT,O,iters=30):
+    for _ in range(iters):
+        p=hung(XV,XT@O); u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt
+        if np.allclose(On,O,atol=1e-10): O=On; break
+        O=On
+    return hung(XV,XT@O)
+for r in (16,6,12):
+    XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r]))
+    t0=time.time(); cand=[]
+    for t in range(200):
+        p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30))))
+        cand.append((energy(p),p))
+    cand.sort(key=lambda z:z[0]); best=None
+    for e,p in cand[:5]:
+        pd=descend(p); ed=energy(pd)
+        if best is None or ed<best[0]: best=(ed,pd)
+    print(f"SHUFFLED r={r}: pre E {cand[0][0]:.4f} acc {acc(cand[0][1]):.3f} -> post E {best[0]:.4f} acc {acc(best[1]):.3f} [{time.time()-t0:.0f}s]",flush=True)