diff --git a/CHANGELOG.md b/CHANGELOG.md index 2ec9a5e3..407f59e0 100644 --- a/CHANGELOG.md +++ b/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, diff --git a/TODO.md b/TODO.md index 3bc18b15..bf3ba813 100644 --- a/TODO.md +++ b/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 --git a/diff_diff/linalg.py b/diff_diff/linalg.py index 9af39a1a..28b75a00 100644 --- a/diff_diff/linalg.py +++ b/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 diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 90d0074d..e6a77e7c 100644 --- a/docs/methodology/REGISTRY.md +++ b/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 diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 1f48eb20..f91738ac 100644 --- a/tests/test_linalg.py +++ b/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,