test+docs: chunk-count invariance is ~1 ULP, not bit-for-bit (CI failure fix)

CI (linux-arm + Windows) failed test_contrast_chunking_bit_identical
at 1-ULP differences: each contrast's B is computed independently, but
the per-cluster GEMM X_g' omega_g[:, c0:c1] runs over a width-c slice
and BLAS kernels (GEMV vs GEMM, platform-dependent) may accumulate a
column differently at width 1 vs width m — exact on Accelerate,
1-ULP drift elsewhere. The documented chunking-reassociation caveat
applies: test renamed to test_contrast_chunking_invariant with
assert_allclose(rtol=1e-13), and the bit-for-bit claims corrected in
the module-constant comment, docstring, CHANGELOG, and REGISTRY.

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

Author: igerber

CHANGELOG.md

@@ -562,7 +562,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   contrast, `O(n²)` memory — 3.2 GB at n=20k). The pairwise matrix now collapses to
   `B = diag(‖ω_g‖²) − P'M_U P` with `P = X'Ω` (disjoint cluster supports), costing
   `O(nk + G²k)` per contrast with no `n×n` allocation (per-cluster product buffers
-  contrast-chunked under a 64 MB cap, bit-identical across chunk counts, so the batched
+  contrast-chunked under a 64 MB cap, chunk-count invariant to ~1 ULP — BLAS kernels
+  can accumulate a GEMM column differently at different slice widths — so the batched
   per-coefficient sweep never allocates `O(G·k²)` at once): ~32x at n=5k/G=50 (0.57s→0.018s);
   n=20k/G=100 completes in 0.12s where the old path would allocate 3.2 GB. Algebraically
   identical — agreement with a frozen pair-loop oracle at rtol 1e-10 (balanced,

diff_diff/linalg.py

@@ -2153,8 +2153,12 @@ def _compute_cr2_bm(
 # (G, k, chunk) per-cluster product buffer so a batched per-coefficient sweep
 # (contrasts=eye(k), i.e. m == k) cannot allocate O(G*k*m) at once on
 # full-dummy / absorbed-FE designs with many clusters and coefficients.
-# Results are chunk-count invariant BIT-FOR-BIT: each contrast's B matrix is
-# computed independently, so chunking never reassociates a contrast's sums.
+# Each contrast's B matrix is computed independently, so chunking never
+# reassociates a contrast's OWN sums — but the per-cluster GEMM
+# `X_g' omega_g[:, c0:c1]` runs over a width-c slice, and BLAS kernels
+# (GEMV vs GEMM, platform-dependent) may accumulate a column differently at
+# different widths: chunk-count invariance holds to ~1 ULP (observed exact
+# on Accelerate, 1-ULP drift on OpenBLAS/arm + Windows), NOT bit-for-bit.
 # Module-level so tests can monkeypatch it to force the multi-chunk path.
 _CR2_BM_CONTRAST_CHUNK_BYTES = 64 * 1024 * 1024
 
@@ -2197,7 +2201,9 @@ def _cr2_bm_dof_inner(
 
     so ``trace_B2 = ||B||_F^2`` costs ``O(n k + G^2 k)`` per contrast with
     ``O(G k)`` memory per contrast (contrast-chunked buffers bounded by
-    ``_CR2_BM_CONTRAST_CHUNK_BYTES``, bit-identical across chunk counts) — the previous form materialized the dense ``n x n``
+    ``_CR2_BM_CONTRAST_CHUNK_BYTES``, chunk-count invariant to ~1 ULP — BLAS
+    kernels may accumulate a GEMM column differently at different slice
+    widths) — the previous form materialized the dense ``n x n``
     ``M`` and looped cluster pairs at ``O(n^2)`` per contrast, the exact
     large-``n`` blowup the TODO row tracked. The two evaluations are
     algebraically identical; floating-point agreement is ~1e-12 relative
@@ -2233,7 +2239,8 @@ def _cr2_bm_dof_inner(
     # Chunk the contrasts so the per-cluster product buffer is (G, k, c) with
     # c bounded by _CR2_BM_CONTRAST_CHUNK_BYTES — a full-m buffer would be
     # O(G*k*m), i.e. O(G*k^2) on the batched per-coefficient sweep. Each
-    # contrast's B is computed independently, so chunking is bit-identical.
+    # contrast's B is computed independently; chunk-count invariance holds
+    # to ~1 ULP (BLAS width-dependent column accumulation), not bit-for-bit.
     chunk = max(1, int(_CR2_BM_CONTRAST_CHUNK_BYTES // max(n_g_clusters * k_X * 8, 1)))
     for c0 in range(0, m, chunk):
         c1 = min(c0 + chunk, m)

docs/methodology/REGISTRY.md

@@ -243,7 +243,7 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients.
   `_compute_cr2_bm_contrast_dof` (DOF-only for arbitrary contrasts) are thin wrappers
   over that shared core, so every CR2 caller routes through one implementation. The
   consolidation is bit-identical to the prior two-call path (proven at atol=0).
-- **Note (unweighted per-coef DOF guard):** the unweighted clustered CR2-BM per-coefficient DOF (`_cr2_bm_dof_inner`, the simple `(tr B)²/tr(B²)` form) carries the same two-part reliability guard as the weighted P-array path. (1) **Noise floor:** for a high-leverage FE-dummy / collinear nuisance column `trace_B2 = Σ B_{g,h}²` collapses to float64 accumulation noise while `trace_B` stays O(1), inflating the ratio to a non-physical DOF (observed ~1e61 on the absorbed-FE golden); a contrast whose `max|B_{g,h}|` sits below the batch-relative (`1e-10×max`, computed on the scale-normalized `max|B|/‖c‖²` since `B ∝ ‖c‖²` while the DOF is scale-invariant) or absolute (`(EPS·n·k·bread_scale)²`) floor is NaN'd. (2) **Cluster-count bound:** the Bell-McCaffrey Satterthwaite DOF is `(tr B)²/tr(B²)` with `B` PSD and cluster-structured, so it is bounded by `rank(B) ≤ G` (number of clusters); the simple unweighted form is numerically less faithful than clubSandwich's P-array form on high-leverage columns and can return a finite-but-inflated DOF above `G` (observed ~32.7 and ~16.3 vs R's 6 and 3, `G=8`), which is NaN'd as non-physical. The well-conditioned contrasts estimators consume — the treatment effect, event-study coefficients, and the compound post-period-average — are unaffected and match R clubSandwich; only the non-user-facing high-leverage nuisance DOFs are suppressed (exact P-array reproduction of those is deferred). A `UserWarning` fires per fit. Regression: `tests/test_estimators_vcov_type.py::TestDiDAbsorbedFERParity::test_unweighted_cr2_bm_per_coef_dof_no_nonphysical`. **Evaluation (2026-07):** the pairwise `B` matrix is computed via the algebraic identity `B = Ω'MΩ = diag(‖ω_g‖²) − P' M_U P` with `P = X'Ω` (the Satterthwaite DOF itself is Pustejovsky-Tipton 2018 §3.1 / Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast, no `n×n` allocation, per-cluster product buffers contrast-chunked under a byte cap (bit-identical across chunk counts); algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test).
+- **Note (unweighted per-coef DOF guard):** the unweighted clustered CR2-BM per-coefficient DOF (`_cr2_bm_dof_inner`, the simple `(tr B)²/tr(B²)` form) carries the same two-part reliability guard as the weighted P-array path. (1) **Noise floor:** for a high-leverage FE-dummy / collinear nuisance column `trace_B2 = Σ B_{g,h}²` collapses to float64 accumulation noise while `trace_B` stays O(1), inflating the ratio to a non-physical DOF (observed ~1e61 on the absorbed-FE golden); a contrast whose `max|B_{g,h}|` sits below the batch-relative (`1e-10×max`, computed on the scale-normalized `max|B|/‖c‖²` since `B ∝ ‖c‖²` while the DOF is scale-invariant) or absolute (`(EPS·n·k·bread_scale)²`) floor is NaN'd. (2) **Cluster-count bound:** the Bell-McCaffrey Satterthwaite DOF is `(tr B)²/tr(B²)` with `B` PSD and cluster-structured, so it is bounded by `rank(B) ≤ G` (number of clusters); the simple unweighted form is numerically less faithful than clubSandwich's P-array form on high-leverage columns and can return a finite-but-inflated DOF above `G` (observed ~32.7 and ~16.3 vs R's 6 and 3, `G=8`), which is NaN'd as non-physical. The well-conditioned contrasts estimators consume — the treatment effect, event-study coefficients, and the compound post-period-average — are unaffected and match R clubSandwich; only the non-user-facing high-leverage nuisance DOFs are suppressed (exact P-array reproduction of those is deferred). A `UserWarning` fires per fit. Regression: `tests/test_estimators_vcov_type.py::TestDiDAbsorbedFERParity::test_unweighted_cr2_bm_per_coef_dof_no_nonphysical`. **Evaluation (2026-07):** the pairwise `B` matrix is computed via the algebraic identity `B = Ω'MΩ = diag(‖ω_g‖²) − P' M_U P` with `P = X'Ω` (the Satterthwaite DOF itself is Pustejovsky-Tipton 2018 §3.1 / Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast, no `n×n` allocation, per-cluster product buffers contrast-chunked under a byte cap (chunk-count invariant to ~1 ULP — BLAS width-dependent column accumulation); algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test).
 - **Note:** `LinearRegression.get_se()` / `get_inference()` clamp the vcov diagonal at 0
   before `sqrt`. A high-leverage / degenerate coefficient (an absorbed-FE dummy
   near-collinear with the treatment, whose Satterthwaite DOF already hits the noise-floor

tests/test_linalg.py

@@ -2808,11 +2808,15 @@ def test_unbalanced_clusters_match_oracle(self):
         oracle = self._oracle_pairloop_dof(X, cl, bread, np.eye(k))
         np.testing.assert_allclose(dof, oracle, rtol=1e-10)
 
-    def test_contrast_chunking_bit_identical(self, monkeypatch):
+    def test_contrast_chunking_invariant(self, monkeypatch):
         """CI-review P2: the per-cluster product buffer is contrast-chunked
         (bounded by _CR2_BM_CONTRAST_CHUNK_BYTES) instead of O(G*k*m).
-        Forcing one-contrast chunks must reproduce the single-chunk DOF
-        bit-for-bit — each contrast's B is computed independently."""
+        Forcing one-contrast chunks reproduces the single-chunk DOF to
+        ~1 ULP: each contrast's B is computed independently, but the
+        per-cluster GEMM runs over a width-c slice and BLAS kernels may
+        accumulate a column differently at width 1 vs width m (observed
+        exact on Accelerate, 1-ULP drift on OpenBLAS/arm + Windows CI —
+        the documented chunking-reassociation caveat)."""
         import diff_diff.linalg as la
 
         rng = np.random.default_rng(23)
@@ -2826,4 +2830,4 @@ def test_contrast_chunking_bit_identical(self, monkeypatch):
         _, dof_one = la._compute_cr2_bm(X, resid, cl, bread)
         monkeypatch.setattr(la, "_CR2_BM_CONTRAST_CHUNK_BYTES", G * k * 8)  # 1 contrast/chunk
         _, dof_many = la._compute_cr2_bm(X, resid, cl, bread)
-        np.testing.assert_array_equal(dof_many, dof_one)
+        np.testing.assert_allclose(dof_many, dof_one, rtol=1e-13)