The missing term was unary: omit-size solved at its ceiling in 1.4s

Adversarial review of the artifacts found three of my numbers to be
artifacts of my own code. All three reproduced here before acceptance:

- anchor_bound presented probe rows in the same index order on both
  sides, so exact twins had their tie broken onto the diagonal. 0.997 ->
  0.920 on omit-size. Fixed by scrambling the T-side presentation.
- The truth is not a strict local minimum: 51 transpositions have
  exactly zero energy delta. fast_pair_descent only looked stationary
  because its break test treats zero as no-improvement.
- scipy's FAQ takes no n_init, so it was swallowed into unknown_options
  and 'FAQ x30 restarts' computed bit-identically to plain FAQ. Replaced
  with a real restart loop over P0='randomized'.

The blind ceiling for omit-size is 0.836, not 1.0: the text field has 51
exact transposition automorphisms, so T[s,s] is bitwise identical to T
and no objective f(V, P T P^T) can separate an orbit at any order. Every
synthetic accuracy was being divided by the wrong denominator.

The fifteen failed solvers share one property -- all purely quadratic or
purely spectral, none with a node-level term. Eight moments of each
node's own field row, blended with the quadratic term through
Frank-Wolfe, reach 0.837 with an energy gap of exactly zero. The term
must stay in the loop: as a seed for pure-quadratic descent it scores
0.21, pinned through the iterations it scores 0.84 -- which is also why
amplification plateaued, being itself pure-quadratic.

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

4 files changed, 333 insertions, 9 deletions

diff --git a/worldalign/anchor_bound.py b/worldalign/anchor_bound.py
index 49f4ed3..2920934 100644
--- a/worldalign/anchor_bound.py
+++ b/worldalign/anchor_bound.py
@@ -79,13 +79,19 @@ def analyse(path: str, label: str, repeats: int, blind: float | None) -> dict:
         order = generator.permutation(size)
         half = size // 2
         anchors, probe = order[:half], order[half:]
+        # The two sides must not be presented in the same index order. Exact
+        # twins have identical anchor-restricted rows, and with matched
+        # ordering `linear_sum_assignment` breaks the tie onto the diagonal --
+        # crediting the bound with information it does not have. On the
+        # caption-omitted field that inflated it from 0.920 to 0.997.
+        order = generator.permutation(len(probe))
         left = rows_against(visual, probe, anchors)
-        right = rows_against(text, probe, anchors)
+        right = rows_against(text, probe[order], anchors)
         similarity = left @ right.T / left.shape[1]
-        truth = np.arange(len(probe))
-        nearest.append(float((similarity.argmax(1) == truth).mean()))
+        truth = order
+        nearest.append(float((order[similarity.argmax(1)] == np.arange(len(probe))).mean()))
         _, columns = linear_sum_assignment(-similarity)
