perf(linalg): opt-in solve_ols normal-equations Cholesky fast path (DIFF_DIFF_SOLVE_OLS_FASTPATH)

solve_ols is the universal OLS entry point; both backends run an equilibrated
thin SVD (gelsd-parity) unconditionally. After #634's stage-0 certification
removed the pivoted-QR cost, the SVD solve itself is the dominant stage on
solver-bound fits (58%/54%/40% of SunAbraham county/firm and CS dr cov40).

Setting DIFF_DIFF_SOLVE_OLS_FASTPATH=1 (positive integer; per-call resolver
mirroring DIFF_DIFF_DEMEAN_CHUNK_COLS) routes certified-well-conditioned
full-rank solves through an equilibrated normal-equations Cholesky:

- numpy twin: reuses the stage-0 rank-certification artifacts (the Gram is
  already built to certify rank), gated by cho_factor + LAPACK dpocon at
  reciprocal condition > 1e-6 (the _IRLS_CHOL_RCOND_GUARD budget: forward
  error ~eps*cond <= 2e-10); self-builds and self-certifies on the
  skip_rank_check route (factorization success alone is not a certificate).
- Rust kernel solve_ols_chol: SEPARATE self-certifying pyfunction (a stale
  extension degrades via the independent-import pattern instead of a
  call-time TypeError): faer Llt with the exact 1-norm rcond whose inverse
  byproduct doubles as the sandwich bread; ONE owned n x k buffer (the
  equilibrated copy, reused in place for the vcov scores - no thin-SVD U
  transient); faer-only heavy ops, so the accelerate/openblas/no-BLAS wheel
  variants behave identically; shared first-appearance cluster indexer with
  the SVD path.
- Fallback contract: any certification decline falls through the UNCHANGED
  legacy chain (Rust SVD, then gelsd) on both dispatch routes, so
  knob-on-decline output equals knob-off exactly; solve_ols gains an
  optional diagnostics_out sink recording which branch ran.

Measured (M4 Max, medians; committed A/B lane with a pristine-base-tree
identity gate): SunAbraham 1.13->0.63 s and 16.77->9.89 s (1.7-1.8x),
CallawaySantAnna dr 40-cov 4.40->3.57 s (1.23x), skip_rank_check micro
1.56x, weighted/survey twin 1.07x; estimate deltas <= 9e-15 ATT / 3e-12
rel SE; rust allocator high-water on the 2.4Mx130 clustered solve
7.27 -> 2.66 GB. Knob-off byte-identical to the base tree on every
benchmark scenario.

Also fixed (found while validating the non-finite-inference contract):
the legacy Rust SVD vcov path silently returned all-Inf vcov on saturated
full-rank designs (n == k) - the (n-k) adjustment divides by zero and the
dispatcher fallback check keyed on NaN only. Both Rust vcov paths now
return the documented all-NaN sentinel (matching the canonical numpy
saturated guard, G >= 2 cluster error keeping precedence), and the
dispatcher rejects any non-finite Rust vcov via a shared helper - on the
skip_rank_check route only Inf reroutes (all-NaN remains the documented
sentinel there; rerouting it to a full-rank-assuming numpy path could
raise on a singular bread). CHANGELOG Fixed entry; sole qualified
exception to the byte-identity claim.

Tests: 57 cargo (LU oracles for beta/hc1/CR1, decline contracts, saturated
NaN, G<2 precedence) + resolver/fast-path/dispatch/kernel/saturated pytest
classes across both backends (property parity on ~20 designs, forced and
natural-guard-trip byte-identity, BETWEEN fixture with self-check,
default-path spy locks, WLS incl. zero weights, hc2/conley spot checks,
20-call clustered bit-determinism, stale-symbol degradation, Inf-injection
reroutes on both routes, NaN-sentinel pass-through).

Docs: CHANGELOG Added+Fixed, REGISTRY Note (thresholds, verbatim-fallback
contract, tol-bounded-not-bit opt-in parity), performance-plan section
with attribution/results/memory tables, doc-deps notes, TODO row swapped
to the default-on-flip follow-up + cross-refs on the parked memory-floor
and QR-reuse rows.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01GPX5Rv8ozQXPdUV23QTfjr

Author: igerber

CHANGELOG.md

@@ -222,6 +222,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   the same dispatch globals the fit consumes (`bootstrap_utils` / `linalg`), with a coherence
   assert between the two (defense in depth).
 
