docs(perf): Phase 2 bootstrap tiled-GEMM numbers + EfficientDiD analytical-stage TODO

CHANGELOG Changed entry with measured before/after; performance-plan.md Phase 2
section (stage split, kernel design, A/B table, rejected sparse variant);
TODO row for the pre-existing EfficientDiD analytical-stage scaling pathology
found while benchmarking; doc-deps note for the shared kernel; refresh a stale
staggered_bootstrap.py line reference in REGISTRY.

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

Author: igerber

CHANGELOG.md

@@ -127,6 +127,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   supersedes it.
 
 ### Changed
+- **CallawaySantAnna multiplier bootstrap rewritten as a fused, column-tiled scatter-GEMM**
+  (EfficientDiD's bootstrap routes through the same kernel). The former loop sliced the
+  `(block × n_units)` weight matrix twice per (g,t) cell per weight block — at a 40-period,
+  10-cohort, 100k-unit panel (390 cells, 999 draws) that was ~95% of bootstrap wall time as
+  pure memory traffic. All perturbation columns (per-cell IFs, overall combined IF,
+  per-event-time combined IFs) are now scattered into a column-tiled influence matrix
+  (byte-capped tiles; each tile replays the bit-identical weight stream via RNG state
+  snapshot/restore) and consumed by one BLAS GEMM per weight block. Measured (5-rep medians,
+  rust backend + Accelerate): bootstrap fits 24.5s → 4.3s (5.7x) at the 4M-row flagship,
+  5.6s → 0.31s (18.3x) no-covariate and 6.5s → 1.2s (5.5x) dr at 2M rows, 4.2s → 1.5s (2.7x)
+  survey-PSU; peak memory flat or lower. Point estimates are bit-identical (the bootstrap
+  never feeds them); bootstrap SEs/CIs/p-values match the previous loop to BLAS
+  reassociation (≤ ~1e-15 relative — far below bootstrap Monte-Carlo error), the same
+  numerics posture as the existing draw-axis chunking. EfficientDiD's unweighted path folds
+  its `1/n` prefactor into the influence column (`W @ (eif/n)` instead of `(W @ eif)/n`) —
+  the same reassociation-level equivalence.
 - **CallawaySantAnna aggregation SE assembly rewritten to O(n_units)** (also inherited by
   `StaggeredTripleDifference` aggregation). The combined influence-function assembly behind
   simple/event-study/group aggregated SEs — previously ~56-85% of analytical fit time at

TODO.md

@@ -47,6 +47,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 
 | Issue | Location | Origin | Effort | Priority |
 |-------|----------|--------|--------|----------|
+| `EfficientDiD` analytical stage has a pathological scaling hot loop: 339s for ONE fit at just 2k units × 20 periods with 5 covariates (15+ CPU-min at 10k units; memory pressure at 100k), while CS handles 100k units in ~2s. Process sampling attributes it to a broadcasted NumPy subtract/multiply inside a Python-level loop (EIF/weights construction). Profile properly and vectorize; until then EfficientDiD is impractical beyond a few thousand units. Found while benchmarking the Phase 2 bootstrap rewrite (its bootstrap stage is 0.018s at 2k units — the analytical stage is the whole problem). | `efficient_did.py` / `efficient_did_weights.py` | CS-scaling | Heavy | Medium |
 | `ImputationDiD` dense `(A0'A0).toarray()` scales `O((U+T+K)^2)` — OOM risk on large panels (only triggers when the sparse solver fails). Needs an alternative dense fallback or richer sparse strategy. | `imputation.py` | #141 | Heavy | Medium |
 | CR2 Bell-McCaffrey DOF uses a naive `O(n²k)` per-coefficient loop over cluster pairs; Pustejovsky-Tipton (2018) Appendix B has a scores-based formulation avoiding the full `n×n` `M`. Switch when a user hits a large-`n` cluster-robust design. | `linalg.py::_compute_cr2_bm` | Phase 1a | Heavy | Low |
 | Rust-backend HC2: the Rust path only supports HC1; HC2 and CR2 Bell-McCaffrey fall through to NumPy. Noticeable for large-`n` fits. | `rust/src/linalg.rs` | Phase 1a | Mid | Low |

docs/doc-deps.yaml

@@ -949,7 +949,7 @@ sources:
       - path: docs/methodology/REGISTRY.md
         section: "Bootstrap, Survey Bootstrap"
         type: methodology
-        note: "Memory-bounded draw-block chunking of multiplier-bootstrap weight generation/consumption (shared by CallawaySantAnna; EfficientDiD/HAD adoption is follow-up). Implementation detail — numerically equivalent (~1 ULP, BLAS reassociation) to the un-chunked path; weight stream bit-identical."
+        note: "Memory-bounded multiplier-bootstrap machinery shared by CallawaySantAnna, EfficientDiD, and HAD: draw-block chunking of weight generation/consumption, plus the fused column-tiled scatter-GEMM perturbation kernel (tiled_if_matmul + ReplayableWeightStream, v3.7 Phase 2). Implementation detail — numerically equivalent (~1 ULP, BLAS reassociation) to the per-column GEMV path; weight stream bit-identical, replayable per column tile."
 
   diff_diff/prep.py:
     drift_risk: low

docs/methodology/REGISTRY.md

@@ -468,7 +468,7 @@ methodology is introduced.
 - `classical`, `hc2`, `hc2_bm`, `conley` — REJECTED at `__init__`. The rejection is **library-architectural, not paper-prescribed**: analytical-sandwich variance families (`classical`, `hc2`, `hc2_bm`) are defined on a single regression's hat matrix, and CS's per-(g,t) doubly-robust / IPW / outcome-regression structure has no equivalent single design matrix to compute hat-matrix leverage or Bell-McCaffrey Satterthwaite DOF on. Spatial-HAC (`conley`) likewise has no defined composition with per-unit IF aggregation today. See ["IF-based variance estimators vs analytical-sandwich estimators"](#if-based-variance-estimators-vs-analytical-sandwich-estimators) above for the structural taxonomy.
 
 *Cluster wiring:*
-Prior to the bare-`cluster=` wiring fix, `CallawaySantAnna(cluster="X")` was a silent no-op — the parameter was stored at `__init__` but never consumed in the fit / aggregator / bootstrap pipeline (users got per-unit IF variance silently, even when they explicitly set `cluster="state"`). The fix synthesizes a minimal `SurveyDesign(psu=X, weight_type="pweight")` when bare `cluster=` is set without an explicit survey design, threading the synthesized PSU through the existing `_compute_stratified_psu_meat` aggregator (`staggered_aggregation.py:735-749`) and PSU-level multiplier bootstrap (`staggered_bootstrap.py:323-347`). Three-branch wiring at `staggered.py:~1500`:
+Prior to the bare-`cluster=` wiring fix, `CallawaySantAnna(cluster="X")` was a silent no-op — the parameter was stored at `__init__` but never consumed in the fit / aggregator / bootstrap pipeline (users got per-unit IF variance silently, even when they explicitly set `cluster="state"`). The fix synthesizes a minimal `SurveyDesign(psu=X, weight_type="pweight")` when bare `cluster=` is set without an explicit survey design, threading the synthesized PSU through the existing `_compute_stratified_psu_meat` aggregator (`staggered_aggregation.py:735-749`) and PSU-level multiplier bootstrap (`staggered_bootstrap.py:334-437`). Three-branch wiring at `staggered.py:~1500`:
   1. Bare `cluster=X` + no `survey_design` → synthesize `SurveyDesign(psu=X)`; refit `_resolve_survey_for_fit` on synthetic; `effective_survey_design = synthetic` so `_validate_unit_constant_survey` runs on it (preventing first-value-wins collapse for movers on panel data).
   2. `survey_design` without PSU + `cluster=X` → call `_inject_cluster_as_psu(resolved_survey, cluster_ids)`.
   3. `survey_design` with PSU + `cluster=X` → PSU wins; `_resolve_effective_cluster` emits `UserWarning` if partitions differ.

docs/performance-plan.md

@@ -183,6 +183,54 @@ precision; analytical SEs to 5 significant digits; bootstrap SEs within ~1% (Mon
 
 ---
 
+## CS multiplier-bootstrap fused tiled-GEMM rewrite (v3.7 Phase 2, 2026-07)
+
+Stage-level instrumentation of `_run_multiplier_bootstrap` at the 40p × 10c × 100k-unit
+flagship (390 cells, 72 event times, 999 draws; 20.9s of a 24.4s fit) split the loop as:
+
+| stage | time | share |
+|---|---:|---:|
+| per-cell fancy-index copies `W[:, idx_t]`, `W[:, idx_c]` | 20.45s | **95%** |
+| per-cell GEMVs | 0.73s | 3% |
+| event-study GEMVs (72 full-width) | 0.27s | 1% |
+| weight generation (rust, 1e8 draws) | 0.04s | 0.2% |
+
+Not FLOPs, not RNG — ~100GB of cache-hostile gather traffic from slicing the
+`(block × n_units)` weight matrix twice per cell (mean cell density 0.33: never-treated
+controls appear in every cell). The fix: scatter every perturbation column (per-cell IFs,
+overall combined IF, per-event-time combined IFs) into a column-tiled influence matrix
+and run one BLAS GEMM per (weight block, column tile) — `bootstrap_chunking.tiled_if_matmul`
++ `ReplayableWeightStream` (each column tile replays the bit-identical weight stream via
+RNG state snapshot/restore; byte-capped tiles keep the kernel's live intermediates bounded
+regardless of n_units or column count). A scipy-sparse cell-column variant was measured and
+rejected (27x slower than dense GEMM at 33% density). EfficientDiD's per-cell full-width
+GEMV bootstrap routes through the same kernel via lazily materialized scaled-EIF columns.
+
+Measured (5-rep medians, M4 Max, rust backend + Accelerate, otherwise-idle machine;
+BEFORE = pristine `origin/main` worktree via PYTHONPATH, AFTER = frozen branch):
+
+| Scenario | Before | After | Speedup | Peak RSS |
+|---|---:|---:|---:|---|
+| 40p × 10c, 5-cov dr, boot-999 (4M rows, 390 cells) | 24.49s | 4.28s | **5.7x** | 2.9→2.9 GB |
+| no-cov boot-999 (2M rows) | 5.64s | 0.31s | **18.3x** | 2.0→1.6 GB |
+| 5-cov dr boot-999 (2M rows) | 6.49s | 1.19s | **5.5x** | 2.3→1.9 GB |
+| survey-PSU 5-cov dr boot-999 (500 PSUs, non-identity map) | 4.22s | 1.54s | **2.7x** | 2.0→1.9 GB |
+| RCS 5-cov dr boot-999 (2M obs) | 42.15s | 40.07s | 1.1x | 3.3→3.2 GB |
+
+ATT deltas exactly 0 on every scenario (the bootstrap never feeds point estimates);
+bootstrap SE relative deltas 3e-16 to 1.2e-15 (BLAS reassociation only). The RCS scenario
+barely moves because its fit time is the analytical repeated-cross-sections path, not the
+bootstrap. EfficientDiD's bootstrap stage went 0.018s → 0.012s but is invisible inside its
+fits: EfficientDiD's **analytical** stage has a pre-existing pathological hot loop
+(339s at just 2k units × 20p with 5 covariates; sampled to a broadcasted subtract/multiply
+in a Python-level loop) — tracked in TODO.md as its own item, out of scope here.
+
+Remaining CS bottlenecks (updated forecast): `_precompute_structures` pandas
+groupby+pivot on analytical no-cov fits, and per-cell DR/IPW nuisance solvers at high
+covariate counts (IRLS `lstsq`, rank-guard QR — Phase 3 candidate).
+
+---
+
 ## Memory scaling at the millions-of-units tail (v3.x, June 2026)
 
 At the scale where the dense working arrays - not compute - are the binding constraint