perf(linalg): scores-based one-way Bell-McCaffrey DOF denominator

CHANGELOG.md

@@ -623,6 +623,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   full fit. Uncertified LSMR (istop outside {0,1,2,4,5} after an uncapped retry) fails
   closed: NaN vcov/SE on the analytical boundary, the established `None` degenerate on
   the bootstrap boundary.
+- **One-way Bell-McCaffrey DOF: scores-based denominator (no dense `n×n` residual-maker).**
+  The non-clustered unweighted BM DOF computed `a'(M∘M)a` against an explicit dense
+  `M = I − H` (`O(n²k)` hat build; ~20 GB at n=50k, documented "practical for n < 10k").
+  The Schur-product expansion `Σ a_i²(1−2h_ii) + tr((B S_a)²)` with `S_a = X'diag(a)X`
+  (k×k) is exact algebra: `O(n k² + k³)` per contrast, n=50k/k=20 in ~16 ms. Frozen
+  dense-oracle parity at rtol 1e-10 (incl. high-leverage and k=40 designs + a compound
+  contrast); R-anchored hc2_bm goldens pass unchanged. New noise-floor cancellation
+  guard NaNs extreme-leverage contrasts whose expanded denominator collapses below
+  float precision instead of reporting an arbitrarily inflated DOF (the same failure
+  mode the clustered scores path guards; the prior dense `den > 0` kept such
+  denominators). Completes the scores-evaluation family: clustered CR2-BM (#656),
+  low-rank A_g (#664), and now the one-way path.
 - **Rust-backend HC2 vcov.** The Rust vcov path supported only HC1/CR1; one-way
   (unclustered, unweighted) HC2 now dispatches to a new `compute_robust_vcov_hc2` kernel
   mirroring the NumPy branch exactly — hat diagonals off the same bread,

TODO.md

@@ -43,7 +43,6 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 |-------|----------|--------|--------|----------|
 | `EfficientDiD` conditional path: the largest remaining O(n) stage is the sieve/nuisance construction outside the tiled pass (~9s at 10k). (The `_ridge_solve_weights` Python-prep shave landed 2026-07-07 — the `omega_stack[rest]` fancy-index copy and tail scatter are skipped when no row is zero-masked, byte-identical outputs; the `zero_mask` abs scan itself remains, needed for correctness.) | `efficient_did_covariates.py` | CS-scaling | Mid | Low |
 | Rust `solve_ols` always runs a full thin SVD (equilibrated, gelsd-parity); at 40 covariates it is the top per-cell solver item on a CS dr fit (1.04s of 3.1s at 2M rows, 95 calls). A Cholesky/QR fast path with SVD fallback — mirroring the Phase 3 Python-side pattern (certify well-conditioned, fall back verbatim) — is the natural lever, but `solve_ols` is the universal OLS entry point (every estimator), so the blast radius needs the full backend-parity treatment (`TestSolveOLSSkipRankCheckParity` posture: fitted-values parity, not beta). | `rust/src/linalg.rs::solve_ols` | CS-scaling | Heavy | Low |
-| One-way (non-clustered) unweighted Bell-McCaffrey DOF still materializes the dense `n×n` residual-maker `M` in `_compute_bm_dof_from_contrasts` (`O(n²k)` hat build + `O(n²m)` per-contrast sums; practical for n < 10k). Switch to a scores-based evaluation like the clustered CR2-BM path now uses (singleton-cluster reduction of the same `B = diag(‖ω‖²) − P'M_U P` identity, row-chunked). | `diff_diff/linalg.py::_compute_bm_dof_from_contrasts` | #656 | Mid | Low |
 
 ### Testing / docs
 

diff_diff/linalg.py

@@ -2746,12 +2746,13 @@ def _compute_bm_dof_from_contrasts(
 
     where ``q = X (X'X)^{-1} c``, ``M = I - H``, ``a(i) = q(i)^2 / (1 - h_ii)``.
     Using the idempotent identity ``M^2 = M``, ``trace(B) = sum_i q(i)^2``
-    matches the numerator. Allocates an ``(n, n)`` temporary for ``M`` so the
-    cost is ``O(n^2 k)`` for the hat build plus ``O(n^2 m)`` for the per-
-    contrast sums. Practical for ``n < 10_000``; larger designs should switch
-    to a scores-based formulation like the clustered CR2-BM path
-    (`_cr2_bm_dof_inner`) now uses — tracked as its own TODO.md row (the
-    original CR2-BM row this note pointed at is resolved).
+    matches the numerator. The denominator ``a'(M∘M)a`` is evaluated via the
+    Schur-product expansion (see the inline derivation) at ``O(n k^2 + k^3)``
+    per contrast with NO dense ``n×n`` residual-maker — the prior form's
+    ``O(n^2 k)`` hat build limited it to ``n < 10_000``. A noise-floor
+    cancellation guard NaNs extreme-leverage contrasts whose expanded
+    denominator collapses below float precision (mirrors the clustered
+    scores path's guard).
 
     **Weighted** (``weights is not None``): dispatches to the clubSandwich
     singleton-cluster CR2 reduction (each observation is its own cluster)
@@ -2779,7 +2780,8 @@ def _compute_bm_dof_from_contrasts(
     Returns
     -------
     ndarray of shape (m,) of Satterthwaite DOF per contrast column. NaN when
-    ``den <= 0`` (degenerate case).
+    the denominator is non-positive or at/below the cancellation noise
+    floor (degenerate / extreme-leverage case; see the inline guard note).
     """
     n, k = X.shape
     if contrasts.ndim != 2 or contrasts.shape[0] != k:
@@ -2811,18 +2813,41 @@ def _compute_bm_dof_from_contrasts(
         raise
     # q has shape (n, m); column j is X @ (bread_inv @ contrasts[:, j]).
     q = X @ bread_inv_c
-    H = X @ np.linalg.solve(bread_matrix, X.T)
-    M = np.eye(n) - H
-    M_sq = M * M  # elementwise square
     one_minus_h = np.maximum(1.0 - h_diag, 1e-10)
+    one_minus_2h = 1.0 - 2.0 * h_diag
     m = contrasts.shape[1]
     dof = np.empty(m)
     for j in range(m):
         qj_sq = q[:, j] * q[:, j]
         num = qj_sq.sum() ** 2
         a_j = qj_sq / one_minus_h
-        den = float(a_j @ M_sq @ a_j)
-        dof[j] = num / den if den > 0 else np.nan
+        # SCORES-BASED EVALUATION (2026-07, evaluation change — exact
+        # algebra): den = a'(M∘M)a with M = I − XBX' expands via the Schur
+        # product as
+        #
+        #   den = sum_i a_i^2 (1 − 2 h_ii) + tr((B S_a)^2),
+        #   S_a = X' diag(a) X   (k x k),
+        #
+        # because (M∘M)_{il} = δ_{il}(1 − 2 h_ii) + (x_i'B x_l)^2 and
+        # sum_{i,l} a_i a_l (x_i'B x_l)^2 = tr(B S_a B S_a). O(n k^2 + k^3)
+        # per contrast, never materializing the dense n×n residual-maker
+        # (the prior form's O(n^2 k) hat build + O(n^2) per contrast).
+        term_diag = float(np.sum(a_j * a_j * one_minus_2h))
+        S_a = X.T @ (X * a_j[:, np.newaxis])
+        BS = np.linalg.solve(bread_matrix, S_a)
+        term_tr = float(np.sum(BS * BS.T))
+        den = term_diag + term_tr
+        # Cancellation guard: the dense den = a'(M∘M)a is >= 0 exactly
+        # (M∘M is PSD — Schur product of the PSD M with itself), but the
+        # expanded difference of same-magnitude terms can collapse to the
+        # float64 noise floor for extreme-leverage contrasts. A den at or
+        # below the floor would inflate the DOF arbitrarily (the same
+        # failure mode the clustered scores path NaN-guards); report NaN
+        # instead of a non-physical DOF. The prior dense path's `den > 0`
+        # kept such noise-floor denominators; on well-conditioned contrasts
+        # the two evaluations agree to ~1e-12 (frozen-oracle parity test).
+        noise_floor = np.finfo(float).eps * (abs(term_diag) + abs(term_tr)) * n
+        dof[j] = num / den if den > noise_floor else np.nan
     return dof
 
 

docs/methodology/REGISTRY.md

@@ -244,7 +244,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. 11; equivalently the `q=1` reduction of the AHT Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast; peak memory = two `O(n k)` input-scale score precomputes plus working buffers capped by a byte cap subject to a one-contrast lower bound (a single contrast intrinsically needs `O(n)` + `O(G k)` buffers) — q vectors, per-cluster omegas, and product buffers contrast-chunked with every width-scaled buffer counted in the chunk denominator, and the `(G, G)` pairwise matrix row-chunked (Frobenius sum and max are row-separable), so none of the `n×n` residual-maker, `O(n m)` score arrays, or `O(G²)` pairwise entries is ever held at once; chunk-count invariant to ~1 ULP — BLAS width-dependent column accumulation; algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test). **Low-rank A_g factorization (2026-07):** the unweighted per-cluster adjustment operator `A_g = (I − H_gg)^{−1/2}` is evaluated from the k×k eigenproblem of `U_g'U_g` (`U_g = X_g M_U^{1/2}`; `A_g = I + (U_g Q) diag(γ) (U_g Q)'`, `γ = ((1−λ)^{−1/2}−1)/λ` via expm1/log1p, Moore-Penrose zeroing at `1−λ ≤ 1e-10` matching `_cr2_adjustment_matrix`'s convention) and only ever applied to skinny matrices — the dense `(n_g, n_g)` `A_g` (an `O(n_g³)` eigh per cluster) is never materialized for `n_g > k`; clusters with `n_g ≤ k` keep the (tiny) dense construction since it is the smaller Gram side. Algebraically identical; vcov/DOF match a frozen dense-eigh oracle at rtol 1e-12/1e-10 incl. the leverage-1 absorbed-cluster-FE and singleton-cluster regimes (`TestCR2BMLowRankAdjustment`). The weighted clubSandwich path keeps the dense construction (its `G_g` carries the `S_W` bias term).
+- **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. 11; equivalently the `q=1` reduction of the AHT Eq. 13) — instead of contracting against an explicit dense `n×n` residual-maker (`O(n k + G² k)` per contrast; peak memory = two `O(n k)` input-scale score precomputes plus working buffers capped by a byte cap subject to a one-contrast lower bound (a single contrast intrinsically needs `O(n)` + `O(G k)` buffers) — q vectors, per-cluster omegas, and product buffers contrast-chunked with every width-scaled buffer counted in the chunk denominator, and the `(G, G)` pairwise matrix row-chunked (Frobenius sum and max are row-separable), so none of the `n×n` residual-maker, `O(n m)` score arrays, or `O(G²)` pairwise entries is ever held at once; chunk-count invariant to ~1 ULP — BLAS width-dependent column accumulation; algebraically identical, ~1e-12 float agreement locked by a frozen-oracle parity test). **One-way scores evaluation (2026-07):** the non-clustered unweighted BM DOF denominator `a'(M∘M)a` (`_compute_bm_dof_from_contrasts`) is likewise evaluated without the dense `n×n` residual-maker, via the Schur-product expansion `Σ a_i²(1−2h_ii) + tr((B S_a)²)` with `S_a = X'diag(a)X` — exact algebra, `O(n k² + k³)` per contrast (was `O(n²k)`), frozen-dense-oracle parity ~1e-12, with a noise-floor cancellation guard NaN-ing extreme-leverage contrasts whose expanded denominator collapses below float precision (`TestOneWayBMScoresDOF`). **Low-rank A_g factorization (2026-07):** the unweighted per-cluster adjustment operator `A_g = (I − H_gg)^{−1/2}` is evaluated from the k×k eigenproblem of `U_g'U_g` (`U_g = X_g M_U^{1/2}`; `A_g = I + (U_g Q) diag(γ) (U_g Q)'`, `γ = ((1−λ)^{−1/2}−1)/λ` via expm1/log1p, Moore-Penrose zeroing at `1−λ ≤ 1e-10` matching `_cr2_adjustment_matrix`'s convention) and only ever applied to skinny matrices — the dense `(n_g, n_g)` `A_g` (an `O(n_g³)` eigh per cluster) is never materialized for `n_g > k`; clusters with `n_g ≤ k` keep the (tiny) dense construction since it is the smaller Gram side. Algebraically identical; vcov/DOF match a frozen dense-eigh oracle at rtol 1e-12/1e-10 incl. the leverage-1 absorbed-cluster-FE and singleton-cluster regimes (`TestCR2BMLowRankAdjustment`). The weighted clubSandwich path keeps the dense construction (its `G_g` carries the `S_W` bias term).
 - **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

@@ -2716,6 +2716,82 @@ def test_return_dropped_mask(self):
         assert nd0 == 3 and dropped0.all()
 
 
+class TestOneWayBMScoresDOF:
+    """The one-way (non-clustered) unweighted Bell-McCaffrey DOF denominator
+    a'(M∘M)a is evaluated via the Schur-product expansion
+    sum a_i^2(1-2h_ii) + tr((B S_a)^2), S_a = X'diag(a)X — never
+    materializing the dense n×n residual-maker (frozen here as the oracle).
+    Exact algebra; parity ~1e-12."""
+
+    @staticmethod
+    def _oracle_dense_dof(X, bread, h_diag, contrasts):
+        """Frozen pre-change dense evaluation (O(n²) M∘M quadratic form)."""
+        n = X.shape[0]
+        bread_inv_c = np.linalg.solve(bread, contrasts)
+        q = X @ bread_inv_c
+        H = X @ np.linalg.solve(bread, X.T)
+        M = np.eye(n) - H
+        M_sq = M * M
+        one_minus_h = np.maximum(1.0 - h_diag, 1e-10)
+        m = contrasts.shape[1]
+        dof = np.empty(m)
+        for j in range(m):
+            qj_sq = q[:, j] * q[:, j]
+            num = qj_sq.sum() ** 2
+            a_j = qj_sq / one_minus_h
+            den = float(a_j @ M_sq @ a_j)
+            dof[j] = num / den if den > 0 else np.nan
+        return dof
+
+    @pytest.mark.parametrize(
+        "kw",
+        [
+            dict(n=300, k=5, seed=3),
+            dict(n=200, k=8, seed=5, high_leverage=True),
+            dict(n=400, k=40, seed=7),
+        ],
+        ids=["basic", "high-leverage", "k40"],
+    )
+    def test_matches_frozen_dense_oracle(self, kw):
+        from diff_diff.linalg import _compute_bm_dof_from_contrasts
+
+        rng = np.random.default_rng(kw["seed"])
+        n, k = kw["n"], kw["k"]
+        X = np.column_stack([np.ones(n), rng.normal(size=(n, k - 1))])
+        if kw.get("high_leverage"):
+            X[0, 1] = 30.0
+        bread = X.T @ X
+        h_diag = np.einsum("ij,ij->i", X @ np.linalg.inv(bread), X)
+        contrasts = np.column_stack([np.eye(k), np.full(k, 1.0 / k)])  # + compound
+        dof = _compute_bm_dof_from_contrasts(X, bread, h_diag, contrasts)
+        oracle = self._oracle_dense_dof(X, bread, h_diag, contrasts)
+        fin = ~np.isnan(oracle)
+        assert fin.all(), "oracle produced NaN on a well-conditioned design"
+        np.testing.assert_allclose(dof[fin], oracle[fin], rtol=1e-10)
+
+    def test_noise_floor_guard_nans_leverage_one_contrast(self):
+        """A dummy column firing on exactly one observation gives that row
+        leverage 1: for the dummy's own coefficient the expanded
+        denominator's two terms cancel at ~1e20 scale down to the float
+        noise floor, so the guard must NaN it (the prior dense den > 0
+        would have kept the noise and inflated the DOF) while every
+        ordinary contrast in the same design stays finite."""
+        from diff_diff.linalg import _compute_bm_dof_from_contrasts
+
+        rng = np.random.default_rng(9)
+        n, k = 120, 4
+        X = np.column_stack([np.ones(n), rng.normal(size=(n, k - 2)), np.zeros(n)])
+        X[0, 3] = 1.0  # single-observation dummy -> h_00 = 1
+        bread = X.T @ X
+        h_diag = np.einsum("ij,ij->i", X @ np.linalg.pinv(bread), X)
+        assert h_diag.max() > 1 - 1e-12
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            dof = _compute_bm_dof_from_contrasts(X, bread, h_diag, np.eye(k))
+        assert np.isnan(dof[3]), "leverage-1 dummy coefficient must NaN"
+        assert np.isfinite(dof[:3]).all(), "ordinary contrasts must stay finite"
+
+
 class TestCR2BMLowRankAdjustment:
     """The low-rank factored A_g apply reproduces the dense per-cluster
     eigendecomposition path (frozen here as the oracle) at ~1e-14: vcov,