diff options
Diffstat (limited to 'artifacts/spectral_frontier_probe')
10 files changed, 551 insertions, 0 deletions
diff --git a/artifacts/spectral_frontier_probe/eval_impure_anchors.log b/artifacts/spectral_frontier_probe/eval_impure_anchors.log new file mode 100644 index 0000000..a4b957a --- /dev/null +++ b/artifacts/spectral_frontier_probe/eval_impure_anchors.log @@ -0,0 +1,60 @@ +=== (a) PURE random anchors: reproduce the published curve === + K= 12 pure accuracy 0.589 +- 0.055 + K= 25 pure accuracy 0.860 +- 0.058 + K= 40 pure accuracy 0.926 +- 0.032 + K= 50 pure accuracy 0.949 +- 0.026 + K= 64 pure accuracy 0.971 +- 0.017 + K= 128 pure accuracy 0.996 +- 0.005 + +=== (b) IMPURE random anchors: c correct, K-c wrong (wrong = derangement among selected) === + K= 12 prec=1.00 (12/12 right) accuracy 0.620 +- 0.093 + K= 12 prec=0.92 (10/12 right) accuracy 0.390 +- 0.101 + K= 12 prec=0.88 (10/12 right) accuracy 0.399 +- 0.109 + K= 12 prec=0.84 (10/12 right) accuracy 0.330 +- 0.098 + K= 12 prec=0.80 (10/12 right) accuracy 0.408 +- 0.120 + K= 12 prec=0.72 (9/12 right) accuracy 0.264 +- 0.068 + K= 12 prec=0.60 (7/12 right) accuracy 0.135 +- 0.064 + + K= 25 prec=1.00 (25/25 right) accuracy 0.852 +- 0.049 + K= 25 prec=0.92 (23/25 right) accuracy 0.766 +- 0.073 + K= 25 prec=0.88 (22/25 right) accuracy 0.746 +- 0.071 + K= 25 prec=0.84 (21/25 right) accuracy 0.704 +- 0.079 + K= 25 prec=0.80 (20/25 right) accuracy 0.677 +- 0.070 + K= 25 prec=0.72 (18/25 right) accuracy 0.527 +- 0.103 + K= 25 prec=0.60 (15/25 right) accuracy 0.371 +- 0.091 + + K= 40 prec=1.00 (40/40 right) accuracy 0.933 +- 0.032 + K= 40 prec=0.92 (37/40 right) accuracy 0.902 +- 0.038 + K= 40 prec=0.88 (35/40 right) accuracy 0.859 +- 0.037 + K= 40 prec=0.84 (34/40 right) accuracy 0.835 +- 0.038 + K= 40 prec=0.80 (32/40 right) accuracy 0.812 +- 0.047 + K= 40 prec=0.72 (29/40 right) accuracy 0.742 +- 0.085 + K= 40 prec=0.60 (24/40 right) accuracy 0.574 +- 0.081 + + K= 50 prec=1.00 (50/50 right) accuracy 0.950 +- 0.024 + K= 50 prec=0.92 (46/50 right) accuracy 0.920 +- 0.026 + K= 50 prec=0.88 (44/50 right) accuracy 0.893 +- 0.033 + K= 50 prec=0.84 (42/50 right) accuracy 0.877 +- 0.033 + K= 50 prec=0.80 (40/50 right) accuracy 0.855 +- 0.031 + K= 50 prec=0.72 (36/50 right) accuracy 0.783 +- 0.055 + K= 50 prec=0.60 (30/50 right) accuracy 0.629 +- 0.102 + + K= 100 prec=1.00 (100/100 right) accuracy 0.992 +- 0.009 + K= 100 prec=0.92 (92/100 right) accuracy 0.953 +- 0.014 + K= 100 prec=0.88 (88/100 right) accuracy 0.933 +- 0.016 + K= 100 prec=0.84 (84/100 right) accuracy 0.915 +- 0.018 + K= 100 prec=0.80 (80/100 right) accuracy 0.895 +- 0.012 + K= 100 prec=0.72 (72/100 right) accuracy 0.852 +- 0.020 + K= 100 prec=0.60 (60/100 right) accuracy 0.787 +- 0.021 + +=== (c) IMPURE with wrong partners drawn from OUTSIDE the anchor set === + K= 25 prec=1.00 accuracy 0.866 +- 0.040 + K= 25 prec=0.88 accuracy 0.729 +- 0.064 + K= 25 prec=0.84 accuracy 0.689 +- 0.070 + K= 25 prec=0.80 accuracy 0.603 +- 0.079 + K= 25 prec=0.72 accuracy 0.501 +- 0.094 + K= 50 prec=1.00 accuracy 0.945 +- 0.025 + K= 50 prec=0.88 accuracy 0.832 +- 0.031 + K= 50 prec=0.84 accuracy 0.808 +- 0.020 + K= 50 prec=0.80 accuracy 0.775 +- 0.031 + K= 50 prec=0.72 accuracy 0.689 +- 0.043 diff --git a/artifacts/spectral_frontier_probe/eval_impure_anchors.py b/artifacts/spectral_frontier_probe/eval_impure_anchors.py new file mode 100644 index 0000000..f6b4397 --- /dev/null +++ b/artifacts/spectral_frontier_probe/eval_impure_anchors.py @@ -0,0 +1,65 @@ +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', weights_only=False) +V = standardise(d['visual_field']); T = standardise(d['text_field']); N = len(V) +truth = np.arange(N) + +def expand(anchorL, anchorR): + """anchorL: vision indices; anchorR: claimed text partner. Returns full accuracy.""" + probe = np.setdiff1d(np.arange(N), anchorL) + # probe scenes on the text side are whatever is not claimed + probeR = np.setdiff1d(np.arange(N), anchorR) + L = V[np.ix_(probe, anchorL)]; R = T[np.ix_(probeR, anchorR)] + L = (L - L.mean(1, keepdims=True)) / L.std(1, keepdims=True).clip(1e-9) + R = (R - R.mean(1, keepdims=True)) / R.std(1, keepdims=True).clip(1e-9) + S = L @ R.T / L.shape[1] + _, cols = linear_sum_assignment(-S) + resolved = probeR[cols] + n_ok = int((resolved == probe).sum()) + int((anchorR == anchorL).sum()) + return n_ok / N + +rng = np.random.default_rng(0) +REP = 20 +print("=== (a) PURE random anchors: reproduce the published curve ===") +for K in (12, 25, 40, 50, 64, 128): + a = [expand(*(lambda s: (s, s))(np.sort(rng.choice(N, K, replace=False)))) for _ in range(REP)] + print(f" K={K:4d} pure accuracy {np.mean(a):.3f} +- {np.std(a):.3f}") + +print("\n=== (b) IMPURE random anchors: c correct, K-c wrong (wrong = derangement among selected) ===") +for K in (12, 25, 40, 50, 100): + for prec in (1.0, 0.92, 0.88, 0.84, 0.80, 0.72, 0.60): + nb = int(round(K * (1 - prec))) + if nb == 1: nb = 2 # a single wrong pair is impossible inside a bijection + accs = [] + for _ in range(REP): + sel = np.sort(rng.choice(N, K, replace=False)) + right = sel.copy() + if nb >= 2: + bad = rng.choice(K, nb, replace=False) + sh = right[bad].copy() + while True: + perm = rng.permutation(nb) + if not (perm == np.arange(nb)).any(): break + right[bad] = sh[perm] + accs.append(expand(sel, right)) + print(f" K={K:4d} prec={prec:.2f} ({K-nb}/{K} right) accuracy {np.mean(accs):.3f} +- {np.std(accs):.3f}") + print() + +print("=== (c) IMPURE with wrong partners drawn from OUTSIDE the anchor set ===") +for K in (25, 50): + for prec in (1.0, 0.88, 0.84, 0.80, 0.72): + nb = int(round(K * (1 - prec))); accs = [] + for _ in range(REP): + sel = np.sort(rng.choice(N, K, replace=False)) + right = sel.copy() + if nb: + bad = rng.choice(K, nb, replace=False) + outside = np.setdiff1d(np.arange(N), sel) + right[bad] = rng.choice(outside, nb, replace=False) + accs.append(expand(sel, right)) + print(f" K={K:4d} prec={prec:.2f} accuracy {np.mean(accs):.3f} +- {np.std(accs):.3f}") diff --git a/artifacts/spectral_frontier_probe/eval_voteanchor_seed7.log b/artifacts/spectral_frontier_probe/eval_voteanchor_seed7.log new file mode 100644 index 0000000..8852244 --- /dev/null +++ b/artifacts/spectral_frontier_probe/eval_voteanchor_seed7.log @@ -0,0 +1,23 @@ +150 ICP solutions in 222s +lowest-8 (E, acc, r): [(0.6733, 0.008, 8), (0.6824, 0.027, 10), (0.6956, 0.031, 10), (0.6988, 0.016, 8), (0.7011, 0.074, 12), (0.7026, 0.023, 10), (0.7052, 0.004, 10), (0.7052, 0.012, 8)] +mean pairwise agreement among top-30: 0.046 + +--- ROUTE A (free readout of the same ensemble) --- + best-E member acc 0.008 -> after exact descent acc 0.012 E=0.6565 [5s] + descend 5 lowest-E, pick by energy: acc 0.043 E=0.6316 + +--- ROUTE B (the proposal) --- + m= 5 K= 12 kept= 12 precision=0.167 -> expansion acc 0.039 -> after descent 0.039 + m= 5 K= 25 kept= 25 precision=0.080 -> expansion acc 0.031 -> after descent 0.051 + m= 5 K= 50 kept= 47 precision=0.064 -> expansion acc 0.031 -> after descent 0.023 + m= 5 K=100 kept= 91 precision=0.044 -> expansion acc 0.031 -> after descent 0.023 + m=10 K= 12 kept= 12 precision=0.250 -> expansion acc 0.043 -> after descent 0.074 + m=10 K= 25 kept= 23 precision=0.130 -> expansion acc 0.035 -> after descent 0.043 + m=10 K= 50 kept= 48 precision=0.104 -> expansion acc 0.059 -> after descent 0.066 + m=10 K=100 kept= 87 precision=0.069 -> expansion acc 0.035 -> after descent 0.035 + m=30 K= 12 kept= 12 precision=0.417 -> expansion acc 0.094 -> after descent 0.180 + m=30 K= 25 kept= 25 precision=0.320 -> expansion acc 0.141 -> after descent 0.125 + m=30 K= 50 kept= 48 precision=0.229 -> expansion acc 0.215 -> after descent 0.816 + m=30 K=100 kept= 80 precision=0.188 -> expansion acc 0.148 -> after descent 0.195 + +E(truth)=0.3396 total wall 617s diff --git a/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.log b/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.log new file mode 100644 index 0000000..c2bd82f --- /dev/null +++ b/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.log @@ -0,0 +1,34 @@ + +##### hidden seed 11: 150 ICP sols in 200s, E(truth)=0.3396, best member E=0.6569 acc=0.172 + ROUTE A (descend 5 lowest-E, blind pick): E=0.5863 acc=0.277 + m=10 K= 12 prec=0.333 expand=0.102 -> descent acc=0.320 E=0.6184 + m=10 K= 25 prec=0.333 expand=0.168 -> descent acc=0.324 E=0.6279 + m=10 K= 50 prec=0.255 expand=0.184 -> descent acc=0.793 E=0.3544 + m=10 K=100 prec=0.193 expand=0.172 -> descent acc=0.770 E=0.3900 + m=30 K= 12 prec=0.455 expand=0.211 -> descent acc=0.535 E=0.4867 + m=30 K= 25 prec=0.292 expand=0.195 -> descent acc=0.754 E=0.4163 + m=30 K= 50 prec=0.227 expand=0.188 -> descent acc=0.695 E=0.4357 + m=30 K=100 prec=0.172 expand=0.184 -> descent acc=0.773 E=0.3837 + m=60 K= 12 prec=0.455 expand=0.086 -> descent acc=0.094 E=0.7738 + m=60 K= 25 prec=0.375 expand=0.164 -> descent acc=0.133 E=0.6840 + m=60 K= 50 prec=0.244 expand=0.098 -> descent acc=0.066 E=0.7911 + m=60 K=100 prec=0.171 expand=0.113 -> descent acc=0.086 E=0.7412 + ROUTE B (blind pick by energy over the 12 cells): m=10 K=50 E=0.3544 acc=0.793 [best cell by accuracy was 0.793] + wall 492s + +##### hidden seed 23: 150 ICP sols in 559s, E(truth)=0.3396, best member E=0.3684 acc=0.691 + ROUTE A (descend 5 lowest-E, blind pick): E=0.3396 acc=0.852 + m=10 K= 12 prec=0.917 expand=0.703 -> descent acc=0.855 E=0.3396 + m=10 K= 25 prec=0.680 expand=0.520 -> descent acc=0.758 E=0.3792 + m=10 K= 50 prec=0.583 expand=0.484 -> descent acc=0.777 E=0.3543 + m=10 K=100 prec=0.505 expand=0.492 -> descent acc=0.848 E=0.3396 + m=30 K= 12 prec=0.917 expand=0.656 -> descent acc=0.812 E=0.3466 + m=30 K= 25 prec=0.920 expand=0.723 -> descent acc=0.809 E=0.3396 + m=30 K= 50 prec=0.673 expand=0.691 -> descent acc=0.816 E=0.3644 + m=30 K=100 prec=0.558 expand=0.547 -> descent acc=0.797 E=0.3593 + m=60 K= 12 prec=0.833 expand=0.504 -> descent acc=0.785 E=0.3468 + m=60 K= 25 prec=0.792 expand=0.734 -> descent acc=0.820 E=0.3396 + m=60 K= 50 prec=0.667 expand=0.641 -> descent acc=0.836 E=0.3433 + m=60 K=100 prec=0.524 expand=0.496 -> descent acc=0.844 E=0.3569 + ROUTE B (blind pick by energy over the 12 cells): m=30 K=25 E=0.3396 acc=0.809 [best cell by accuracy was 0.855] + wall 752s diff --git a/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.py b/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.py new file mode 100644 index 0000000..d97d239 --- /dev/null +++ b/artifacts/spectral_frontier_probe/eval_voteanchor_vs_bestE.py @@ -0,0 +1,78 @@ +"""Replication with BLIND energy selection over the (m,K) grid, two fresh hidden perms.""" +import sys, time, numpy as np, torch +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:0' +def standardise(M): + M = np.asarray(M, dtype=np.float64); m = ~np.eye(len(M), dtype=bool); v = M[m] + o = (M - v.mean())/v.std(); np.fill_diagonal(o, 0.0); return o +d = torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt', map_location='cpu', weights_only=False) +V = standardise(d['visual_field']); T0 = standardise(d['text_field']); N = len(V) +Vt = torch.tensor(V, dtype=torch.float32, device=dev) +wV, UV = np.linalg.eigh(V); wV = wV[::-1]; UV = UV[:, ::-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) + +for SEED in (11, 23): + rng = np.random.default_rng(SEED); sigma = rng.permutation(N) + T = T0[np.ix_(sigma, sigma)]; truth = np.argsort(sigma) + acc = lambda p: float((np.asarray(p) == truth).mean()) + 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.ascontiguousarray(p), dtype=torch.long, device=dev) + return fast_pair_descent(Tt, Vt, P, steps).cpu().numpy() + wT, UT = np.linalg.eigh(T); wT = wT[::-1]; UT = UT[:, ::-1] + t0 = time.time(); sols = [] + for r in (8, 10, 12, 14, 16): + XV = UV[:, :r]*np.sqrt(np.abs(wV[:r])); XT = UT[:, :r]*np.sqrt(np.abs(wT[:r])) + for _ in range(30): + p = icp(XV, XT, ortho_group.rvs(r, random_state=int(rng.integers(1 << 30)))) + sols.append((energy(p), p, r)) + sols.sort(key=lambda z: z[0]) + print(f"\n##### hidden seed {SEED}: {len(sols)} ICP sols in {time.time()-t0:.0f}s, " + f"E(truth)={energy(truth):.4f}, best member E={sols[0][0]:.4f} acc={acc(sols[0][1]):.3f}", flush=True) + A = [descend(p) for _, p, _ in sols[:5]] + bA = min(A, key=energy) + print(f" ROUTE A (descend 5 lowest-E, blind pick): E={energy(bA):.4f} acc={acc(bA):.3f}", flush=True) + + def expand(aL, aR): + pr = np.setdiff1d(np.arange(N), aL); prR = np.setdiff1d(np.arange(N), aR) + L = V[np.ix_(pr, aL)]; R = T[np.ix_(prR, aR)] + L = (L-L.mean(1, keepdims=True))/L.std(1, keepdims=True).clip(1e-9) + R = (R-R.mean(1, keepdims=True))/R.std(1, keepdims=True).clip(1e-9) + S = L@R.T/L.shape[1]; _, c = linear_sum_assignment(-S) + full = np.zeros(N, dtype=int); full[aL] = aR; full[pr] = prR[c]; return full + grid = [] + for m in (10, 30, 60): + votes = np.zeros((N, N)) + for _, p, _ in sols[:m]: votes[np.arange(N), p] += 1 + conf = votes.max(1) - np.partition(votes, -2, axis=1)[:, -2] + order = np.argsort(-conf); pm = votes.argmax(1) + for K in (12, 25, 50, 100): + sel = order[:K]; part = pm[sel]; keep = np.zeros(len(sel), bool); seen = set() + for i, c in enumerate(part): + if c not in seen: seen.add(c); keep[i] = True + sel, part = sel[keep], part[keep] + prec = float((part == truth[sel]).mean()) + e = expand(sel, part); pd = descend(e) + grid.append((energy(pd), acc(pd), m, K, prec, acc(e))) + print(f" m={m:2d} K={K:3d} prec={prec:.3f} expand={acc(e):.3f} -> descent acc={acc(pd):.3f} E={energy(pd):.4f}", flush=True) + grid.sort() + print(f" ROUTE B (blind pick by energy over the 12 cells): m={grid[0][2]} K={grid[0][3]} " + f"E={grid[0][0]:.4f} acc={grid[0][1]:.3f} [best cell by accuracy was {max(g[1] for g in grid):.3f}]", flush=True) + print(f" wall {time.time()-t0:.0f}s", flush=True) diff --git a/artifacts/spectral_frontier_probe/nmf_adv_oracle.log b/artifacts/spectral_frontier_probe/nmf_adv_oracle.log new file mode 100644 index 0000000..8176ba3 --- /dev/null +++ b/artifacts/spectral_frontier_probe/nmf_adv_oracle.log @@ -0,0 +1,26 @@ +E(truth) = 0.339635 +REF: GRAMPA on raw fields -> raw 0.0039 refined 0.0117 E=0.8380 +REF: field correlation (aligned) = 0.8302 + +=== r=24 [6.2s] resid V 0.0400 T 0.0331 V(seed2) 0.0447 + [uniqueness] symNMF(V,seed1) vs symNMF(V,seed2) matched cos: median 0.885 >0.9 10/24 >0.5 24/24 min 0.541 + [oracle cols] cos: >0.9 10/24 >0.5 23/24 median 0.842 min 0.427 + sorted: [0.98 0.98 0.97 0.96 0.95 0.94 0.94 0.92 0.91 0.91 0.88 0.86 0.82 0.8 + 0.79 0.69 0.68 0.68 0.68 0.68 0.68 0.63 0.55 0.43] + identity-map fraction (shared-seed artifact check): 0.083 + [oracle-col row descriptors] Hungarian acc 0.6719 -> descent 0.8320 + [prune thr=0.5] kept 23/24 rho(cleanV,T)=0.8466 (was 0.8302) GRAMPA 0.0078 -> clean-descent 0.0547 -> true-descent 0.0547 E=0.6471 + [prune thr=0.7] kept 15/24 rho(cleanV,T)=0.7684 (was 0.8302) GRAMPA 0.0039 -> clean-descent 0.0391 -> true-descent 0.0391 E=0.6817 + [prune thr=0.8] kept 13/24 rho(cleanV,T)=0.7804 (was 0.8302) GRAMPA 0.0000 -> clean-descent 0.0000 -> true-descent 0.0000 E=0.7292 + +=== r=42 [32.2s] resid V 0.0241 T 0.0141 V(seed2) 0.0231 + [uniqueness] symNMF(V,seed1) vs symNMF(V,seed2) matched cos: median 0.832 >0.9 12/42 >0.5 42/42 min 0.541 + [oracle cols] cos: >0.9 4/42 >0.5 41/42 median 0.766 min 0.416 + sorted: [0.97 0.94 0.93 0.9 0.89 0.88 0.88 0.87 0.86 0.86 0.86 0.84 0.84 0.84 + 0.82 0.81 0.81 0.8 0.79 0.78 0.77 0.76 0.76 0.76 0.75 0.75 0.75 0.74 + 0.74 0.73 0.72 0.72 0.72 0.71 0.7 0.7 0.69 0.66 0.65 0.63 0.62 0.42] + identity-map fraction (shared-seed artifact check): 0.000 + [oracle-col row descriptors] Hungarian acc 0.8125 -> descent 0.8828 + [prune thr=0.5] kept 41/42 rho(cleanV,T)=0.8337 (was 0.8302) GRAMPA 0.0273 -> clean-descent 0.0547 -> true-descent 0.0430 E=0.6804 + [prune thr=0.7] kept 35/42 rho(cleanV,T)=0.7976 (was 0.8302) GRAMPA 0.0391 -> clean-descent 0.0625 -> true-descent 0.0859 E=0.6966 + [prune thr=0.8] kept 17/42 rho(cleanV,T)=0.6036 (was 0.8302) GRAMPA 0.0000 -> clean-descent 0.0039 -> true-descent 0.0078 E=0.6982 diff --git a/artifacts/spectral_frontier_probe/nmf_adv_oracle.py b/artifacts/spectral_frontier_probe/nmf_adv_oracle.py new file mode 100644 index 0000000..2702e29 --- /dev/null +++ b/artifacts/spectral_frontier_probe/nmf_adv_oracle.py @@ -0,0 +1,133 @@ +"""Decisive ORACLE test of the symNMF de-rotation / component-pruning candidate. + +If the ORACLE version (perfect column match, perfect pruning) does not beat the +current blind state of the art (~0.83 on omit-size), the candidate is dead. +""" +import numpy as np, torch, time, warnings +warnings.filterwarnings('ignore') +from scipy.optimize import linear_sum_assignment + +dev = torch.device('cuda:0') +DT = torch.float64 + +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 standardise(M): + M = np.asarray(M, float); mask = ~np.eye(len(M), dtype=bool); v = M[mask] + out = (M - v.mean()) / v.std(); np.fill_diagonal(out, 0.0); return out + + +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_gpu(M, r, iters=3000, seed=0): + Mg = torch.tensor(np.clip(M, 0, None), device=dev, dtype=DT) + g = torch.Generator(device='cpu').manual_seed(seed) + W = torch.abs(torch.randn(len(M), r, generator=g, dtype=DT)).to(dev) * float(np.sqrt(M.mean() / r)) + for _ in range(iters): + W = W * (0.5 + 0.5 * (Mg @ W) / (W @ (W.T @ W) + 1e-12)) + return W.cpu().numpy() + + +# ---- energy + exact steepest descent (closed-form swap table) ---- +Vs = standardise(V0); Ts = standardise(T0s) +CONST = float((Vs * Vs).sum() + (Ts * Ts).sum()); DEN = N * (N - 1) +Ag = torch.tensor(Vs, device=dev, dtype=DT); Tg = torch.tensor(Ts, device=dev, dtype=DT) + + +def energy(p): + return (CONST - 2.0 * float((Ts[np.ix_(p, p)] * Vs).sum())) / DEN + + +def descend(p0, max_steps=4000, A=None, B=None): + A = Ag if A is None else A + Bfull = Tg if B is None else B + p = torch.tensor(np.asarray(p0), device=dev, dtype=torch.long) + iu = torch.triu_indices(N, N, offset=1, device=dev) + for _ in range(max_steps): + B = Bfull[p][:, p] + C = A @ B + dg = torch.diagonal(C) + G = C + C.T - dg[:, None] - dg[None, :] + 2 * A * B + vals = G[iu[0], iu[1]] + k = int(vals.argmax()) + if float(vals[k]) <= 1e-12: + break + u, v = int(iu[0][k]), int(iu[1][k]) + p[u], p[v] = p[v].clone(), p[u].clone() + return p.cpu().numpy() + + +def grampa(Vm, Tm, eta=0.2): + lam, u = np.linalg.eigh(Vm); mu, v = np.linalg.eigh(Tm) + ones = np.ones(len(Vm)); left = u.T @ ones; right = v.T @ ones + w = np.outer(left, right) / ((lam[:, None] - mu[None, :]) ** 2 + eta**2) + s = u @ w @ v.T + r, c = linear_sum_assignment(-s); return c + + +print(f"E(truth) = {energy(truth):.6f}") +print("REF: GRAMPA on raw fields ->", end=" ") +p = grampa(Vs, Ts); pd = descend(p) +print(f"raw {acc(p):.4f} refined {acc(pd):.4f} E={energy(pd):.4f}") +rho_raw = np.corrcoef(Vs[~np.eye(N, dtype=bool)], Ts[np.ix_(truth, truth)][~np.eye(N, dtype=bool)])[0, 1] +print(f"REF: field correlation (aligned) = {rho_raw:.4f}") + +for r in (24, 42): + t0 = time.time() + WV = sym_nmf_gpu(V0, r, seed=1) + WT = sym_nmf_gpu(T0s, r, seed=1) + WV2 = sym_nmf_gpu(V0, r, seed=2) # stability check: same matrix, new seed + print(f"\n=== r={r} [{time.time()-t0:.1f}s] resid V {np.linalg.norm(V0-WV@WV.T)/np.linalg.norm(V0):.4f}" + f" T {np.linalg.norm(T0s-WT@WT.T)/np.linalg.norm(T0s):.4f}" + f" V(seed2) {np.linalg.norm(V0-WV2@WV2.T)/np.linalg.norm(V0):.4f}") + + nrm = lambda X: X / (np.linalg.norm(X, axis=0, keepdims=True) + 1e-12) + An, Bn, A2n = nrm(WV), nrm(WT), nrm(WV2) + + # (i) UNIQUENESS: same matrix, two seeds + S = An.T @ A2n; rr, cc = linear_sum_assignment(-S) + cs = np.sort(S[rr, cc])[::-1] + print(f" [uniqueness] symNMF(V,seed1) vs symNMF(V,seed2) matched cos: " + f"median {np.median(cs):.3f} >0.9 {(cs>0.9).sum()}/{r} >0.5 {(cs>0.5).sum()}/{r} min {cs.min():.3f}") + + # (ii) ORACLE cross-modal column match -- full distribution + Sx = An.T @ Bn[truth]; rr, cc = linear_sum_assignment(-Sx) + ocs = Sx[rr, cc]; colmap = cc[np.argsort(rr)] + srt = np.sort(ocs)[::-1] + print(f" [oracle cols] cos: >0.9 {(ocs>0.9).sum()}/{r} >0.5 {(ocs>0.5).sum()}/{r}" + f" median {np.median(ocs):.3f} min {ocs.min():.3f}") + print(f" sorted: {np.round(srt,2)}") + print(f" identity-map fraction (shared-seed artifact check): {(colmap==np.arange(r)).mean():.3f}") + + # (iii) row descriptors from oracle-matched columns + Ar = WV / np.linalg.norm(WV, axis=1, keepdims=True).clip(1e-9) + Br = WT[:, colmap] / np.linalg.norm(WT[:, colmap], axis=1, keepdims=True).clip(1e-9) + p_desc = hung(Ar, Br) + print(f" [oracle-col row descriptors] Hungarian acc {acc(p_desc):.4f} -> descent {acc(descend(p_desc)):.4f}") + + # (iv) THE PROPOSAL, ORACLE FORM: prune unmatched V components, rebuild, re-solve + for thr in (0.5, 0.7, 0.8): + keep = np.where(ocs >= thr)[0] + kv = rr[np.isin(np.arange(r), np.arange(r))] # rr is identity-ordered by lsa + keepV = rr[ocs >= thr] + if len(keepV) < 4: + print(f" [prune thr={thr}] only {len(keepV)} kept, skip"); continue + Vc = WV[:, keepV] @ WV[:, keepV].T + Vcs = standardise(Vc) + rho = np.corrcoef(Vcs[~np.eye(N, dtype=bool)], Ts[np.ix_(truth, truth)][~np.eye(N, dtype=bool)])[0, 1] + # solve on the cleaned vision field against the original text field + pg = grampa(Vcs, Ts) + Acg = torch.tensor(Vcs, device=dev, dtype=DT) + pgd = descend(pg, A=Acg) # descent on the CLEANED objective + pgd_true = descend(pgd) # then polish on the TRUE objective + print(f" [prune thr={thr}] kept {len(keepV)}/{r} rho(cleanV,T)={rho:.4f} (was {rho_raw:.4f})" + f" GRAMPA {acc(pg):.4f} -> clean-descent {acc(pgd):.4f} -> true-descent {acc(pgd_true):.4f}" + f" E={energy(pgd_true):.4f}") diff --git a/artifacts/spectral_frontier_probe/nmf_adv_prune_control.log b/artifacts/spectral_frontier_probe/nmf_adv_prune_control.log new file mode 100644 index 0000000..347318c --- /dev/null +++ b/artifacts/spectral_frontier_probe/nmf_adv_prune_control.log @@ -0,0 +1,21 @@ +--- symNMF seed 1, r=42 + A. symNMF reconstruction, NO pruning (blind) GW 0.0312 -> clean-desc 0.0391 -> true-desc 0.0391 E=0.6405 + B. ORACLE prune worst 1 (thr .5) GW 0.8594 -> clean-desc 0.8594 -> true-desc 0.8672 E=0.3396 + C0. RANDOM single-column drop (col 32, blind) GW 0.8359 -> clean-desc 0.8281 -> true-desc 0.8359 E=0.3396 + C1. RANDOM single-column drop (col 13, blind) GW 0.0469 -> clean-desc 0.0430 -> true-desc 0.0469 E=0.6486 + C2. RANDOM single-column drop (col 18, blind) GW 0.0195 -> clean-desc 0.0156 -> true-desc 0.0156 E=0.6661 + D. both sides symNMF-reconstructed (blind) GW 0.0234 -> clean-desc 0.0234 -> true-desc 0.0234 E=0.6020 +--- symNMF seed 2, r=42 + A. symNMF reconstruction, NO pruning (blind) GW 0.0156 -> clean-desc 0.0156 -> true-desc 0.0156 E=0.6419 + B. ORACLE prune worst 4 (thr .5) GW 0.8203 -> clean-desc 0.8125 -> true-desc 0.8203 E=0.3396 + C0. RANDOM single-column drop (col 32, blind) GW 0.0078 -> clean-desc 0.0078 -> true-desc 0.0078 E=0.6720 + C1. RANDOM single-column drop (col 13, blind) GW 0.0195 -> clean-desc 0.0195 -> true-desc 0.0195 E=0.7034 + C2. RANDOM single-column drop (col 18, blind) GW 0.0039 -> clean-desc 0.0039 -> true-desc 0.0039 E=0.6649 + D. both sides symNMF-reconstructed (blind) GW 0.0156 -> clean-desc 0.0156 -> true-desc 0.0156 E=0.7342 +--- symNMF seed 3, r=42 + A. symNMF reconstruction, NO pruning (blind) GW 0.1328 -> clean-desc 0.1523 -> true-desc 0.1641 E=0.5919 + B. ORACLE prune worst 2 (thr .5) GW 0.0000 -> clean-desc 0.0000 -> true-desc 0.0000 E=0.6669 + C0. RANDOM single-column drop (col 32, blind) GW 0.1211 -> clean-desc 0.1172 -> true-desc 0.1211 E=0.5922 + C1. RANDOM single-column drop (col 13, blind) GW 0.8125 -> clean-desc 0.8203 -> true-desc 0.8203 E=0.3396 + C2. RANDOM single-column drop (col 18, blind) GW 0.0078 -> clean-desc 0.0117 -> true-desc 0.0117 E=0.6752 + D. both sides symNMF-reconstructed (blind) GW 0.8398 -> clean-desc 0.8477 -> true-desc 0.8477 E=0.3396 diff --git a/artifacts/spectral_frontier_probe/nmf_adv_prune_control.py b/artifacts/spectral_frontier_probe/nmf_adv_prune_control.py new file mode 100644 index 0000000..a3b23ef --- /dev/null +++ b/artifacts/spectral_frontier_probe/nmf_adv_prune_control.py @@ -0,0 +1,61 @@ +"""CONTROL: is the 0.86 from PRUNING (needs oracle) or from symNMF RECONSTRUCTION alone (blind)?""" +import numpy as np, torch, warnings, ot +warnings.filterwarnings('ignore') +from scipy.optimize import linear_sum_assignment +dev = torch.device('cuda:0'); DT = torch.float64 +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 standardise(M): + M = np.asarray(M, float); mask = ~np.eye(len(M), dtype=bool); v = M[mask] + out = (M - v.mean()) / v.std(); np.fill_diagonal(out, 0.0); return out +def sym_nmf_gpu(M, r, iters=3000, seed=0): + Mg = torch.tensor(np.clip(M, 0, None), device=dev, dtype=DT) + g = torch.Generator(device='cpu').manual_seed(seed) + W = torch.abs(torch.randn(len(M), r, generator=g, dtype=DT)).to(dev) * float(np.sqrt(M.mean()/r)) + for _ in range(iters): + W = W * (0.5 + 0.5 * (Mg @ W) / (W @ (W.T @ W) + 1e-12)) + return W.cpu().numpy() +Vs = standardise(V0); Ts = standardise(T0s) +CONST = float((Vs*Vs).sum() + (Ts*Ts).sum()); DEN = N*(N-1) +Ag = torch.tensor(Vs, device=dev, dtype=DT); Tg = torch.tensor(Ts, device=dev, dtype=DT) +energy = lambda p: (CONST - 2.0*float((Ts[np.ix_(p, p)]*Vs).sum()))/DEN +def descend(p0, A=None, max_steps=4000): + A = Ag if A is None else A + p = torch.tensor(np.asarray(p0), device=dev, dtype=torch.long) + iu = torch.triu_indices(N, N, offset=1, device=dev) + for _ in range(max_steps): + B = Tg[p][:, p]; C = A @ B; dg = torch.diagonal(C) + G = C + C.T - dg[:, None] - dg[None, :] + 2*A*B + vals = G[iu[0], iu[1]]; k = int(vals.argmax()) + if float(vals[k]) <= 1e-12: break + u, v = int(iu[0][k]), int(iu[1][k]); p[u], p[v] = p[v].clone(), p[u].clone() + return p.cpu().numpy() +def gw(A, B): + q = np.ones(N)/N + G, _ = ot.gromov.gromov_wasserstein(A, B, q, q, 'square_loss', log=True, max_iter=200) + _, c = linear_sum_assignment(-G); return c +def run(tag, Vm, Tm=Ts): + p = gw(Vm, Tm); pd = descend(p, A=torch.tensor(Vm, device=dev, dtype=DT)); pd2 = descend(pd) + print(f" {tag:46s} GW {acc(p):.4f} -> clean-desc {acc(pd):.4f} -> true-desc {acc(pd2):.4f} E={energy(pd2):.4f}", flush=True) + return acc(pd2) + +r = 42 +for seed in (1, 2, 3): + print(f"--- symNMF seed {seed}, r={r}") + WV = sym_nmf_gpu(V0, r, seed=seed); WT = sym_nmf_gpu(T0s, r, seed=seed) + nrm = lambda X: X/(np.linalg.norm(X, axis=0, keepdims=True)+1e-12) + Sx = nrm(WV).T @ nrm(WT)[truth]; rr, cc = linear_sum_assignment(-Sx); ocs = Sx[rr, cc] + worst = rr[int(np.argmin(ocs))] + Vrec = standardise(WV @ WV.T) # BLIND: no pruning at all + run("A. symNMF reconstruction, NO pruning (blind)", Vrec) + keepV = rr[ocs >= 0.5] + run(f"B. ORACLE prune worst {r-len(keepV)} (thr .5)", standardise(WV[:, keepV] @ WV[:, keepV].T)) + for t in range(3): # BLIND: drop a random column + j = np.random.default_rng(100+t).integers(r) + kk = np.setdiff1d(np.arange(r), [j]) + run(f"C{t}. RANDOM single-column drop (col {j}, blind)", standardise(WV[:, kk] @ WV[:, kk].T)) + # also: text side reconstructed the same way (symmetric treatment) + run("D. both sides symNMF-reconstructed (blind)", Vrec, standardise(WT @ WT.T)) diff --git a/artifacts/spectral_frontier_probe/nmf_adv_prune_gw.py b/artifacts/spectral_frontier_probe/nmf_adv_prune_gw.py new file mode 100644 index 0000000..2ecdfed --- /dev/null +++ b/artifacts/spectral_frontier_probe/nmf_adv_prune_gw.py @@ -0,0 +1,50 @@ +"""Second solver on the ORACLE-pruned field, so the refutation is not GRAMPA-specific.""" +import numpy as np, torch, warnings, ot +warnings.filterwarnings('ignore') +from scipy.optimize import linear_sum_assignment +dev = torch.device('cuda:0'); DT = torch.float64 +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 standardise(M): + M = np.asarray(M, float); mask = ~np.eye(len(M), dtype=bool); v = M[mask] + out = (M - v.mean()) / v.std(); np.fill_diagonal(out, 0.0); return out +def sym_nmf_gpu(M, r, iters=3000, seed=0): + Mg = torch.tensor(np.clip(M, 0, None), device=dev, dtype=DT) + g = torch.Generator(device='cpu').manual_seed(seed) + W = torch.abs(torch.randn(len(M), r, generator=g, dtype=DT)).to(dev) * float(np.sqrt(M.mean() / r)) + for _ in range(iters): + W = W * (0.5 + 0.5 * (Mg @ W) / (W @ (W.T @ W) + 1e-12)) + return W.cpu().numpy() +Vs = standardise(V0); Ts = standardise(T0s) +CONST = float((Vs*Vs).sum() + (Ts*Ts).sum()); DEN = N*(N-1) +Ag = torch.tensor(Vs, device=dev, dtype=DT); Tg = torch.tensor(Ts, device=dev, dtype=DT) +energy = lambda p: (CONST - 2.0*float((Ts[np.ix_(p, p)]*Vs).sum()))/DEN +def descend(p0, A=None, max_steps=4000): + A = Ag if A is None else A + p = torch.tensor(np.asarray(p0), device=dev, dtype=torch.long) + iu = torch.triu_indices(N, N, offset=1, device=dev) + for _ in range(max_steps): + B = Tg[p][:, p]; C = A @ B; dg = torch.diagonal(C) + G = C + C.T - dg[:, None] - dg[None, :] + 2*A*B + vals = G[iu[0], iu[1]]; k = int(vals.argmax()) + if float(vals[k]) <= 1e-12: break + u, v = int(iu[0][k]), int(iu[1][k]); p[u], p[v] = p[v].clone(), p[u].clone() + return p.cpu().numpy() +def gw(A, B, seed=0): + q = np.ones(N)/N + G, _ = ot.gromov.gromov_wasserstein(A, B, q, q, 'square_loss', log=True, max_iter=200) + r_, c_ = linear_sum_assignment(-G); return c_ +r = 42 +WV = sym_nmf_gpu(V0, r, seed=1); WT = sym_nmf_gpu(T0s, r, seed=1) +nrm = lambda X: X/(np.linalg.norm(X, axis=0, keepdims=True)+1e-12) +Sx = nrm(WV).T @ nrm(WT)[truth]; rr, cc = linear_sum_assignment(-Sx); ocs = Sx[rr, cc] +print("baseline GW on raw:", end=" ") +p = gw(Vs, Ts); pd = descend(p); print(f"raw {acc(p):.4f} refined {acc(pd):.4f} E={energy(pd):.4f}") +for thr in (0.5, 0.65, 0.75): + keepV = rr[ocs >= thr] + Vc = standardise(WV[:, keepV] @ WV[:, keepV].T) + p = gw(Vc, Ts); pd = descend(p, A=torch.tensor(Vc, device=dev, dtype=DT)); pd2 = descend(pd) + print(f"prune thr={thr} kept {len(keepV)}/{r}: GW raw {acc(p):.4f} -> clean-desc {acc(pd):.4f} -> true-desc {acc(pd2):.4f} E={energy(pd2):.4f}") |