-        assigned.append(float((columns == truth).mean()))
+        assigned.append(float((order[columns] == np.arange(len(probe))).mean()))
 
     row = {
         "label": label,
diff --git a/worldalign/automorphism_probe.py b/worldalign/automorphism_probe.py
new file mode 100644
index 0000000..46949f4
--- /dev/null
+++ b/worldalign/automorphism_probe.py
@@ -0,0 +1,96 @@
+"""Exact automorphisms of a relation field, and the blind-recovery ceiling.
+
+If sigma is an exact automorphism of the text field T (P_sigma T P_sigma^T = T
+bit-for-bit) then the blind matching problem cannot distinguish the truth from
+truth o sigma: the input (V, T) is literally unchanged.  Every functional of
+(V, P T P^T) -- the pairwise energy, the third-order tr(M^3) term, any moment
+in the ladder, the TLB linear term -- is invariant, so no optimiser and no
+relaxation can break the tie.  The posterior over the truth is uniform on the
+orbit, and the expected accuracy of *any* blind estimator is bounded by
+
+    ceiling = (1/n) * sum_over_orbits  |orbit| * (1/|orbit|)
+            = 1 - (moved - #orbits) / n
+
+Run:  python -m worldalign.automorphism_probe artifacts/synth_v1/omit_size.pt
+"""
+
+from __future__ import annotations
+
+import math
+import sys
+
+import numpy as np
+import torch
+
+
+def standardise(matrix: np.ndarray) -> np.ndarray:
+    matrix = np.asarray(matrix, dtype=np.float64)
+    mask = ~np.eye(len(matrix), dtype=bool)
+    out = (matrix - matrix[mask].mean()) / matrix[mask].std()
+    np.fill_diagonal(out, 0.0)
+    return out
+
+
+def automorphic_transpositions(field: np.ndarray, tol: float = 1e-9) -> np.ndarray:
+    """adj[i, j] iff swapping i and j leaves the field exactly invariant.
+
+    Rows i and j must agree on every coordinate outside {i, j}; the two
+    excluded coordinates are exactly the ones the swap moves.
+    """
+    size = len(field)
+    mismatch = np.full((size, size), np.inf)
+    for i in range(size):
+        gap = np.abs(field - field[i])
+        gap[:, i] = 0.0
+        gap[np.arange(size), np.arange(size)] = 0.0
+        mismatch[i] = gap.max(axis=1)
+    np.fill_diagonal(mismatch, np.inf)
+    return mismatch < tol
+
+
+def orbits(adjacency: np.ndarray) -> list[list[int]]:
+    size = len(adjacency)
+    seen = np.zeros(size, dtype=bool)
+    found: list[list[int]] = []
+    for start in range(size):
+        if seen[start]:
+            continue
+        stack = [start]
+        seen[start] = True
+        component = [start]
+        while stack:
+            node = stack.pop()
+            for nxt in np.nonzero(adjacency[node] & ~seen)[0]:
+                seen[nxt] = True
+                stack.append(nxt)
+                component.append(int(nxt))
+        found.append(sorted(component))
+    return [c for c in found if len(c) > 1]
+
+
+def report(field: np.ndarray, name: str) -> float:
+    size = len(field)
+    groups = orbits(automorphic_transpositions(field))
+    moved = sum(len(c) for c in groups)
+    log_order = sum(math.lgamma(len(c) + 1) for c in groups) / math.log(10)
+    ceiling = 1.0 - (moved - len(groups)) / size
+    print(
+        f"{name}: {len(groups)} nontrivial orbits covering {moved}/{size} items, "
+        f"sizes {sorted((len(c) for c in groups), reverse=True)[:12]}, "
+        f"|Aut| >= 10^{log_order:.1f}"
+    )
+    print(f"{name}: blind accuracy ceiling for any energy-only method = {ceiling:.4f}")
+    return ceiling
+
+
+def main() -> None:
+    path = sys.argv[1] if len(sys.argv) > 1 else "artifacts/synth_v1/omit_size.pt"
+    state = torch.load(path, map_location="cpu", weights_only=False)
+    visual = standardise(state["visual_field"].double().numpy())
+    text = standardise(state["text_field"].double().numpy())
+    report(text, "text field")
+    report(visual, "visual field")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/fused_match.py b/worldalign/fused_match.py
new file mode 100644
index 0000000..aad3360
--- /dev/null
+++ b/worldalign/fused_match.py
@@ -0,0 +1,214 @@
+"""The unary term the search was missing.
+
+Fifteen solvers returned under 5% on the caption-omitted field: Umeyama,
+GRAMPA, five Gromov-Wasserstein variants, FAQ, PATH convex-concave, a moment
+ladder, consensus voting. They have one thing in common, and it is the thing
+that was wrong with all of them at once. **Every one is purely quadratic or
+purely spectral.** None carries a node-level term. The instance was built by
+deleting one factor from the caption side, which flattens exactly the coarse
+second-order statistics those methods read, while leaving scenes distinguishable
+in the shape of their individual similarity profiles.
+
+So give the objective a unary term. Each node gets a permutation-invariant
+descriptor -- moments of its own row of the relation field, which is a
+row-order statistic and not a spectral one -- and the descriptor distance
+becomes a linear cost blended with the quadratic one. The blend must stay in
+the loop: using the descriptor only to seed a pure-quadratic descent scores
+0.21, because the descent immediately discards it and wanders back into the
+decoy region. Kept pinned through Frank-Wolfe, the same descriptor scores 0.84.
+
+That 0.84 is the exact answer. The caption-omitted text field has 51 exact
+transposition automorphisms, so `T[s,s]` is bitwise identical to `T` for each
+and no objective built from `(V, P T P^T)` can distinguish the members of an
+orbit at any order. The blind ceiling is 0.836, and this reaches it with an
+energy gap of exactly zero.
+
+Alpha is chosen by energy, never by accuracy, so the procedure stays blind.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+from scipy.optimize import linear_sum_assignment
+
+from .common import write_json
+from .synth_fast_gate import ClosedFormEnergy, fast_pair_descent
+from .synth_triangle_gate import standardized
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", nargs="+", required=True)
+    parser.add_argument("--labels", nargs="+", default=None)
+    parser.add_argument("--trials", type=int, default=3)
+    parser.add_argument("--iterations", type=int, default=40)
+    parser.add_argument("--device", default="cuda:3")
+    parser.add_argument("--output", default="artifacts/synth_v1/fused_match.json")
+    return parser.parse_args()
+
+
+def standardise(matrix: np.ndarray) -> np.ndarray:
+    mask = ~np.eye(len(matrix), dtype=bool)
+    values = matrix[mask]
+    out = (matrix - values.mean()) / values.std()
+    np.fill_diagonal(out, 0.0)
+    return out
+
+
+def row_descriptors(field: np.ndarray) -> np.ndarray:
+    """Moments of each node's own similarity profile.
+
+    Invariant to how the other nodes are ordered, so it survives the unknown
+    permutation, and it reads the distribution of a node's relations rather
+    than the field's leading eigenvectors -- which is why it sees what the
+    spectral family cannot.
+    """
+    size = len(field)
+    rows = []
+    for index in range(size):
+        row = np.delete(field[index], index)
+        rows.append([
+            row.mean(), row.std(), np.mean(row ** 3), np.mean(row ** 4),
+            *np.percentile(row, [10, 50, 90]), row.max(),
+        ])
+    return np.array(rows)
+
+
+def unary_cost(visual: np.ndarray, text: np.ndarray) -> np.ndarray:
+    """Squared descriptor distance, z-scored jointly over both node sets."""
+    left, right = row_descriptors(visual), row_descriptors(text)
+    both = np.vstack([left, right])
+    centre, spread = both.mean(0), both.std(0).clip(1e-9)
+    left, right = (left - centre) / spread, (right - centre) / spread
+    return ((left[:, None, :] - right[None, :, :]) ** 2).sum(-1)
+
+
+def blind_automorphism_ceiling(field: np.ndarray) -> tuple[float, int]:
+    """The best any blind method can do, given exact automorphisms of the field.
+
+    Two scenes whose rows agree as multisets may be exchangeable: if swapping
+    them leaves the field bitwise unchanged, every objective of the form
+    f(V, P T P^T) is blind to which is which, at every order. Members of an
+    orbit are then guessable only up to one correct pick per orbit.
+    """
+    size = len(field)
+    buckets: dict[bytes, list[int]] = {}
+    for index in range(size):
+        row = np.sort(np.delete(field[index], index))
+        buckets.setdefault(row.tobytes(), []).append(index)
+    orbits = [group for group in buckets.values() if len(group) > 1]
+    automorphisms = 0
+    for group in orbits:
+        for position, first in enumerate(group):
+            for second in group[position + 1:]:
+                swap = np.arange(size)
+                swap[[first, second]] = swap[[second, first]]
+                if np.array_equal(field[np.ix_(swap, swap)], field):
+                    automorphisms += 1
+    ambiguous = sum(len(group) for group in orbits)
+    ceiling = (size - ambiguous + len(orbits)) / size
+    return float(ceiling), automorphisms
+
+
+def fused_frank_wolfe(
+    visual: np.ndarray, text: np.ndarray, cost: np.ndarray,
+    alpha: float, iterations: int,
+) -> np.ndarray:
+    """Maximise alpha*<V, P T P^T> - (1-alpha)*<cost, P> over the Birkhoff polytope."""
+    size = len(visual)
+    coupling = np.full((size, size), 1.0 / size)
+    scale = abs(visual @ coupling @ text).mean() / max(abs(cost).mean(), 1e-12)
+    scaled = cost * scale
+    for _ in range(iterations):
+        gradient = alpha * 2.0 * (visual @ coupling @ text) - (1.0 - alpha) * scaled
+        _, columns = linear_sum_assignment(-gradient)
+        vertex = np.zeros_like(coupling)
+        vertex[np.arange(size), columns] = 1.0
+        direction = vertex - coupling
+        quadratic = alpha * float((direction * (visual @ direction @ text)).sum())
+        linear = float((direction * gradient).sum())
+        if quadratic >= 0:
+            step = 1.0
+        else:
+            step = float(np.clip(-linear / (2.0 * quadratic), 0.0, 1.0))
+        if step <= 1e-12:
+            break
+        coupling = coupling + step * direction
+    _, columns = linear_sum_assignment(-coupling)
+    return columns
+
+
+def main() -> None:
+    args = parse_args()
+    labels = args.labels or [p.split("/")[-1] for p in args.fields]
+    device = torch.device(args.device)
+    rows = []
+
+    for path, label in zip(args.fields, labels):
+        state = torch.load(path, map_location="cpu", weights_only=False)
+        visual = standardise(state["visual_field"].double().numpy())
+        text = standardise(state["text_field"].double().numpy())
+        size = len(visual)
+        ceiling, automorphisms = blind_automorphism_ceiling(text)
+        visual_gpu = standardized(torch.from_numpy(visual).to(device)).double()
+        print(f"\n=== {label}  N={size}  automorphisms={automorphisms}  "
+              f"blind ceiling={ceiling:.4f}", flush=True)
+
+        accuracies, gaps = [], []
+        for trial in range(args.trials):
+            hidden = np.random.default_rng(trial).permutation(size)
+            shuffled = text[np.ix_(hidden, hidden)]
+            text_gpu = standardized(torch.from_numpy(shuffled).to(device)).double()
+            energy = ClosedFormEnergy(text_gpu, visual_gpu, 1.0, 0.0, 256)
+            truth_energy = float(energy.energy(
+                torch.from_numpy(np.argsort(hidden).copy()).to(device)[None]
+            )[0])
+            cost = unary_cost(visual, shuffled)
+
+            best = (np.inf, 0.0, 0.0)
+            for alpha in np.arange(0.70, 0.99, 0.04):
+                columns = fused_frank_wolfe(
+                    visual, shuffled, cost, float(alpha), args.iterations
+                )
+                final = fast_pair_descent(
+                    text_gpu, visual_gpu,
+                    torch.from_numpy(np.ascontiguousarray(columns)).to(device), 2000,
+                )
+                value = float(energy.energy(final[None])[0])
+                if value < best[0]:                       # selection by energy only
+                    best = (value, float(
+                        (hidden[final.cpu().numpy()] == np.arange(size)).mean()
+                    ), float(alpha))
+            accuracies.append(best[1])
+            gaps.append(best[0] - truth_energy)
+            print(f"  trial {trial}: accuracy={best[1]:.4f}  alpha={best[2]:.2f}  "
+                  f"dE={best[0] - truth_energy:+.2e}", flush=True)
+
+        row = {
+            "field": label, "size": size, "chance": 1.0 / size,
+            "automorphisms": automorphisms, "blind_ceiling": ceiling,
+            "accuracy": float(np.mean(accuracies)),
+            "fraction_of_ceiling": float(np.mean(accuracies)) / ceiling,
+            "energy_gap": float(np.mean(gaps)),
+        }
+        rows.append(row)
+        print(f"  MEAN accuracy={row['accuracy']:.4f} = "
+              f"{row['fraction_of_ceiling']:.3f} of the blind ceiling", flush=True)
+
+    write_json(args.output, {
+        "protocol": (
+            "Blind. Hidden permutation generated per trial and read only for "
+            "scoring; alpha selected by energy. The ceiling is computed from "
+            "exact automorphisms of the text field and bounds any blind method."
+        ),
+        "rows": rows,
+    })
+    print(json.dumps({"done": True}))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/match_battery.py b/worldalign/match_battery.py
index 9035e62..63a4bf1 100644
--- a/worldalign/match_battery.py
+++ b/worldalign/match_battery.py
@@ -142,11 +142,19 @@ def solve_faq_multi(visual, shuffled, **_):
     """FAQ from many random doubly stochastic starts, best by its own objective."""
     from scipy.optimize import quadratic_assignment
 
-    result = quadratic_assignment(
-        visual, shuffled, method="faq",
-        options={"maximize": True, "n_init": 30, "rng": 0},
-    )
-    return result.col_ind
+    # scipy's FAQ takes no `n_init`; passing one is silently swallowed into
+    # unknown_options, which is why an earlier version of this function
+    # returned bit-identical results to plain FAQ. Restart explicitly.
+    best, best_value = None, -np.inf
+    for seed in range(15):
+        result = quadratic_assignment(
+            visual, shuffled, method="faq",
+            options={"maximize": True, "rng": seed, "P0": "randomized"},
+        )
+        value = float((visual * shuffled[np.ix_(result.col_ind, result.col_ind)]).sum())
+        if value > best_value:
+            best, best_value = result.col_ind, value
+    return best
 
 
 def solve_2opt(visual, shuffled, **_):