fix(rust): deterministic clustered vcov — first-appearance cluster aggregation

The cluster-score aggregation in compute_robust_vcov (also reached via
solve_ols(return_vcov=True)) built its (G,k) cluster-scores matrix in
HashMap iteration order, which is SipHash-randomized per call:
mathematically identical, but the GEMM accumulation order changed on
every invocation, wobbling the clustered vcov at ~1e-14 (3 distinct
values across 8 identical calls; the Python backend is bit-stable).
Rows now accumulate in first-appearance order — ascending for the
factorized 0..G-1 ids the Python dispatcher passes, matching NumPy's
groupby order. Verified 1 distinct value across 50 identical calls;
NumPy parity unchanged (~1e-14 cross-backend GEMM tolerance).
Bit-identity regression tests cover contiguous and non-contiguous
unsorted cluster ids. Resolves the TODO row currently in flight on
PR #647 (row deletion follows once both merge).

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

Author: igerber

CHANGELOG.md

@@ -35,6 +35,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
     `did_had_pretest_workflow`, ...) are unchanged in this release and removed separately.
 
 ### Fixed
+- **Rust clustered vcov is now run-to-run deterministic.** The cluster-score aggregation
+  in `compute_robust_vcov` (also reached via `solve_ols(return_vcov=True)`) built its
+  (G, k) cluster-scores matrix in `HashMap` iteration order, which is SipHash-randomized
+  per call — mathematically identical, but the GEMM accumulation order changed on every
+  invocation, wobbling the vcov at ~1e-14 (3 distinct values observed across 8 identical
+  calls; the Python backend was bit-stable). Rows now accumulate in first-appearance
+  order — ascending for the factorized 0..G-1 ids the Python dispatcher passes, matching
+  NumPy's groupby order. Verified: 1 distinct value across 50 identical calls (was 3/8);
+  NumPy parity unchanged (~1e-14, the normal cross-backend GEMM tolerance); bit-identity
+  regression tests cover both contiguous and non-contiguous unsorted cluster ids.
 - **`HonestDiD` Δ^SD optimal-FLCI center parity with R (SE-audit B2b).** The optimal
   Fixed-Length CI optimizer was a flat Nelder-Mead over slope weights that landed on a
   different affine estimator than R `HonestDiD::findOptimalFLCI` at intermediate smoothness

rust/src/linalg.rs

@@ -349,18 +349,22 @@ fn compute_robust_vcov_internal(
             let residuals_col = residuals.insert_axis(Axis(1)); // (n, 1)
             let scores = x * &residuals_col; // (n, k) - broadcasts residuals across columns
 
-            // Aggregate scores by cluster using HashMap
-            let mut cluster_sums: HashMap<i64, Array1<f64>> = HashMap::new();
+            // Aggregate scores by cluster DETERMINISTICALLY. HashMap
+            // iteration order is SipHash-randomized per map instance, which
+            // reordered the cluster rows on every call — mathematically
+            // identical, but the GEMM accumulation order changed, making
+            // the clustered vcov run-to-run nondeterministic at ~1e-14
+            // (distinct values across identical calls; the Python backend
+            // is bit-stable). Rows accumulate in first-appearance order
+            // instead: for the factorized 0..G-1 ids the Python dispatcher
+            // passes this is ascending id order, matching NumPy's groupby.
+            let mut cluster_index: HashMap<i64, usize> = HashMap::new();
             for i in 0..n_obs {
-                let cluster = clusters[i];
-                let row = scores.row(i).to_owned();
-                cluster_sums
-                    .entry(cluster)
-                    .and_modify(|sum| *sum = &*sum + &row)
-                    .or_insert(row);
+                let next = cluster_index.len();
+                cluster_index.entry(clusters[i]).or_insert(next);
             }
 
-            let n_clusters = cluster_sums.len();
+            let n_clusters = cluster_index.len();
 
             if n_clusters < 2 {
                 return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(format!(
@@ -369,10 +373,13 @@ fn compute_robust_vcov_internal(
                 )));
             }
 
-            // Build cluster scores matrix (G, k)
+            // Build cluster scores matrix (G, k) in first-appearance order
             let mut cluster_scores = Array2::<f64>::zeros((n_clusters, k));
-            for (idx, (_cluster_id, sum)) in cluster_sums.iter().enumerate() {
-                cluster_scores.row_mut(idx).assign(sum);
+            for i in 0..n_obs {
+                let idx = cluster_index[&clusters[i]];
+                cluster_scores
+                    .row_mut(idx)
+                    .zip_mut_with(&scores.row(i), |a, b| *a += *b);
             }
 
             // Compute meat: Σ_g (X_g' e_g)(X_g' e_g)'

tests/test_rust_backend.py

@@ -3622,3 +3622,36 @@ def _unstable(*a, **k):
             v = la.compute_robust_vcov(X, resid, vcov_type="hc2")
         v_py = la._compute_robust_vcov_numpy(X, resid, None, vcov_type="hc2")
         np.testing.assert_allclose(v, v_py, rtol=1e-12, atol=1e-15)
+class TestClusterVcovDeterminism:
+    """Clustered vcov is bit-identical across repeated identical calls.
+
+    The cluster-score aggregation previously built rows in HashMap iteration
+    order (SipHash-randomized per call): mathematically identical, but the
+    GEMM accumulation order changed run-to-run, wobbling the vcov at ~1e-14
+    (3 distinct values observed in 8 identical calls) while the Python
+    backend was bit-stable. Rows now accumulate in first-appearance order
+    (ascending for the factorized ids the dispatcher passes)."""
+
+    def test_solve_ols_cluster_vcov_bit_identical_across_calls(self):
+        from diff_diff.linalg import solve_ols
+
+        rng = np.random.default_rng(0)
+        X = rng.normal(size=(400, 3))
+        y = rng.normal(size=400)
+        cl = np.repeat(np.arange(10), 40)
+        baseline = solve_ols(X, y, cluster_ids=cl, return_vcov=True)[2]
+        for _ in range(20):
+            v = solve_ols(X, y, cluster_ids=cl, return_vcov=True)[2]
+            np.testing.assert_array_equal(v, baseline)
+
+    def test_compute_robust_vcov_cluster_bit_identical_across_calls(self):
+        from diff_diff.linalg import compute_robust_vcov
+
+        rng = np.random.default_rng(1)
+        X = rng.normal(size=(300, 4))
+        resid = rng.normal(size=300)
+        # Non-contiguous, unsorted ids exercise the first-appearance remap.
+        cl = np.repeat(np.array([7, 3, 11, 5, 42, 3, 7, 11, 5, 42]), 30)
+        baseline = compute_robust_vcov(X, resid, cluster_ids=cl)
+        for _ in range(20):
+            np.testing.assert_array_equal(compute_robust_vcov(X, resid, cluster_ids=cl), baseline)