From ef104fe4f07713bf11f266d2954b7446f176f8ae Mon Sep 17 00:00:00 2001 From: YurenHao0426 Date: Sat, 1 Aug 2026 19:32:58 -0500 Subject: 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 --- artifacts/spectral_frontier_probe/band.py | 26 ++++++++++ artifacts/spectral_frontier_probe/confirm.log | 7 +++ artifacts/spectral_frontier_probe/confirm.py | 49 ++++++++++++++++++ artifacts/spectral_frontier_probe/diag10.py | 41 +++++++++++++++ artifacts/spectral_frontier_probe/diag11.log | 8 +++ artifacts/spectral_frontier_probe/diag11.py | 52 +++++++++++++++++++ artifacts/spectral_frontier_probe/diag3.py | 21 ++++++++ artifacts/spectral_frontier_probe/diag4.py | 29 +++++++++++ artifacts/spectral_frontier_probe/diag6.py | 46 +++++++++++++++++ artifacts/spectral_frontier_probe/diag7.py | 74 +++++++++++++++++++++++++++ artifacts/spectral_frontier_probe/seed.log | 7 +++ artifacts/spectral_frontier_probe/seed.py | 49 ++++++++++++++++++ artifacts/spectral_frontier_probe/shufctl.log | 4 ++ artifacts/spectral_frontier_probe/shufctl.py | 47 +++++++++++++++++ 14 files changed, 460 insertions(+) create mode 100644 artifacts/spectral_frontier_probe/band.py create mode 100644 artifacts/spectral_frontier_probe/confirm.log create mode 100644 artifacts/spectral_frontier_probe/confirm.py create mode 100644 artifacts/spectral_frontier_probe/diag10.py create mode 100644 artifacts/spectral_frontier_probe/diag11.log create mode 100644 artifacts/spectral_frontier_probe/diag11.py create mode 100644 artifacts/spectral_frontier_probe/diag3.py create mode 100644 artifacts/spectral_frontier_probe/diag4.py create mode 100644 artifacts/spectral_frontier_probe/diag6.py create mode 100644 artifacts/spectral_frontier_probe/diag7.py create mode 100644 artifacts/spectral_frontier_probe/seed.log create mode 100644 artifacts/spectral_frontier_probe/seed.py create mode 100644 artifacts/spectral_frontier_probe/shufctl.log create mode 100644 artifacts/spectral_frontier_probe/shufctl.py (limited to 'artifacts/spectral_frontier_probe') 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 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("|| top-16 block (rows=V index, cols=T index):") +print(X[:16,:16]) +print("\ndiagonal || 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= ", 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) +# = (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 post E {best[0]:.4f} acc {acc(best[1]):.3f} [{time.time()-t0:.0f}s]",flush=True) -- cgit v1.2.3