+### Fixed
+- **Rust backend: saturated full-rank designs (`n == k`) no longer leak `Inf`
+  variance-covariance silently.** `solve_ols`'s Rust HC1/CR1 path computed the
+  `n/(n-k)` (clustered `(n-1)/(n-k)`) adjustment with zero residual degrees of
+  freedom, producing all-`Inf` vcov that bypassed the dispatcher's Python-fallback
+  check (which keyed on NaN only) and reached users without a warning. The Rust
+  kernel now returns the documented non-finite-inference contract — an all-NaN
+  vcov, matching the canonical numpy saturated guard — with the `G >= 2` cluster
+  error keeping precedence, and the dispatcher's fallback check now rejects any
+  non-finite vcov (NaN or Inf). Point estimates and residuals are unchanged;
+  only the undefined-inference sentinel changes (`Inf` → `NaN` + warning).
+  Found while building the opt-in Cholesky fast path below, which honors the
+  same contract.
+
 ### Added
 - **`SyntheticControl` conformal extensions: one-sided alternatives + covariates in the
   proxy (Chernozhukov-Wüthrich-Zhu 2021).** `conformal_test` / `conformal_confidence_intervals`
@@ -236,6 +250,26 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   T*-block structure. Locked by a hand-rolled signed-statistic permutation oracle, directional
   rejection/half-line/composition tests, a proxy-weight mechanism check, and a permutation-floor
   invariant; all 31 pre-existing conformal tests pass unchanged.
+- **Opt-in `DIFF_DIFF_SOLVE_OLS_FASTPATH` normal-equations Cholesky fast path for
+  `solve_ols` (both backends).** Set the env var to `1` (a positive integer; boolean
+  spellings like `true`/`on` are silently ignored per the resolver convention, mirroring
+  `DIFF_DIFF_DEMEAN_CHUNK_COLS`) to route certified-well-conditioned full-rank OLS solves
+  through an equilibrated normal-equations Cholesky instead of the default SVD/gelsd:
+  the Python twin reuses the stage-0 rank-certification Gram and gates on LAPACK `dpocon`
+  reciprocal condition > 1e-6; the Rust side is a new self-certifying faer `Llt` kernel
+  (`solve_ols_chol`, faer-only heavy ops — identical behavior across the
+  accelerate/openblas/no-BLAS wheels) whose exact 1-norm rcond inverse doubles as the
+  sandwich bread. Default OFF = byte-identical legacy behavior on both backends (locked
+  by dispatch spies, exact-equality tests, and a pristine-base-tree identity gate; sole
+  exception: the saturated n == k vcov contract fix above); any
+  certification decline falls back verbatim to the SVD path, so a knob-on decline equals
+  knob-off exactly. Within the opt-in path, parity vs the default is tol-bounded (fitted
+  ~1e-8 abs / SE ~1e-6 rel; certified forward-error budget ~eps·cond ≤ 2e-10). Measured
+  (M4 Max, medians): SunAbraham 1.13→0.63 s and 16.8→9.9 s (1.7-1.8x), CallawaySantAnna
+  dr 40-covariate 4.40→3.57 s (1.23x), `skip_rank_check` micro-solve 1.56x; rust-side
+  allocator high-water on a 2.4M×130 clustered solve 7.27→2.66 GB (no thin-SVD U
+  transient). See docs/performance-plan.md § "Opt-in solve_ols normal-equations Cholesky
+  fast path".
 - **Opt-in `df_convention="cluster"` inference-df knob (DiD / TWFE / MultiPeriodDiD +
   `LinearRegression`).** Clustered analytical fits historically compute t-statistics,
   p-values, and CIs at the fitted **residual df** (`n − K_full`), while `fixest`/Stata use

TODO.md

