linalg: address CI codex R8 P2 (single-contrast noise-floor guard)

CI codex on PR #475 (✅ verdict) flagged a real P2: the noise-floor NaN-
guard in `_cr2_bm_dof_inner_weighted` was batch-relative only — for a
single-contrast call to `_compute_cr2_bm_contrast_dof`, `max|P|_overall`
equals the contrast's own max|P|, so the `1e-10 × max|P|_overall` rule
could never classify it as degenerate. That left direct single-contrast
weighted callers (e.g., MPD avg_att) unprotected: they could still emit
BLAS-implementation-dependent finite DOF on noise-floor contrasts even
though the registry/changelog said the helper was guarded.

Fix: union the batch-relative criterion with an absolute floor scaled to
the bread matrix's magnitude: `(EPS × n × k × max(bread_inv_scale, 1))²`.
This covers the worst-case dgemm accumulation roundoff floor for
`H1/H2/H3 @ contrast` products. A single-contrast call now correctly
fires the NaN-guard on a high-leverage FE-dummy contrast.

New regression tests in `tests/test_methodology_wls_cr2.py::
TestWLSCR2SingleContrastNoiseFloor` (2 tests): single weighted FE-dummy
contrast triggers NaN-guard + warning; single non-noise contrast still
returns finite DOF matching clubSandwich at atol=1e-10.

