diff --git a/CHANGELOG.md b/CHANGELOG.md index 2f7c8dab..1b3e59dc 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/rust/src/linalg.rs b/rust/src/linalg.rs index c9940277..c808a967 100644 --- a/rust/src/linalg.rs +++ b/rust/src/linalg.rs @@ -281,18 +281,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> = 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 = 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::(format!( @@ -301,10 +305,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::::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)' diff --git a/tests/test_rust_backend.py b/tests/test_rust_backend.py index 210ebb84..279edb40 100644 --- a/tests/test_rust_backend.py +++ b/tests/test_rust_backend.py @@ -3507,3 +3507,39 @@ def test_shape_validation_errors(self): _batched_chol_symbol(np.zeros((2, 3, 4)), np.zeros(2)) with pytest.raises(ValueError, match="ridge length"): _batched_chol_symbol(np.zeros((2, 3, 3)), np.zeros(5)) + + +@pytest.mark.skipif(not HAS_RUST_BACKEND, reason="Rust backend not available") +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)