diff --git a/CHANGELOG.md b/CHANGELOG.md index 09b45f9c..245bd4ef 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 29db7c80..246b2ab2 100644 --- a/rust/src/linalg.rs +++ b/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> = 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!( @@ -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::::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 1275ae24..485f4f6f 100644 --- a/tests/test_rust_backend.py +++ b/tests/test_rust_backend.py @@ -3622,3 +3622,55 @@ 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) + + def test_raw_kernel_noncontiguous_ids_bit_identical(self): + """Review P3: the public wrapper factorizes cluster ids before Rust, + so the test above never hands the RAW kernel non-contiguous labels. + Exercise _rust_backend.compute_robust_vcov directly with unsorted, + non-contiguous int64 ids (the first-appearance remap path).""" + from diff_diff._rust_backend import compute_robust_vcov as raw_vcov + + rng = np.random.default_rng(5) + X = np.ascontiguousarray(rng.normal(size=(240, 3))) + resid = np.ascontiguousarray(rng.normal(size=240)) + cl = np.ascontiguousarray( + np.repeat(np.array([907, 3, -11, 42, 7, 3000, -1, 12], dtype=np.int64), 30) + ) + baseline = raw_vcov(X, resid, cl) + for _ in range(20): + np.testing.assert_array_equal(raw_vcov(X, resid, cl), baseline)