CI codex P3 (perf): LinearRegression.fit() pays CR2 twice on the new
weighted hc2_bm path (solve_ols + compute_robust_vcov). Added as a TODO
follow-up row (PR #475 follow-up, Low priority).

All 198 linalg+methodology+estimators-vcov-type+TWFE+SA tests pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: igerber

TODO.md

@@ -108,6 +108,7 @@ Deferred items from PR reviews that were not addressed before merge.
      scenarios in benchmarks/data/clubsandwich_cr2_golden.json. -->
 | ~~Weighted one-way Bell-McCaffrey (`vcov_type="hc2_bm"` + `weights`, no cluster)~~ — LIFTED via clubSandwich WLS-CR2 port. Routes through `_compute_cr2_bm_contrast_dof` with singleton-cluster reduction (each obs is its own cluster). See REGISTRY.md Phase 1a row for the algebra. | `linalg.py` | LIFTED | — |
 | ~~Weighted CR2 Bell-McCaffrey cluster-robust (`vcov_type="hc2_bm"` + `cluster_ids` + `weights`)~~ — LIFTED via clubSandwich WLS-CR2 port. `_compute_cr2_bm` now accepts `weights=` and threads through clubSandwich's specific WLS-CR2 algebra. | `linalg.py::_compute_cr2_bm` | LIFTED | — |
+| `LinearRegression.fit()` pays the CR2 cost twice on the weighted `hc2_bm` path: once inside `solve_ols(..., return_vcov=True)` and again via `compute_robust_vcov(..., return_dof=True)` to populate `_bm_dof`. Correct but redundant. Fix: thread `return_dof` through `solve_ols` so the same CR2 computation produces both vcov + DOF, or cache the per-cluster `A_g` / `MUWTWUM` precomputes between calls. CI codex P3 on PR #475. | `linalg.py::LinearRegression.fit`, `linalg.py::solve_ols` | PR #475 follow-up | Low |
 | `TwoWayFixedEffects(vcov_type in {"hc2","hc2_bm"})` with replicate-weight survey designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate (re-demeaning depends on the per-replicate weight vector), which doesn't compose with the full-dummy HC2/HC2-BM build — a correct implementation would need per-replicate full-dummy refit. Workaround: use `vcov_type="hc1"` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Low |
 | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper (or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch, with TWFE-specific cluster default threading) to reduce drift risk on FE naming, survey behavior, and result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Low |
 | Unify Rust local-method `estimate_model` solver path to `solve_wls_svd` (the same SVD helper used by the global-method since PR #348) for sub-1e-14 bootstrap SE parity. Current local-method bootstrap parity test (`tests/test_rust_backend.py::TestTROPRustEdgeCaseParity::test_bootstrap_seed_reproducibility_local`) passes at `atol=1e-5` — the residual ~1e-7 gap is roundoff between Rust's `estimate_model` matrix factorization and numpy's `lstsq`, which accumulates differently across per-replicate bootstrap fits. Main-fit ATT parity is regime-dependent (`atol=1e-14` for `lambda_nn=inf`, `atol=1e-10` for finite `lambda_nn` — see `test_local_method_main_fit_parity`); the bootstrap gap is a same-solver-path roundoff concern and not a user-visible correctness bug. | `rust/src/trop.rs::estimate_model`, `rust/src/linalg.rs::solve_wls_svd` | follow-up | Low |

diff_diff/linalg.py

@@ -1859,27 +1859,46 @@ def _factor_psd(M):
     # different specific values in this regime due to its own BLAS reduction
     # order; we return NaN with a warning rather than ship arbitrarily-
     # different small-sample DOF on the user-facing surface.
+    # The noise-floor detector has two criteria, applied union-wise:
+    #   1. Batch-relative: a contrast's max(|P|) < 1e-10 × the largest contrast's
+    #      max(|P|). Useful when computing per-coefficient DOF for an entire design
+    #      (`contrasts=np.eye(k)`) — a single noise-floor coefficient stands out.
+    #   2. Absolute (single-contrast safe): max(|P|) < eps_floor based on the
+    #      bread matrix scale. Necessary because batch-relative reduces to
+    #      max|P| < 1e-10 × max|P| (i.e. False) for a single contrast, leaving
+    #      direct `_compute_cr2_bm_contrast_dof(...)` callers (e.g., MPD avg_att)
+    #      unprotected. Threshold: `(EPS × n × k × bread_inv_scale)²` covers the
+    #      worst-case dgemm accumulation roundoff floor for `H1/H2/H3 @ contrast`
+    #      products. CI codex flagged this as P2 (R8 round); regression test in
+    #      tests/test_methodology_wls_cr2.py::TestWLSCR2FEDoFNoiseGuard.
+    _EPS = np.finfo(np.float64).eps
+    n_obs, k_X = X.shape
+    bread_inv_scale = float(np.max(np.abs(bread_inv))) if bread_inv.size else 1.0
+    abs_noise_floor = (_EPS * n_obs * k_X * max(bread_inv_scale, 1.0)) ** 2
+
     if m > 1:
         max_P_overall = float(np.max(max_abs_P_arr))
-        if max_P_overall > 0:
-            noise_floor = 1e-10 * max_P_overall
-            degenerate = max_abs_P_arr < noise_floor
-            n_degenerate = int(np.sum(degenerate))
-            if n_degenerate > 0:
-                dof_vec[degenerate] = np.nan
-                warnings.warn(
-                    f"Satterthwaite DOF for {n_degenerate} of {m} contrast(s) "
-                    f"is at the float64 noise floor (max|P| < 1e-10 × "
-                    f"max|P|_overall); reporting NaN. This typically affects "
-                    f"high-leverage FE-dummy coefficients whose contrast "
-                    f"vector projects to near-zero on the design — the "
-                    f"resulting DOF varies across BLAS implementations and is "
-                    f"unreliable. The coefficient SEs remain valid; only the "
-                    f"Satterthwaite DOF (and any t-test or CI that depends on "
-                    f"it) is suppressed.",
-                    UserWarning,
-                    stacklevel=3,
-                )
+        relative_floor = 1e-10 * max_P_overall if max_P_overall > 0 else 0.0
+    else:
+        relative_floor = 0.0
+    noise_floor = max(relative_floor, abs_noise_floor)
+    degenerate = max_abs_P_arr < noise_floor
+    n_degenerate = int(np.sum(degenerate))
+    if n_degenerate > 0:
+        dof_vec[degenerate] = np.nan
+        warnings.warn(
+            f"Satterthwaite DOF for {n_degenerate} of {m} contrast(s) "
+            f"is at the float64 noise floor (max|P| < noise_floor = "
+            f"max({relative_floor:.3e}, {abs_noise_floor:.3e})); reporting "
+            f"NaN. This typically affects high-leverage FE-dummy "
+            f"coefficients whose contrast vector projects to near-zero on "
+            f"the design — the resulting DOF varies across BLAS "
+            f"implementations and is unreliable. The coefficient SEs "
+            f"remain valid; only the Satterthwaite DOF (and any t-test or "
+            f"CI that depends on it) is suppressed.",
+            UserWarning,
+            stacklevel=3,
+        )
 
     return dof_vec
 

tests/test_methodology_wls_cr2.py

@@ -533,6 +533,96 @@ def test_did_absorbed_fe_nan_guards_treated_unit_dummies(self, goldens):
         ), f"Expected 1 noise-floor UserWarning, got {len(noise_warnings)}"
 
 
+class TestWLSCR2SingleContrastNoiseFloor:
+    """Regression for CI codex R8 P2: the noise-floor NaN-guard must trigger
+    on single-contrast calls to `_compute_cr2_bm_contrast_dof` too, not just
+    batched per-coefficient calls. The earlier batch-relative-only criterion
+    couldn't catch single contrasts because `max|P|_overall == contrast's own
+    scale` for m=1. The absolute (eps × n × k × bread_scale)² floor closes
+    the gap."""
+
+    def test_single_fe_dummy_contrast_triggers_nan_guard(self, goldens):
+        """A direct single-contrast call for a high-leverage FE-dummy unit
+        dummy must return NaN, not a BLAS-dependent finite DOF."""
+        import warnings
+
+        d = goldens["weighted_did_absorbed_fe"]
+        n = len(d["y"])
+        unit_int = np.asarray(d["unit"])
+        period_int = np.asarray(d["period"])
+        treat_post = np.asarray(d["treat_post"])
+        np.asarray(d["y"])
+        w = np.asarray(d["weights"])
+        cluster = unit_int
+
+        cols = [np.ones(n), treat_post.astype(float)]
+        col_names = ["(Intercept)", "treat_post"]
+        for u in sorted(np.unique(unit_int))[1:]:
+            cols.append((unit_int == u).astype(float))
+            col_names.append(f"unit{u}")
+        for p in sorted(np.unique(period_int))[1:]:
+            cols.append((period_int == p).astype(float))
+            col_names.append(f"period{p}")
+        X = np.column_stack(cols)
+        bread = X.T @ (X * w[:, None])
+
+        # Single contrast = unit2 dummy (the noise-floor case from the
+        # weighted_did_absorbed_fe fixture).
+        unit2_idx = col_names.index("unit2")
+        c_unit2 = np.zeros(X.shape[1])
+        c_unit2[unit2_idx] = 1.0
+        contrasts = c_unit2[:, None]  # (k, 1) — single contrast
+
+        with warnings.catch_warnings(record=True) as ws:
+            warnings.simplefilter("always")
+            dof = _compute_cr2_bm_contrast_dof(X, cluster, bread, contrasts, weights=w)
+        assert np.isnan(dof[0]), (
+            "Single weighted FE-dummy contrast must trigger absolute noise-"
+            "floor guard; got finite DOF instead."
+        )
+        noise_warnings = [w for w in ws if "noise floor" in str(w.message)]
+        assert (
+            len(noise_warnings) == 1
+        ), f"Expected 1 noise-floor UserWarning, got {len(noise_warnings)}"
+
+    def test_single_non_noise_contrast_returns_finite_dof(self, goldens):
+        """Sanity check: a single contrast that's NOT at the noise floor
+        (e.g., treat_post) must return finite DOF, not NaN."""
+        d = goldens["weighted_did_absorbed_fe"]
+        n = len(d["y"])
+        unit_int = np.asarray(d["unit"])
+        period_int = np.asarray(d["period"])
+        treat_post = np.asarray(d["treat_post"])
+        np.asarray(d["y"])
+        w = np.asarray(d["weights"])
+        cluster = unit_int
+
+        cols = [np.ones(n), treat_post.astype(float)]
+        col_names = ["(Intercept)", "treat_post"]
+        for u in sorted(np.unique(unit_int))[1:]:
+            cols.append((unit_int == u).astype(float))
+            col_names.append(f"unit{u}")
+        for p in sorted(np.unique(period_int))[1:]:
+            cols.append((period_int == p).astype(float))
+            col_names.append(f"period{p}")
+        X = np.column_stack(cols)
+        bread = X.T @ (X * w[:, None])
+
+        # Single contrast = treat_post (NOT a noise-floor contrast).
+        treat_post_idx = col_names.index("treat_post")
+        c_treat = np.zeros(X.shape[1])
+        c_treat[treat_post_idx] = 1.0
+        contrasts = c_treat[:, None]
+
+        dof = _compute_cr2_bm_contrast_dof(X, cluster, bread, contrasts, weights=w)
+        assert np.isfinite(dof[0]) and dof[0] > 0, (
+            f"Single non-noise contrast should return finite positive DOF; " f"got {dof[0]}"
+        )
+        # Sanity: should match R's df_satt for treat_post (3.754) at atol=1e-10.
+        r_dof_treat_post = float(np.asarray(d["dof_cr2"])[treat_post_idx])
+        np.testing.assert_allclose(dof[0], r_dof_treat_post, atol=1e-10)
+
+
 class TestLinearRegressionFENanGuardEndToEnd:
     """Regression for R7 P0: NaN BM DOF from the noise-floor guard must
     propagate to NaN inference fields in `LinearRegression.get_inference()`,