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>

1 files changed, 74 insertions, 0 deletions

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])