diff options
Diffstat (limited to 'artifacts/spectral_frontier_probe')
| -rw-r--r-- | artifacts/spectral_frontier_probe/band.py | 26 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/confirm.log | 7 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/confirm.py | 49 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag10.py | 41 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag11.log | 8 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag11.py | 52 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag3.py | 21 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag4.py | 29 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag6.py | 46 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/diag7.py | 74 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/seed.log | 7 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/seed.py | 49 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/shufctl.log | 4 | ||||
| -rw-r--r-- | artifacts/spectral_frontier_probe/shufctl.py | 47 |
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) |