fix(linalg): rank-guard structural (non-covariate) matrix inverses

CHANGELOG.md

@@ -118,6 +118,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   numerical, or public-API change.
 
 ### Fixed
+- **Structural (non-covariate) matrix inverses are now rank-guarded.** The internal design-Gram
+  bread inversions in `ContinuousDiD` (ACRT-variance `Psi'WPsi`), `TwoStageDiD` (Stage-2
+  `X_2'WX_2`, both the analytical and multiplier-bootstrap surfaces), `SpilloverDiD` (Wave D
+  `A_22`), and the Conley spatial-HAC variance (`X'WX`) sat on a `LinAlgError`-only fallback:
+  `np.linalg.inv`/`solve` raise only on an *exactly* singular matrix, so a **near**-singular
+  Gram returned a garbage inverse (~1e13) straight into the SE (and `ContinuousDiD`'s exact-
+  singular fallback was a *silent* minimum-norm `pinv`; Conley *raised* `ValueError`). All now
+  route through the shared `_rank_guarded_inv` (`diff_diff/linalg.py`) — the same generalized
+  inverse already used for the covariate IF SEs — which truncates redundant directions on the
+  equilibrated Gram to give a finite SE on the identified subspace (the well-conditioned
+  near-collinear limit, not minimum-norm; NaN only at rank 0) and emits a `UserWarning` when a
+  direction is dropped. For the per-coefficient reporters (TwoStageDiD event/group, SpilloverDiD
+  rings, Conley), a dropped (unidentified) named coefficient is reported with **NaN** SE (not the
+  zero-filled `0`); linear-combination reporters (ContinuousDiD dose curves) are unaffected, since
+  a dropped direction correctly contributes 0 there. These invert *internal* bases users cannot
+  perturb with `covariates=`.
+  **Behavior change:** a rank-deficient Conley design no longer raises — it rank-reduces with a
+  warning. Well-conditioned designs are unchanged (the fast path is `np.linalg.solve(A, I)`,
+  R-parity preserved). (`HeterogeneousAdoptionDiD`'s non-symmetric IV bread and `ImputationDiD`
+  — whose vcov is already rank-guarded upstream via `solve_ols` — were assessed and excluded;
+  see `TODO.md`.)
 - **`CallawaySantAnna` / `StaggeredTripleDifference` covariate outcome-regression is now
   scale-robust.** The covariate OR nuisance fits — `_compute_all_att_gt_covariate_reg` and
   `_doubly_robust` (CS) and `_compute_or` (StaggeredTripleDifference) — previously used an

TODO.md

@@ -30,7 +30,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 |-------|----------|--------|--------|----------|
 | `SyntheticControl` cv: thread an `"infeasible"` reason-code from `_outer_solve_V_cv()` / `_placebo_fit_unit()` so `in_space_placebo()` / `leave_one_out()` distinguish a structural cv-refit exclusion (donor-indistinguishable re-aggregated window) from a genuine inner-solver non-convergence — mirror the split `in_time_placebo()` already emits. Warnings already distinguish the two causes; only the machine-readable status/count is missing. | `synthetic_control.py`, `synthetic_control_results.py` | follow-up | Mid | Low |
 | `CallawaySantAnna`: materialize NaN entries for non-estimable `(g,t)` cells in `group_time_effects` (currently omitted with a consolidated warning); requires updating downstream consumers (event study, `balance_e`, aggregation). | `staggered.py` | #256 | Mid | Low |
-| Adopt the shared `_rank_guarded_inv` for the *structural* (non-covariate) matrix inverses sharing the `LinAlgError`-only fallback that can go near-singular: `continuous_did.py:1056` (dose B-spline), `spillover.py:3371` (ring-solve, partially guarded), `two_stage.py:3154` (TSL Stage-2 variance), `imputation.py:2403`, `had.py:2413`, `conley.py:1109`. These invert internal bases users cannot perturb with `covariates=` (distinct from the already-fixed covariate-triggered SE bug); the helper is the seam. | `continuous_did.py`, `spillover.py`, `two_stage.py`, `imputation.py`, `had.py`, `conley.py` | dr-or-se-rank-guard | Mid | Low |
 | Survey-design resolution / collapse patterns are inconsistent across panel estimators — `ContinuousDiD` rebuilds unit-level design in SE code, `EfficientDiD` builds once in `fit()`, `StackedDiD` re-resolves on stacked data. Extract shared helpers for panel-to-unit collapse, post-filter re-resolution, and metadata recomputation. | `continuous_did.py`, `efficient_did.py`, `stacked_did.py` | #226 | Mid | Low |
 | `SyntheticControl` remaining ADH-2015 §4 items: the regression-weight `W^reg = X_0'(X_0 X_0')^{-1} X_1` extrapolation diagnostic (flag implied OLS weights outside `[0,1]`) and sparse-SC subset search (`l < J`, holding `V` fixed). LOO, in-time placebo, CV `V`-selection, and inverse-variance `V` have landed; these two are the deferred tail. | `synthetic_control.py`, `synthetic_control_results.py` | ADH-2015 | Mid | Low |
 | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
