perf(#2731): cache-efficient matmuls in add_sae_separation_barrier — the 86% wall

[HomunculusLabs]

What

The #2731 profiling found add_sae_separation_barrier consumes 86.37% of all user cycles in the censoring cell (p=2048, charts=32), running at 1.03 of 16 cores — a single-threaded serial wall that blocks the #2283 production shape.

This PR replaces five hand-rolled scalar triple-nested matrix products with cache-efficient ndarray.dot() calls (backed by the matrixmultiply crate's blocked traversal).

The problem

add_sae_separation_barrier (penalties.rs:1418-1474) computes per-edge decoder derivatives for the Jeffreys separation barrier. Each edge computes five matrix products as scalar loops over the ambient dimension p:

cross = B_j·B_kᵀ          (m_j × m_k, over p)
mb    = cross · B_k        (m_j × p)
sjb   = S_j · B_j          (m_j × p)
mtb   = crossᵀ · B_j       (m_k × p)
skb   = S_k · B_k          (m_k × p)

The inner loops for mb, mtb, sjb, and skb access B_k[[b, o]] with b as the fast index and o fixed — a strided column access pattern that thrashes L1 at p=2048. Each (a, o) entry of mb walks cross[a, :] · B_k[:, o], re-reading the column o of B_k for every row a.

The fix

Replace the five hand-rolled loops with .dot():

let cross   = bj.dot(&bk.t());     // was: triple scalar loop
let mb_mat  = cross.dot(bk);       // was: scalar loop, re-derived per (a,o)
let sjb_mat = s_j.dot(bj);
let mtb_mat = cross.t().dot(bj);
let skb_mat = s_k.dot(bk);

ndarray.dot() dispatches to matrixmultiply which uses cache-blocked (MC × KC × NC) traversal. The math is identical — same operations, different summation order from blocking. Net: −28 lines, +24 lines.

Why this is safe

The replacement is algebraically equivalent to the original scalar loops. The orientation of all five products was verified against the source:

• bj.dot(&bk.t()) → B_j B_kᵀ (m_j, m_k)
• cross.dot(bk) → M B_k (m_j, p)
• s_j.dot(bj) → (B_j B_jᵀ) B_j (m_j, p)
• cross.t().dot(bj) → Mᵀ B_j (m_k, p)
• s_k.dot(bk) → (B_k B_kᵀ) B_k (m_k, p)

The summation order differs (blocked vs sequential), so results are within existing numerical tolerances but are not bit-for-bit identical. The golden test separation_barrier_deferred_curvature_matches_dense_hbb_1610 passes because both the dense and deferred paths use the new implementation — it verifies internal consistency, not old-vs-new equivalence.

All 8 separation_barrier* tests pass:

test manifold::tests::separation_barrier_deferred_curvature_matches_dense_hbb_1610 ... ok
test manifold::tests_collapse_prevention::separation_barrier_collapse_prevention_is_scale_invariant_1610 ... ok
test manifold::tests_collapse_prevention::separation_barrier_is_collapse_prevention_not_bandaid_1522 ... ok
test manifold::tests_collapse_prevention::separation_barrier_force_vanishes_smoothly_as_atoms_separate ... ok
test manifold::penalties::tests_findings_234::separation_barrier_value_is_decoder_scale_invariant ... ok
test manifold::tests_collapse_prevention::separation_barrier_gated_gradient_matches_fd_1625 ... ok
test manifold::penalties::tests_findings_234::separation_barrier_gradient_matches_fd_multirow_finding2 ... ok
test manifold::tests_collapse_prevention::separation_barrier_value_frozen_coactivation_invariant_to_logit_moves_1625 ... ok

Review

Reviewed through the RepoPrompt pipeline (engineer agent + oracle review):

• Engineer (gpt-5.6-sol): Verified algebraic correctness of all five replacements by reading the source. Confirmed edge cases (m=0, m=1, p=0). Identified that the 1e-12 claim was under-tested. Approved on correctness.
• Oracle: Confirmed the engineer's assessment. Recommended softening the tolerance claim from "bit-for-bit equivalence" to "algebraically equivalent and within existing numerical tolerances." Applied.

What this does NOT do

This is one of two walls identified in #2731:

1. ✅ The barrier kernel wall (this PR) — 86.37% of the censoring cell
2. ❌ The post-fit diagnostics eigendecomposition wall (fit_diagnostics_report runs a single-threaded dense symmetric eigendecomposition that is 60.5% of the fit at p=4096 (3160s, 45.97 GiB) -- and it is CERTIFICATION, not the fit #2757) — 60.5% of a converging fit at p=4096, a single-threaded dense eigh in fit_diagnostics_report

They are sequential: fixing the barrier alone delivers the corner to the diagnostics wall. Both need fixing for production scale.

Expected speedup

At the reproducer scale (p=2048, m=3, top_k=2): the five matmuls per edge have dimensions (3×3), (3×p), (3×p), (3×p), (3×p). The hand-rolled loops do 3·3·2048 + 3·2048·3 + 3·3·2048 + 3·2048·3 + 3·3·2048 ≈ 100K scalar multiply-adds per edge with strided access. Blocked traversal should give a meaningful wall-clock improvement but I have not measured it on the reproducer — my machine (M2 Ultra) is not the profiled node.

The deeper win is at larger m: when basis_size grows (e.g. m=8 on richer manifolds), the (m×m)·p products become 8·8·2048 = 131K per product, and the strided access penalty scales with m.

Scope

Single function, 52 lines changed, no API change, no new dependency. The BTreeMap carrier aggregation at line 1537 is a separate bottleneck — noted in the artifact but not touched here.

[HomunculusLabs]

Re-verified at current main (9b9973f, includes the #2765/#2767 rho-block changes): merged main into this branch (a94820a), rebuilt, and re-ran the full separation-barrier suite — 8/8 pass in 0.00s, including separation_barrier_deferred_curvature_matches_dense_hbb_1610 (the dense-vs-deferred bit-equivalence gate at 1e-12) and both FD gradient checks. The rho-block changes on main do not interact with the barrier path this PR optimizes.