@@ -42,7 +42,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
 | `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 |
+| Evaluate flipping `DIFF_DIFF_SOLVE_OLS_FASTPATH` default-ON after an opt-in soak (the 2026-07 certified normal-equations Cholesky fast path, both backends). A flip needs: golden/parity-suite recapture at the tol-bounded posture (fitted ~1e-8 abs / SE ~1e-6 rel — the default today is byte-pinned in several benchmark conventions), certification-rate telemetry across real workloads (any decline is silent-correct but forfeits the speedup), and the staged default-flip protocol used for `df_convention` (v4-class change). | `diff_diff/linalg.py::_resolve_solve_ols_fastpath`, `rust/src/linalg.rs::solve_ols_chol` | CS-scaling | Mid | Low |
 
 ### Testing / docs
 
@@ -114,8 +114,8 @@ Doable in principle, but no current caller and/or explicitly out of paper scope.
 | CallawaySantAnna **unbalanced-panel R parity — LANDED** via `allow_unbalanced_panel=True` (matches R `did::att_gt(allow_unbalanced_panel=TRUE)` / `DRDID::reg_did_rc`: ATT bit-exact on cells AND dynamic aggregation via fixed unit-cohort-mass `pg` + a per-unit WIF; SE up to the documented CR1 `sqrt(G/(G-1))` factor). The earlier "weighting" framing was a mis-diagnosis — on unbalanced panels the dominant divergence from R is the *estimator* (within-cell differencing vs RC-on-pooled-obs), not only the weighting; both are resolved by the flag. The DEFAULT path keeps within-cell differencing as a documented design choice and now emits a `UserWarning` on unbalanced input (no-silent-failures). **Remaining deferred:** `survey_design=` × `allow_unbalanced_panel=` (per-obs vs per-unit weight resolution — currently fail-closed `NotImplementedError`); and covariate / ipw / dr × the flag R-parity verification (the RC path supports them; the committed golden covers `reg` no-cov). | `staggered.py`, `staggered_aggregation.py` | SE-audit D3 | Low |
 | CallawaySantAnna event-study bucket/weight construction is duplicated between the analytical aggregator (`staggered_aggregation.py::_aggregate_event_study`) and the multiplier bootstrap (`staggered_bootstrap.py`): both group (g,t) by `e = t - g`, apply the finite/NaN/reference masks, and read cohort weights. Both already consume the same source-materialized universal reference cells (so they agree), but the bucket logic is copy-pasted. Extract one shared helper returning per-event-time buckets (finite cells, NaN cells, reference flags, cohort weights, combined-IF inputs) used by both. Pure refactor; gate on byte-identical analytical + bootstrap output. | `staggered_aggregation.py`, `staggered_bootstrap.py` | SE-audit D3 | Low |
 | `StackedDiD` survey re-resolution intra-file dedup (raw-weight extraction ×3, compose-normalize ×3, resolve-on-stacked ×2). The cross-estimator ContinuousDiD/EfficientDiD panel-to-unit collapse consolidation LANDED (#226 shared helpers `ResolvedSurveyDesign.subset_to_units_by_row_idx` / `build_unit_first_row_index`); StackedDiD is deliberately NOT on that path (control units are duplicated across sub-experiments, so it re-resolves at stacked granularity rather than collapsing to one row per unit). The residual is stacked-specific, low value, and touches the numerically-sensitive composed-weight renormalization. Post-filter re-resolution / metadata-recompute unification across the three estimators was assessed and is not warranted — they use genuinely different mechanisms and already delegate to shared `_resolve_survey_for_fit` / `compute_survey_metadata`. | `stacked_did.py` | #226 | Low |
-| `solve_ols` residual memory floor is inherent to SVD-based lstsq after the PR-E marshalling slim (rust-side allocator high-water 15.32 -> 7.81 GB on the 2.4M x 130 clustered solve; remaining = one fused equilibrated input copy + thin-SVD U transient + vcov scores block + LAPACK gelsd internals on the python path). Further reduction requires the tall-skinny-QR algorithm change - same trade-off as the parked QR-reuse row below (library-wide last-digit perturbation + loses the second independent collinearity detector); the two land together if ever reopened. Diagnosis + measurements in docs/performance-plan.md. | `rust/src/linalg.rs::solve_ols`, `diff_diff/linalg.py` | PR-E | Low |
-| SunAbraham-family fits are SOLVER-bound after the v3.6.x demeaning speedups: on the county-class shape (3.1k units x 60 months, ~100 interaction columns) `solve_ols` = ~60% of fit time (faer SVD 1.24s + pivoted-QR rank detection 0.70s of 3.27s), and the share grows further under the Rust demean kernel; the pivoted QR alone is ~20% of the whole fit (measured attribution in docs/performance-plan.md). Candidate fix is QR-reuse (one factorization for rank detection + solve). **Parked per the 2026-07 correctness-first decision**: the change perturbs every coefficient in the library at the last digits (golden/parity/bit-identity recapture) and removes the SVD's independent singular-value truncation - a second, different near-collinearity detector that is deliberate defense-in-depth for the no-silent-failures contract (the FE-span guard episode showed borderline cases are real). Re-open only on practitioner demand for wide event studies; any implementation gates on fixest/R end-to-end parity, not just internal consistency. | `linalg.py::solve_ols`, `linalg.py::_detect_rank_deficiency` | PR-C attribution | Medium |
+| `solve_ols` residual memory floor is inherent to SVD-based lstsq after the PR-E marshalling slim (rust-side allocator high-water 15.32 -> 7.81 GB on the 2.4M x 130 clustered solve; remaining = one fused equilibrated input copy + thin-SVD U transient + vcov scores block + LAPACK gelsd internals on the python path). Further reduction requires the tall-skinny-QR algorithm change - same trade-off as the parked QR-reuse row below (library-wide last-digit perturbation + loses the second independent collinearity detector); the two land together if ever reopened. Diagnosis + measurements in docs/performance-plan.md. (The 2026-07 opt-in `DIFF_DIFF_SOLVE_OLS_FASTPATH` Cholesky path sidesteps the floor when enabled — no thin-SVD U transient, scores reuse the one equilibrated buffer; this row tracks the DEFAULT path only.) | `rust/src/linalg.rs::solve_ols`, `diff_diff/linalg.py` | PR-E | Low |
+| SunAbraham-family fits are SOLVER-bound after the v3.6.x demeaning speedups: on the county-class shape (3.1k units x 60 months, ~100 interaction columns) `solve_ols` = ~60% of fit time (faer SVD 1.24s + pivoted-QR rank detection 0.70s of 3.27s), and the share grows further under the Rust demean kernel; the pivoted QR alone is ~20% of the whole fit (measured attribution in docs/performance-plan.md). Candidate fix is QR-reuse (one factorization for rank detection + solve). **Parked per the 2026-07 correctness-first decision**: the change perturbs every coefficient in the library at the last digits (golden/parity/bit-identity recapture) and removes the SVD's independent singular-value truncation - a second, different near-collinearity detector that is deliberate defense-in-depth for the no-silent-failures contract (the FE-span guard episode showed borderline cases are real). Re-open only on practitioner demand for wide event studies; any implementation gates on fixest/R end-to-end parity, not just internal consistency. (2026-07 updates: the pivoted-QR share quoted above predates #634's stage-0 Gram certification — rank detection is now ~0.03s on the county shape — and the opt-in `DIFF_DIFF_SOLVE_OLS_FASTPATH` Cholesky path is the certification-gated alternative that captures the solver win WITHOUT the library-wide perturbation: the default stays byte-identical and keeps the SVD truncation as the second collinearity detector.) | `linalg.py::solve_ols`, `linalg.py::_detect_rank_deficiency` | PR-C attribution | Medium |
 | dCDH parity-test SE/CI assertions only cover pure-direction scenarios; mixed-direction SE comparison is structurally apples-to-oranges (cell-count vs obs-count weighting). | `test_chaisemartin_dhaultfoeuille_parity.py` | #294 | Low |
 | `HeterogeneousAdoptionDiD` joint cross-horizon covariance / sup-t bands: per-horizon SEs use independent sandwiches (paper-faithful pointwise CIs per Pierce-Schott Fig 2). Follow-ups (low demand): IF-based stacking for joint cross-horizon inference; analytical H×H covariance on the weighted ES path. (A sup-t band on the unweighted ES path shipped via the clustered band — `cluster=` fires the simultaneous band on the unweighted path.) | `had.py::_fit_event_study` | Phase 2b / 4.5 B | Low |
 | `HeterogeneousAdoptionDiD` event-study staggered timing beyond the last cohort: Phase 2b auto-filters to the last cohort (paper App B.2); earlier-cohort effects aren't HAD-identified (redirect to dCDH). Full staggered HAD needs a different identification path (out of paper scope). | `had.py::_validate_had_panel_event_study` | Phase 2b | Low |

benchmarks/speed_review/baselines/solve_ols_fastpath_after.json

@@ -0,0 +1,135 @@
+{
+  "metadata": {
+    "platform": "macOS-26.5.2-arm64-arm-64bit-Mach-O",
+    "arms": [
+      "off",
+      "on"
+    ]
+  },
+  "county_policy": {
+    "off": {
+      "times": [
+        1.127356459001021,
+        1.1081956670022919,
+        1.1512349999975413
+      ],
+      "median": 1.127356459001021,
+      "att": 0.28037854139903057,
+      "se": 0.031244864610215983
+    },
+    "on": {
+      "times": [
+        0.6245858750044135,
+        0.6251335830020253,
+        0.6231642499988084
+      ],
+      "median": 0.6245858750044135,
+      "att": 0.28037854139903007,
+      "se": 0.031244864610213023
+    }
+  },
+  "firm_churn": {
+    "off": {
+      "times": [
+        17.31286900000123,
+        16.22713416699844
+      ],
+      "median": 16.770001583499834,
+      "att": 0.314432718449715,
+      "se": 0.005837281701253922
+    },
+    "on": {
+      "times": [
+        9.93284683299862,
+        9.850229166004283
+      ],
+      "median": 9.891537999501452,
+      "att": 0.3144327184497242,
+      "se": 0.005837281701238088
+    }
+  },
+  "scanner_twfe": {
+    "off": {
+      "times": [
+        0.5774591250010417,
+        0.5748795420004171,
+        0.5841960840043612
+      ],
+      "median": 0.5774591250010417,
+      "att": 0.25116413149499167,
+      "se": 0.002414938742160611
+    },
+    "on": {
+      "times": [
+        0.5765144159959164,
+        0.5650162499950966,
+        0.5698209999973187
+      ],
+      "median": 0.5698209999973187,
+      "att": 0.25116413149499867,
+      "se": 0.0024149387421607276
+    }
+  },
+  "cs_cov40": {
+    "off": {
+      "times": [
+        4.3890298749975045,
+        4.401376791000075,
+        4.414817625001888
+      ],
+      "median": 4.401376791000075,
+      "att": 2.0008128010664223,
+      "se": 0.005610568442727097
+    },
+    "on": {
+      "times": [
+        3.700359333997767,
+        3.5647381669987226,
+        3.5419092910015024
+      ],
+      "median": 3.5647381669987226,
+      "att": 2.0008128010664223,
+      "se": 0.005610568442727097
+    }
+  },
+  "survey_absorb": {
+    "off": {
+      "times": [
+        2.7668485419999342,
+        2.7689070839987835
+      ],
+      "median": 2.767877812999359,
+      "att": 0.24554568144130481,
+      "se": 0.007662698712656388
+    },
+    "on": {
+      "times": [
+        2.5810511670060805,
+        2.595715333998669
+      ],
+      "median": 2.5883832505023747,
+      "att": 0.24554568144130481,
+      "se": 0.00766269871265628
+    }
+  },
+  "skiprank_micro": {
+    "off": {
+      "times": [
+        6.819884582997474,
+        6.934963792002236
+      ],
+      "median": 6.877424187499855,
+      "att": 2.6266509653021792,
+      "se": 0.0006455844619141583
+    },
+    "on": {
+      "times": [
+        4.4150471670000115,
+        4.423518374998821
+      ],
+      "median": 4.419282770999416,
+      "att": 2.626650965303164,
+      "se": 0.0006455844619141326
+    }
+  }
+}

benchmarks/speed_review/baselines/solve_ols_fastpath_before.json

@@ -0,0 +1,77 @@
+{
+  "metadata": {
+    "platform": "macOS-26.5.2-arm64-arm-64bit-Mach-O",
+    "arms": [
+      "off"
+    ]
+  },
+  "county_policy": {
+    "off": {
+      "times": [
+        1.0828709579946008,
+        1.1249838750009076,
+        1.0858678750009858
+      ],
+      "median": 1.0858678750009858,
+      "att": 0.28037854139903057,
+      "se": 0.031244864610215983
+    }
+  },
+  "firm_churn": {
+    "off": {
+      "times": [
+        17.221773457997188,
+        16.04917508300423
+      ],
+      "median": 16.63547427050071,
+      "att": 0.314432718449715,
+      "se": 0.005837281701253922
+    }
+  },
+  "scanner_twfe": {
+    "off": {
+      "times": [
+        0.575249708002957,
+        0.5752674159957678,
+        0.5818307080044178
+      ],
+      "median": 0.5752674159957678,
+      "att": 0.25116413149499167,
+      "se": 0.002414938742160611
+    }
+  },
+  "cs_cov40": {
+    "off": {
+      "times": [
+        4.368334792001406,
+        4.373289832998125,
+        4.364916250000533
+      ],
+      "median": 4.368334792001406,
+      "att": 2.0008128010664223,
+      "se": 0.005610568442727097
+    }
+  },
+  "survey_absorb": {
+    "off": {
+      "times": [
+        2.8271069999973406,
+        2.837882542000443
+      ],
+      "median": 2.832494770998892,
+      "att": 0.24554568144130481,
+      "se": 0.007662698712656388
+    }
+  },
+  "skiprank_micro": {
+    "off": {
+      "times": [
+        7.0087420839990955,
+        7.03886975000205
+      ],
+      "median": 7.023805917000573,
+      "att": 2.6266509653021792,
+      "se": 0.0006455844619141583
+    }
+  }
+}