@@ -150,6 +149,8 @@ Doable in principle, but no current caller and/or explicitly out of paper scope.
 |----------|----------|----------|
 | **`bread_inv` reuse not bit-identically achievable.** "Factor `(X'WX)` once, reuse across HC2/HC2-BM" can't be done bit-identically (the bar for a pure perf refactor of the inference path). Internal bread ops solve against *different* RHS (`X.T`, `eye`, `meat`+`temp.T`, `contrasts`); only same-RHS results are bit-reusable. Measured: `lu_solve(lu_factor(A),B)` differs from `solve(A,B)` up to 6.4e-15; the `inv(A)@meat@inv(A)` sandwich differs up to 1.24e-14 — both nonzero and *below* the affected goldens' tolerances (1e-12/1e-10), so a broad reuse would silently shift SEs without tripping the suite. The one genuine bit-identical redundancy (a duplicated `solve(bread, X.T)` in the unweighted one-way `hc2_bm`+`return_dof` path) is dwarfed by that path's dense `M=I−H` build, so the saving is negligible. | `linalg.py::compute_robust_vcov` | 2026-06-01 |
 | **R-script-per-test consolidation has no CI impact.** No CI workflow installs R, so every R-parity test skips in CI behind a per-file availability gate — consolidating `Rscript` spawns yields zero CI speedup. `test_methodology_twfe.py` already session-caches its R fits. The only residual is a LOCAL-dev micro-opt for `test_methodology_continuous_did.py` / `test_methodology_callaway.py` (re-spawn `library(...)` per call). Low value; retained as a local-dev note. | `tests/test_methodology_continuous_did.py`, `tests/test_methodology_callaway.py` | #139 / 2026-06-07 |
