Cluster matching and Gaussian null kernels: reuse assessment and reference gate

Status: ABI 2.20.0 (unreleased), n4m_cluster_jaccard_matrix, n4m_linear_sum_assignment, n4m_random_permutations and n4m_random_multivariate_normal in n4m/metrics.h. Requested by nirs4all-preflight lot PF-05 (cluster appearance frequency and the unimodal Gaussian null model). The kernels are generic; the matching threshold, the covariance regularization and the null protocol belong to the caller.

Reuse assessment (2026-10-13)

Need	Candidates	Decision
One-to-one assignment	SciPy 1.18.1 scipy/optimize/rectangular_lsap/rectangular_lsap.cpp (BSD-3-Clause, C++, about 250 lines, Crouse’s shortest augmenting path) is the reference itself. gatagat/lap LAPJV (BSD-2-Clause) solves the same problem with another tie order. dlib max_cost_assignment (Boost) is integer-only. The Rust pathfinding Kuhn–Munkres (MIT/Apache-2.0) cannot be linked into C++. Methods had no assignment solver.	Port SciPy into cpp/src/core/optimization/_vendored/, keeping its control flow, tie rules and BSD-3 header (LICENSE.scipy.txt, THIRD_PARTY_LICENSES.md, NOTICE.md).
Jaccard between clusters	scikit-learn contingency_matrix (Python); no C/C++ library worth a dependency.	New native code (about 30 lines); scikit-learn is the test reference.
Subsampling plans	Methods PCG64 already reproduced NumPy Generator.permutation (n4m_pcg64_engine_permutation, used by the factor-association test).	Expose it as n4m_random_permutations: successive permutations from one seeded stream.
Gaussian draws	Methods PCG64 with the NumPy Ziggurat (n4m_pcg64_engine_standard_normal_fill, bit-identical to default_rng(seed).standard_normal, frozen by test_rng_pcg64.cpp). The Cholesky routines in gpr_pls.cpp, gp.cpp and emsc.c are file-local to those methods.	Reuse the PCG64 Ziggurat stream; a 15-line lower Cholesky lives with the kernel. No new normal sampler.

Frozen semantics

• Domains of validity are stated in each header. Real-valued inputs (assignment costs, Gaussian mean and covariance) follow the rule-13 domain; Gaussian draws factor a scaled covariance; Jaccard and ARI take integer labels (n ≤ 2^53, n ≤ 2^23). n4m_cluster_robustness_tests checks scaling, translation and constant-column invariance for these kernels, HDBSCAN and ARI.

• n4m_cluster_jaccard_matrix: negative labels are unassigned and belong to no cluster; labels must be below the declared cluster counts. J_ij = |A_i ∩ B_j| / (|A_i| + |B_j| - |A_i ∩ B_j|), NaN when both clusters are empty. Callers restrict both labelings to the elements they share before the call.

• n4m_linear_sum_assignment: the algorithm and tie rules of scipy.optimize.linear_sum_assignment, evaluated in exact arithmetic. Inside the rules 13/15 domain (nonzero |cost| in [2^-256, 2^256], max/min nonzero |cost| ≤ 2^32, +inf forbidden pairs excluded) every cost is an integer multiple of 2^E, E the last mantissa bit of the smallest nonzero |cost|, below 2^87 units; duals, reduced and path costs are kept as signed 192-bit integers in those units, with a guard far above any reachable value. The result is an exact optimum. On random 2–4 × 2–6 matrices with signed costs whose exponents span ±60 to ±600, SciPy’s double version is suboptimal in 2–9 % of cases; n4m_cluster_robustness_tests and the JS/WASM test check the exact optimum by enumeration with independent exact sums (2240-bit integers, BigInt) on in-domain matrices with full 53-bit mantissas and exponents mixed over 32 binades, for minimisation, maximisation and transposed shapes. Costs are never centred; inputs outside the domain are refused before computing. Work is bounded by rows·cols·min(rows, cols) ≤ 1e8 fixed-size operations and memory by 16·rows·cols + 128·(rows + cols) bytes ≤ 256 MiB, both checked before any allocation (a 1 × 32,000,000 matrix is refused). A tall matrix is solved on its transpose and the pairs are sorted by row; maximize negates the costs; +inf forbids a pair; NaN or -inf (after negation) is invalid; a missing feasible assignment is a numerical failure. Ties: in each augmenting-path search the columns are scanned from the last to the first, a strictly smaller reduced cost wins, and on equal reduced cost an unassigned column wins. A constant matrix therefore gives the identity assignment.

• n4m_random_permutations: one PCG64(seed) stream; row b of the output is the b-th rng.permutation(n) (Fisher–Yates from the last position over NumPy random_interval).

• n4m_random_multivariate_normal: z holds n × d standard normals from PCG64(seed) (NumPy SeedSequence expansion) in row-major order, L is the lower Cholesky factor of the lower triangle of cov, and x_ij = sum_k z_ik L_jk + mean_j. This is NumPy Generator.multivariate_normal(mean, cov, size=n, method="cholesky"); NumPy’s default method="svd" uses another factor and gives other draws. A covariance that is not numerically positive definite fails instead of being repaired. Storage and cost are checked before any allocation, independently of n: d·d and n·d at most 5e7 values, d^3/6 and n·d·d at most 2.5e9 operations.

Reference gate

parity/generators/cluster_matching_null_v1.py (test-only; pins NumPy 2.4.6, SciPy 1.18.1, scikit-learn 1.9.1 and hdbscan 0.8.44) writes parity/fixtures/cluster_matching_null_v1.json. Its SHA-256 at generation is ff70ef7e34cdb15d1d018e1adf88cee4a04941e032d7123fd8a1d7d018de32a1. The cases are:

• twelve assignment matrices: square, wide, tall (transposed path), small-integer ties, constant (identity), tall constant maximized, forbidden +inf pairs, a maximized Jaccard matrix, a single cell, and the NaN, maximized +inf and infeasible errors;

• four Jaccard pairs with noise, merged and new clusters, split clusters and an empty cluster slot, against scikit-learn contingency_matrix;

• four permutation streams, including 50 permutations of 40 units, a seed above 2^63, a single unit and zero permutations;

• four NumPy normal draws, including a full 3×3 covariance, seeds above 2^63 and 2^64 - 1, a ridge-regularized diagonal covariance and a univariate case;

• three Gaussian datasets (n = 80, d = 5, min_cluster_size = min_samples = 5) drawn by NumPy and clustered by the stable-tie scikit-learn reference of hdbscan_v1.py. One yields three clusters with 59 noise points and two yield none; the stock scikit-learn ARI is recorded for information (1, 1 and 0: stock scikit-learn finds clusters in the third dataset because of its unspecified MST tie order).

Three large-n adjusted Rand indices (up to 3 million elements, values near 0) are checked within 1e-14 relative against exact-integer scikit-learn. Five in-domain multi-magnitude HDBSCAN cases (see docs/dev/hdbscan.md) use scikit-learn HDBSCAN(metric="precomputed") on exact math.dist distances.

n4m_cluster_matching_null_tests checks the assignments, Jaccard values and permutations exactly, the draws within 1e-12 relative, the HDBSCAN partitions of the Methods draws up to label permutation, column-major cost storage, the lower-triangle-only covariance read, the argument errors, and the review’s overlapping 65536 × 65536 covariance view with zero draws, which is refused before any allocation (natively and in WASM). The JS/WASM test replays the same fixture through linearSumAssignment, clusterJaccardMatrix, randomPermutations, randomMultivariateNormal and clusterHdbscan.