+| **`HeterogeneousAdoptionDiD` mass-point IV bread is non-symmetric — `_rank_guarded_inv` inapplicable.** The structural rank-guard sweep (continuous_did / two_stage / spillover / conley) excluded `had.py`'s `ZtWX = Zd'WX` ([1, instrument]' × [1, endogenous]): it is a non-symmetric 2×2 Wald-IV bread (`V = ZtWX_inv @ Omega @ ZtWX_inv.T`), and `_rank_guarded_inv` assumes a **symmetric PSD** Gram (symmetric `D=diag(A)` equilibration + eigendecomposition), so applying it would be methodologically wrong. The existing fallback already returns NaN SE on a singular bread; an IV-appropriate near-singular guard would need a different mechanism. | `had.py` | structural-rank-guard / 2026-06-28 |
+| **`ImputationDiD` SE vcov is already rank-guarded upstream.** Excluded from the structural rank-guard sweep: the lead/effect vcov comes from `solve_ols(..., return_vcov=True, rank_deficient_action=...)` at the OLS fit (`imputation.py:~2316`), which already drops rank-deficient columns. The only raw inverse (`solve(V_gamma, gamma)`, `imputation.py:~2530`) is the pretrends **Wald F-test statistic** with a safe `NaN` fallback — a test statistic, not a sandwich bread — so there is no garbage-SE exposure. No structural rank-guard needed. | `imputation.py` | structural-rank-guard / 2026-06-28 |
 
 ---
 

diff_diff/conley.py

@@ -1104,16 +1104,37 @@ def _compute_conley_vcov(
         _conley_sparse=_conley_sparse,
     )
 
-    # Sandwich via two solves (mirrors _compute_cr2_bm pattern in linalg.py)
-    try:
-        temp = np.linalg.solve(bread_matrix, meat)
-        vcov = np.linalg.solve(bread_matrix, temp.T).T
-    except np.linalg.LinAlgError as e:
-        if "Singular" in str(e):
-            raise ValueError(
-                "Design matrix is rank-deficient (singular X'X matrix). "
-                "Cannot compute Conley spatial HAC variance."
-            ) from e
-        raise
+    # Sandwich via the shared rank-guarded inverse of the design Gram.
+    # np.linalg.solve only raises on an *exactly* singular bread, so a *near*-
+    # singular X'WX would otherwise flow a garbage inverse (~1e13) straight into
+    # the spatial-HAC variance. `_rank_guarded_inv` truncates redundant
+    # directions on the equilibrated Gram -> a finite SE on the identified
+    # subspace (NaN only at rank 0), matching the covariate IF rank-guard and the
+    # other structural bread inversions (ContinuousDiD / TwoStageDiD /
+    # SpilloverDiD). Lazy import: `linalg` imports this module, so a top-level
+    # `from diff_diff.linalg import ...` would be circular; resolving at call time
+    # is safe (linalg is already loaded by the time this runs).
+    from diff_diff.linalg import _rank_guarded_inv
+
+    bread_inv, n_dropped, _, dropped = _rank_guarded_inv(bread_matrix, return_dropped=True)
+    if n_dropped:
+        warnings.warn(
+            "Conley spatial HAC variance: the design Gram (X'WX) is "
+            f"rank-deficient ({n_dropped} redundant direction(s) dropped); "
+            "rank-reducing to a finite SE on the identified subspace "
+            "(NaN if rank 0). This usually indicates collinear regressors.",
+            UserWarning,
+            stacklevel=2,
+        )
+    # vcov = bread^{-1} @ meat @ bread^{-1}; algebraically identical to the prior
+    # two symmetric solves given `bread` symmetric (holds for any meat).
+    vcov = bread_inv @ meat @ bread_inv
+    # A dropped (unidentified) coefficient is zero-filled in bread_inv, which would
+    # otherwise report se=0 for that named coefficient. NaN its row/col in the
+    # FINAL vcov so per-coefficient SE extraction yields NaN (not 0) for the
+    # unidentified directions, while the identified coefficients stay finite.
+    if dropped.any():
+        vcov[dropped, :] = np.nan
+        vcov[:, dropped] = np.nan
 
     return vcov

diff_diff/continuous_did.py

@@ -29,7 +29,7 @@
     ContinuousDiDResults,
     DoseResponseCurve,
 )
-from diff_diff.linalg import solve_ols
+from diff_diff.linalg import _rank_guarded_inv, solve_ols
 from diff_diff.survey import (
     ResolvedSurveyDesign,
     _resolve_survey_for_fit,
@@ -1047,21 +1047,32 @@ def _compute_dose_response_gt(
 
         # Store bootstrap info for influence function computation
         # bread = (Psi'WPsi / n_treated)^{-1} when survey, (Psi'Psi / n_treated)^{-1} otherwise
+        # Bread = (Psi'WPsi / mass)^{-1} via the shared rank-guarded inverse:
+        # np.linalg.inv only raises on an *exactly* singular Gram, so a *near*-
+        # singular B-spline design (clustered doses / near-duplicate knots)
+        # previously returned a garbage inverse (~1e13) -> garbage SE. The prior
+        # `pinv` fallback was both minimum-norm (not the column-drop / near-
+        # collinear limit) and *silent*. `_rank_guarded_inv` truncates redundant
+        # directions on the equilibrated Gram -> finite SE on the identified
+        # subspace (NaN only at rank 0), matching the covariate IF rank-guard.
         if w_treated is not None:
             w_treated_sum = float(np.sum(w_treated))
             PtWP = Psi.T @ (Psi * w_treated[:, np.newaxis])
             # Normalize bread by weighted mass (not raw count) for consistency
             # with downstream IF score denominators that also use weighted mass
-            try:
-                bread = np.linalg.inv(PtWP / w_treated_sum)
-            except np.linalg.LinAlgError:
-                bread = np.linalg.pinv(PtWP / w_treated_sum)
+            bread, n_dropped, _ = _rank_guarded_inv(PtWP / w_treated_sum)
         else:
             PtP = Psi.T @ Psi
-            try:
-                bread = np.linalg.inv(PtP / n_treated)
-            except np.linalg.LinAlgError:
-                bread = np.linalg.pinv(PtP / n_treated)
+            bread, n_dropped, _ = _rank_guarded_inv(PtP / n_treated)
+        if n_dropped:
+            warnings.warn(
+                "ContinuousDiD ACRT variance: the B-spline design Gram is "
+                f"rank-deficient ({n_dropped} redundant direction(s) dropped); "
+                "rank-reducing to a finite SE on the identified subspace. "
+                "Analytical SEs reflect the reduced rank (NaN if rank 0).",
+                UserWarning,
+                stacklevel=2,
+            )
 
         # ee_treated: per-unit estimating equation vectors (K-vector per unit)
         # For WLS (survey weights), the score is w_i * X_i * u_i to match the

diff_diff/linalg.py

@@ -312,12 +312,33 @@ def _equilibrated_lstsq(X: np.ndarray, y: np.ndarray) -> np.ndarray:
     return coef_scaled / safe_norms
 
 
+@overload
+def _rank_guarded_inv(
+    A: np.ndarray,
+    *,
+    rcond: float = ...,
+    tracker: Optional[list] = ...,
+    return_dropped: Literal[False] = ...,
+) -> Tuple[np.ndarray, int, int]: ...
+
+
+@overload
+def _rank_guarded_inv(
+    A: np.ndarray,
+    *,
+    rcond: float = ...,
+    tracker: Optional[list] = ...,
+    return_dropped: Literal[True],
+) -> Tuple[np.ndarray, int, int, np.ndarray]: ...
+
+
 def _rank_guarded_inv(
     A: np.ndarray,
     *,
     rcond: float = 1e-10,
     tracker: Optional[list] = None,
-) -> Tuple[np.ndarray, int, int]:
+    return_dropped: bool = False,
+) -> Union[Tuple[np.ndarray, int, int], Tuple[np.ndarray, int, int, np.ndarray]]:
     """Rank-guarded (generalized) inverse of a symmetric PSD Gram matrix.
 
     Influence-function standard errors invert a covariate Gram matrix
@@ -384,8 +405,19 @@ def _rank_guarded_inv(
     unchanged, so well-conditioned fits are numerically unaffected.
     """
     k = A.shape[0]
+
+    def _ret(inv, n_dropped, n_keep, dropped):
+        # ``dropped`` is a length-k boolean mask of the truncated (unidentified)
+        # coordinates. Callers that report a PER-COEFFICIENT SE from ``vcov``
+        # diagonals must NaN the dropped coordinates in the FINAL vcov — the
+        # zero-filled inverse would otherwise report ``se=0`` for an unidentified
+        # coefficient. (Linear-combination consumers — e.g. an ATT or dose
+        # prediction — keep the default 3-tuple: the dropped direction correctly
+        # contributes 0 to an identified linear combination.)
+        return (inv, n_dropped, n_keep, dropped) if return_dropped else (inv, n_dropped, n_keep)
+
     if k == 0:
-        return np.zeros((0, 0), dtype=float), 0, 0
+        return _ret(np.zeros((0, 0), dtype=float), 0, 0, np.zeros(0, dtype=bool))
 
     # Symmetric equilibration: scale row/col i by sqrt(A[i, i]) so the
     # eigenvalue threshold is scale-invariant. Zero/negative diagonal -> 1.0
@@ -406,7 +438,7 @@ def _rank_guarded_inv(
 
     # Fast path: full rank -> exact solve (bit-identical to the prior code).
     if max_eig > 0.0 and n_keep == k:
-        return np.linalg.solve(A, np.eye(k)), 0, k
+        return _ret(np.linalg.solve(A, np.eye(k)), 0, k, np.zeros(k, dtype=bool))
 
     # Rank-deficient: record one condition-number sample for the aggregate
     # fallback warning (the helper is the sole owner of this append).
@@ -415,7 +447,7 @@ def _rank_guarded_inv(
             tracker.append(float(np.linalg.cond(A)))
 
     if n_keep == 0:
-        return np.full((k, k), np.nan), k, 0
+        return _ret(np.full((k, k), np.nan), k, 0, np.ones(k, dtype=bool))
 
     # Column-drop generalized inverse: keep the n_keep most-independent columns
     # (pivoted QR on the equilibrated Gram), invert that principal submatrix, and
@@ -431,7 +463,9 @@ def _rank_guarded_inv(
     A_eq_ginv[np.ix_(kept, kept)] = np.linalg.inv(A_eq[np.ix_(kept, kept)])
     A_ginv = A_eq_ginv * inv_scales[:, None] * inv_scales[None, :]
     n_dropped = k - n_keep
-    return A_ginv, n_dropped, n_keep
+    dropped_mask = np.ones(k, dtype=bool)
+    dropped_mask[kept] = False
+    return _ret(A_ginv, n_dropped, n_keep, dropped_mask)
 
 
 def _solve_ols_rust(

diff_diff/spillover.py

@@ -38,7 +38,7 @@
     _haversine_km,
     _validate_callable_metric_result,
 )
-from diff_diff.linalg import solve_ols
+from diff_diff.linalg import _rank_guarded_inv, solve_ols
 from diff_diff.results import SpilloverDiDResults
 from diff_diff.two_stage import _compute_gmm_corrected_meat
 from diff_diff.utils import safe_inference
@@ -3357,27 +3357,38 @@ def fit(
             score_pad_mask=score_pad_mask_arg,
         )
 
-        # Bread sandwich: A_22^{-1} = (X_2' W X_2)^{-1} via `np.linalg.solve`
-        # with dense lstsq fallback + UserWarning (mirrors the bread-fallback
-        # pattern at `two_stage.py:1763-1788`). Wave E.1 adds the W diagonal
-        # under the survey path so the bread aligns with the WLS gamma /
-        # weighted Psi construction in the meat helper.
+        # Bread sandwich: A_22^{-1} = (X_2' W X_2)^{-1} via the shared rank-guarded
+        # generalized inverse `_rank_guarded_inv` (column-drop on a near-singular
+        # Gram + UserWarning; dropped coordinates are NaN'd in the vcov below).
+        # Wave E.1 adds the W diagonal under the survey path so the bread aligns
+        # with the WLS gamma / weighted Psi construction in the meat helper.
         if survey_weights_fit is not None:
             A_22_kept = X_2_kept.T @ (X_2_kept * survey_weights_fit[:, None])
         else:
             A_22_kept = X_2_kept.T @ X_2_kept
-        eye_kept = np.eye(A_22_kept.shape[0])
-        try:
-            bread_kept = np.linalg.solve(A_22_kept, eye_kept)
-        except np.linalg.LinAlgError:
+        # np.linalg.solve only raises on an *exactly* singular Gram; a *near*-
+        # singular A_22 would otherwise flow a garbage inverse (~1e13) into the
+        # SE. `_rank_guarded_inv` truncates redundant directions on the
+        # equilibrated Gram -> finite SE on the identified subspace (NaN at
+        # rank 0), matching the covariate IF rank-guard. A_22_kept is already
+        # column-dropped upstream; this is the within-kept near-singular guard
+        # (its rank-0 all-NaN return composes with the (k,k) re-inflation below).
+        bread_kept, n_dropped, _, dropped = _rank_guarded_inv(A_22_kept, return_dropped=True)
+        if n_dropped:
             warnings.warn(
-                "SpilloverDiD Wave D bread: A_22 = X_2' X_2 is singular; "
-                "falling back to dense lstsq. SE may be unreliable.",
+                "SpilloverDiD Wave D bread: A_22 = X_2' W X_2 is rank-deficient; "
+                "rank-reducing to a finite SE on the identified subspace "
+                f"({n_dropped} redundant direction(s) dropped, NaN if rank 0).",
                 UserWarning,
                 stacklevel=2,
             )
-            bread_kept = np.linalg.lstsq(A_22_kept, eye_kept, rcond=None)[0]
         vcov_kept = bread_kept @ meat_kept @ bread_kept
+        # A within-kept dropped (unidentified) coefficient is zero-filled in
+        # bread_kept, which would report se=0; NaN its row/col so per-coef SE is
+        # NaN, not 0. These ride along the (k, k) re-inflation below.
+        if dropped.any():
+            vcov_kept[dropped, :] = np.nan
+            vcov_kept[:, dropped] = np.nan
 
         # Re-inflate to (k, k) with NaN at rank-deficient column positions
         # so downstream code (which indexes vcov[i, i] for per-coef SE) sees