From 5c895253cb33eae965b83b3c17739ab97589d7a1 Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 11:09:49 +0200
Subject: [PATCH 1/6] exact-CL2 default within the CL2 path (A18) + hardening
 test

Implemented by a codex subagent (gpt-5.6-terra, high), independently
re-verified: exact Bell-McCaffrey block B_g = I - 2 H_gg + (Ht^2)_gg with
(1-cap)^2 eigenvalue floor; cl2_adjustment = auto/exact/shortcut on pffr()
and coef.pffr() (auto: exact for G <= 100 and a dense-cost bound, visible
fallback to shortcut); diagnostics n_adjusted/min_block_eig/max_block_kappa
in fit metadata; docs + NEWS.

A18 hardening gate: 16 cells x {gaussian, poisson, binomial, Gamma} x
k {4,7} x unbalanced/influential-cluster designs - all SEs finite/positive,
3 influential Poisson/Gamma blocks hit the floor and match the uncapped
analogue to 8.8e-6, median SE ratio to shortcut 1.00-1.02
(hardening-exactcl2-results.csv). test-pffr-exactcl2.R: 19 PASS (re-run
independently); full test-pffr.R + sandwich/satterthwaite files passed.

Open: exact Bell-McCaffrey Satterthwaite df remains separate (documented).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md    |  84 ++++++++++++
 NEWS.md                             |   7 +
 R/pffr-core.R                       | 196 +++++++++++++++++++++++-----
 R/pffr-methods.R                    |  15 ++-
 R/pffr.R                            |  16 ++-
 hardening-exactcl2-results.csv      |  17 +++
 hardening-exactcl2.R                | 141 ++++++++++++++++++++
 man/apply_sandwich_correction.Rd    |   6 +-
 man/coef.pffr.Rd                    |   7 +
 man/gam_sandwich_cluster_cl2.Rd     |  12 +-
 man/pffr.Rd                         |  13 +-
 man/pffr_build_metadata.Rd          |   3 +-
 man/pffr_compute_sandwich.Rd        |   6 +-
 man/pffr_vcov.Rd                    |   6 +-
 man/resolve_cl2_adjustment.Rd       |  35 +++++
 tests/testthat/test-pffr-exactcl2.R |  60 +++++++++
 tests/testthat/test-pffr.R          |  14 +-
 17 files changed, 591 insertions(+), 47 deletions(-)
 create mode 100644 EXACTCL2-IMPLEMENTATION-NOTES.md
 create mode 100644 hardening-exactcl2-results.csv
 create mode 100644 hardening-exactcl2.R
 create mode 100644 man/resolve_cl2_adjustment.Rd
 create mode 100644 tests/testthat/test-pffr-exactcl2.R

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
new file mode 100644
index 00000000..940bdca4
--- /dev/null
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -0,0 +1,84 @@
+# Exact-CL2 implementation notes
+
+Date: 2026-07-21
+
+## Change made
+
+`gam_sandwich_cluster_cl2()` now implements the exact Bell--McCaffrey
+leverage block
+
+```
+B_g = I - 2 H_gg + (H_t^2)_gg
+```
+
+where `(H_t^2)_gg` is evaluated as
+`(Xw_g Vp) (Xw' Xw) (Xw_g Vp)'`. Its eigendecomposition floors eigenvalues at
+`(1 - leverage_cap)^2`; the historical shortcut continues to cap the
+eigenvalues of `H_gg` and uses `(I - H_gg)^(-1/2)`.
+
+`pffr()` and `coef.pffr()` have a new `cl2_adjustment` argument:
+`"auto"` (default), `"exact"`, or `"shortcut"`. Fit metadata records the
+resolved adjustment plus `n_adjusted`, `min_block_eig`, and
+`max_block_kappa` for exact CL2. Legacy shortcut diagnostics remain available.
+
+## Default rule
+
+Within the existing CL2 path only, `"auto"` selects exact CL2 when:
+
+* `G <= 100`; and
+* `G * max(D_g) * p^2 <= 5e8`.
+
+The first bound confines exact CL2 to the small-to-moderate-cluster regime
+where its finite-sample correction is relevant. The second is a conservative
+dense-matrix proxy for the per-block multiplication
+`(Xw_g Vp) (Xw'Xw)`. Outside either bound the package visibly records and uses
+the shortcut. The package-level `sandwich = "auto"` policy is unchanged.
+
+## Tests run
+
+All R invocations used `OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1` and
+`nice -n 15`.
+
+* `Rscript hardening-exactcl2.R` — passed; wrote
+  `hardening-exactcl2-results.csv`.
+* Targeted test files run with `NOT_CRAN=true` — passed:
+  `test-pffr-exactcl2.R`, `test-pffr-s2-autopolicy.R`,
+  `test-pffr-sandwich-storage.R`, `test-pffr-sandwich-refit.R`,
+  `test-pffr-satterthwaite.R`, and `test-pffr.R`.
+* `devtools::load_all()` — passed.
+* `devtools::document()` — completed and regenerated the affected Rd files.
+  It printed pre-existing roxygen warnings about unrelated tags/method exports
+  in `pffr-core.R`, `poridge.R`, `pffr-methods.R`, and plot methods.
+
+## Hardening output summary
+
+The 16-cell table covers families Gaussian, Poisson, binomial, and Gamma;
+small/moderate response bases `k = 4, 7`; and unbalanced and
+influential-cluster designs. Each has five unbalanced clusters with curve
+counts `1, 1, 2, 3, 5`; the influential design gives one curve an extreme
+covariate value while retaining an ordinary-range outcome.
+
+Actual results:
+
+* All 16 exact covariances and coefficient SE vectors were finite and strictly
+  positive.
+* The median exact/shortcut SE ratio ranged from 1.00 to 1.02.
+* Three influential Poisson/Gamma blocks hit the shipped exact floor
+  (`n_adjusted = 1`); all other cells reported zero adjusted blocks.
+* The smallest observed pre-floor block eigenvalue was `3.20e-07`.
+* Exact versus the `cl2_exact_nocap` analogue (floor `tol = 1e-8`) differed by
+  at most `8.76e-06`; the nonzero differences occur in the floored cells.
+* The largest reported block condition proxy was `2.98e+06`.
+
+## Open questions
+
+The pointwise Satterthwaite degrees-of-freedom routine remains the established
+working-iid approximation based on the shortcut adjustment; this change
+updates the CL2 covariance block only. Exact Bell--McCaffrey df remains a
+separate methodological task, as already documented in the package.
+
+The requested local commits could not be created in this execution environment:
+Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
+Read-only file system`. The branch `exact-cl2-default` exists and all changes
+remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
+change was touched.
diff --git a/NEWS.md b/NEWS.md
index 72b2fbab..84073aa3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -50,6 +50,13 @@
   report CL2 leverage diagnostics: the returned covariance carries
   `max_leverage` and `n_capped_clusters` attributes, and a warning is emitted
   when one or more clusters hit the leverage cap.
+* CL2 now uses the exact Bell--McCaffrey leverage block by default when the
+  finite-sample correction is relevant and the dense-block calculation is
+  affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <=
+  5e8`). At larger/prohibitive problems it uses the historical shortcut.
+  `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
+  explicitly selects either variant. Exact-CL2 diagnostics include the number
+  of eigenvalue-floored blocks (`n_adjusted`).
 * AR(1) support improvements: `pffr()` now automatically switches to
   `algorithm = "bam"` and `method = "fREML"` when `rho` is supplied, and
   sets `discrete = TRUE` for non-Gaussian families.
diff --git a/R/pffr-core.R b/R/pffr-core.R
index 61c8567a..e225eba0 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -675,7 +675,8 @@ pffr_build_metadata <- function(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   list(
     call = call,
@@ -698,6 +699,7 @@ pffr_build_metadata <- function(
     sandwich = sandwich,
     dof_correction = dof_correction,
     edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment,
     # Covariance storage contract version: format 2 keeps $Vp/$Vc/$Ve
     # model-based ALWAYS; the robust covariance lives in $pffr$Vsandwich.
     cov_format = PFFR_COV_STORAGE_FORMAT,
@@ -1616,6 +1618,10 @@ gam_sandwich_cluster <- function(
 #'   If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).
 #' @param tol Eigenvalue floor for numerical stability.
 #' @param leverage_cap Cap for cluster leverage eigenvalues (< 1).
+#' @param cl2_adjustment CL2 leverage adjustment: `"auto"` (default) selects
+#'   the exact Bell--McCaffrey block where it is relevant and affordable,
+#'   `"exact"` forces that block, and `"shortcut"` uses the historical
+#'   per-cluster shortcut `(I - H_gg)^(-1/2)`.
 #' @param b2 Internal ablation switch (default `TRUE` = current behavior). When
 #'   `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 #'   V_e} (X5).
@@ -1623,8 +1629,10 @@ gam_sandwich_cluster <- function(
 #'   behavior). When `TRUE`, center the leverage-adjusted per-cluster
 #'   contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 #'   (X6).
-#' @returns A p x p covariance matrix with attributes `n_capped_clusters` and
-#'   `max_leverage` for the CL2 leverage diagnostic.
+#' @returns A p x p covariance matrix with leverage diagnostics. Exact CL2
+#'   returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+#'   `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
+#'   `n_capped_clusters` and `max_leverage` attributes.
 #' @keywords internal
 gam_sandwich_cluster_cl2 <- function(
   b,
@@ -1632,6 +1640,7 @@ gam_sandwich_cluster_cl2 <- function(
   freq = FALSE,
   tol = 1e-8,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 ) {
@@ -1673,41 +1682,78 @@ gam_sandwich_cluster_cl2 <- function(
   G <- n_clusters_checked(cluster_id_work)
   groups <- unique(cluster_id_work)
 
-  Vp <- b$Vp
-  B2 <- if (freq) 0 else b$Vp - b$Ve
+  Vp_raw <- b$Vp
+  Vp <- 0.5 * (Vp_raw + t(Vp_raw))
   p <- ncol(Xw)
+  cl2_adjustment <- resolve_cl2_adjustment(
+    cl2_adjustment,
+    G = G,
+    maxDg = max(table(cluster_id_work)),
+    p = p
+  )
+  B2 <- if (freq) {
+    0
+  } else if (cl2_adjustment == "exact") {
+    Vp - b$Ve
+  } else {
+    Vp_raw - b$Ve
+  }
+  Cmat <- if (cl2_adjustment == "exact") crossprod(Xw) else NULL
   meat <- matrix(0, nrow = p, ncol = p)
   Usum <- numeric(p)
   n_capped_clusters <- 0L
   max_leverage <- NA_real_
+  n_adjusted <- 0L
+  min_block_eig <- Inf
+  max_block_kappa <- 0
 
   for (g in groups) {
     idx <- which(cluster_id_work == g)
     Xwg <- Xw[idx, , drop = FALSE]
     zg <- z[idx]
 
-    Hgg <- Xwg %*% Vp %*% t(Xwg)
+    Tg <- Xwg %*% Vp
+    Hgg <- Tg %*% t(Xwg)
     Hgg <- 0.5 * (Hgg + t(Hgg))
 
-    ee_H <- eigen(Hgg, symmetric = TRUE)
-    cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
-    if (is.finite(cluster_max_leverage)) {
-      max_leverage <- if (is.na(max_leverage)) {
-        cluster_max_leverage
-      } else {
-        max(max_leverage, cluster_max_leverage)
+    if (cl2_adjustment == "shortcut") {
+      ee_H <- eigen(Hgg, symmetric = TRUE)
+      cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
+      if (is.finite(cluster_max_leverage)) {
+        max_leverage <- if (is.na(max_leverage)) {
+          cluster_max_leverage
+        } else {
+          max(max_leverage, cluster_max_leverage)
+        }
       }
+      if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
+        n_capped_clusters <- n_capped_clusters + 1L
+        ee_H$values <- pmin(ee_H$values, leverage_cap)
+        Hgg <- ee_H$vectors %*%
+          diag(ee_H$values, nrow = length(ee_H$values)) %*%
+          t(ee_H$vectors)
+      }
+      Ag <- sym_inv_sqrt(diag(length(idx)) - Hgg, tol = tol)
+    } else {
+      # The exact Bell--McCaffrey block is B_g = I - 2 H_gg +
+      # (H_t^2)_gg, with (H_t^2)_gg = (Xw_g Vp) (Xw' Xw)
+      # (Xw_g Vp)'. This is a squared operator, so the gain-matched analogue
+      # of capping H_gg at `leverage_cap` floors B_g at (1 - cap)^2.
+      Bg <- diag(length(idx)) - 2 * Hgg + Tg %*% Cmat %*% t(Tg)
+      ee_B <- eigen(0.5 * (Bg + t(Bg)), symmetric = TRUE)
+      block_min <- min(ee_B$values)
+      min_block_eig <- min(min_block_eig, block_min)
+      max_block_kappa <- max(
+        max_block_kappa,
+        abs(max(ee_B$values)) / max(abs(block_min), 1e-300)
+      )
+      floor_val <- (1 - leverage_cap)^2
+      if (block_min < floor_val) n_adjusted <- n_adjusted + 1L
+      ee_B$values <- pmax(ee_B$values, floor_val)
+      Ag <- ee_B$vectors %*%
+        diag(1 / sqrt(ee_B$values), nrow = length(ee_B$values)) %*%
+        t(ee_B$vectors)
     }
-    if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
-      n_capped_clusters <- n_capped_clusters + 1L
-      ee_H$values <- pmin(ee_H$values, leverage_cap)
-      Hgg <- ee_H$vectors %*%
-        diag(ee_H$values, nrow = length(ee_H$values)) %*%
-        t(ee_H$vectors)
-    }
-
-    Mg <- diag(length(idx)) - Hgg
-    Ag <- sym_inv_sqrt(Mg, tol = tol)
     Ug <- crossprod(Xwg, Ag %*% zg)
     meat <- meat + Ug %*% t(Ug)
     Usum <- Usum + as.vector(Ug)
@@ -1723,9 +1769,14 @@ gam_sandwich_cluster_cl2 <- function(
   V <- hc1 * Vp %*% meat %*% Vp
   if (isTRUE(b2)) V <- V + B2
   V <- 0.5 * (V + t(V))
+  attr(V, "cl2_adjustment") <- cl2_adjustment
   attr(V, "n_capped_clusters") <- n_capped_clusters
   attr(V, "max_leverage") <- max_leverage
-  if (n_capped_clusters > 0) {
+  attr(V, "n_adjusted") <- n_adjusted
+  attr(V, "min_block_eig") <- if (is.finite(min_block_eig)) min_block_eig else
+    NA_real_
+  attr(V, "max_block_kappa") <- max_block_kappa
+  if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) {
     max_leverage_label <- if (is.finite(max_leverage)) {
       sprintf("%.3f", max_leverage)
     } else {
@@ -1750,6 +1801,41 @@ gam_sandwich_cluster_cl2 <- function(
   V
 }
 
+#' Resolve the CL2 leverage adjustment
+#'
+#' Exact CL2 is the default when the finite-sample correction can matter
+#' (`G <= 100`) and the dense block multiplication is affordable. Its
+#' conservative cost proxy `G * maxDg * p^2` covers the dominant
+#' `(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
+#' the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
+#' CL2 uses the historical shortcut because the correction tends to zero and
+#' the exact block has no material finite-sample benefit there.
+#'
+#' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`.
+#' @param G Number of clusters.
+#' @param maxDg Largest cluster size.
+#' @param p Coefficient-space dimension.
+#' @returns `"exact"` or `"shortcut"`.
+#' @keywords internal
+resolve_cl2_adjustment <- function(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+) {
+  cl2_adjustment <- match.arg(cl2_adjustment)
+  if (cl2_adjustment != "auto") return(cl2_adjustment)
+
+  cost <- G * maxDg * p^2
+  if (
+    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e8
+  ) {
+    "exact"
+  } else {
+    "shortcut"
+  }
+}
+
 #' Resolve a pffr fit's covariance matrices to a canonical representation
 #'
 #' The current storage contract (format 2) keeps `Vp`/`Vc`/`Ve` model-based
@@ -1899,6 +1985,8 @@ restore_model_cov <- function(object) {
 #' @param center_scores Internal ablation switch (default `FALSE`); center the
 #'   per-cluster score sums before the meat when `TRUE` (X6). Applies to
 #'   `"cluster"`/`"cl2"` only.
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"`, `"exact"`, or
+#'   `"shortcut"`), used only for `type = "cl2"`.
 #' @returns A covariance matrix (with CL2 leverage attributes for `type =
 #'   "cl2"`).
 #' @keywords internal
@@ -1910,7 +1998,8 @@ pffr_compute_sandwich <- function(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 ) {
   switch(
     type,
@@ -1928,7 +2017,8 @@ pffr_compute_sandwich <- function(
       cluster_id,
       freq = freq,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     ),
     hc = mgcv::vcov.gam(b, sandwich = TRUE, freq = freq),
     none = if (freq) b$Ve else (b$Vc %||% b$Vp),
@@ -1960,6 +2050,8 @@ pffr_compute_sandwich <- function(
 #' @param center_scores Internal ablation switch (default `FALSE` = current
 #'   behavior). When `TRUE`, center the per-cluster score sums before forming
 #'   the cluster/CL2 meat (X6). Same cache semantics as `b2`.
+#' @param cl2_adjustment CL2 leverage adjustment. `NULL` (default) inherits
+#'   the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.
 #' @returns A covariance matrix.
 #' @keywords internal
 pffr_vcov <- function(
@@ -1970,7 +2062,8 @@ pffr_vcov <- function(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 ) {
   # Ablation variants (X5/X6) bypass the fit-time and recompute caches entirely,
   # so they never overwrite or shadow the standard cached matrices.
@@ -1990,6 +2083,13 @@ pffr_vcov <- function(
 
   dof_correction <- dof_correction %||% (object$pffr$dof_correction %||% "none")
   edf_type <- edf_type %||% (object$pffr$edf_type %||% "trace")
+  cl2_adjustment_explicit <- !is.null(cl2_adjustment)
+  cl2_adjustment <- cl2_adjustment %||%
+    (object$pffr$cl2_adjustment %||% "auto")
+  cl2_adjustment <- match.arg(
+    cl2_adjustment,
+    c("auto", "exact", "shortcut")
+  )
 
   # Serve the cached fit-time robust matrix when the request matches it exactly.
   opts_match <- if (requested == "cluster") {
@@ -1999,11 +2099,20 @@ pffr_vcov <- function(
   } else {
     TRUE
   }
+  adjustment_match <- requested != "cl2" ||
+    !cl2_adjustment_explicit ||
+    identical(
+      cl2_adjustment,
+      object$pffr$sandwich_info$cl2_adjustment %||%
+        object$pffr$cl2_adjustment %||%
+        "shortcut"
+    )
   if (
     !ablation &&
       is.null(cluster) &&
       identical(requested, canon$fit_type) &&
       opts_match &&
+      adjustment_match &&
       !is.null(canon$Vsandwich)
   ) {
     return(
@@ -2024,6 +2133,17 @@ pffr_vcov <- function(
   key <- if (!ablation && is.null(cluster)) {
     if (requested == "cluster" && (freq || dof_correction != "none")) {
       paste(requested, freq, dof_correction, edf_type, sep = "|")
+    } else if (requested == "cl2") {
+      if (cl2_adjustment == "auto") {
+        if (freq) paste(requested, "freq", sep = "|") else requested
+      } else {
+        paste(
+          requested,
+          cl2_adjustment,
+          if (freq) "freq" else "bayes",
+          sep = "|"
+        )
+      }
     } else if (freq) {
       paste(requested, "freq", sep = "|")
     } else {
@@ -2054,7 +2174,8 @@ pffr_vcov <- function(
     dof_correction = dof_correction,
     edf_type = edf_type,
     b2 = b2,
-    center_scores = center_scores
+    center_scores = center_scores,
+    cl2_adjustment = cl2_adjustment
   )
   if (!is.null(cache) && !is.null(key)) {
     cache[[key]] <- V
@@ -2956,6 +3077,8 @@ pffr_upgrade_fit <- function(object) {
 #'   warning) for `"cl2"`, which already corrects per-cluster leverage.
 #' @param edf_type Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 #'   `"basis"`).
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"` (default),
+#'   `"exact"`, or `"shortcut"`), used only for `type = "cl2"`.
 #' @returns Model with the robust covariance stored in `$pffr$Vsandwich`.
 #' @keywords internal
 apply_sandwich_correction <- function(
@@ -2963,7 +3086,8 @@ apply_sandwich_correction <- function(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   gam_obj <- if (as.character(algorithm) %in% c("gamm4", "gamm")) m$gam else m
 
@@ -2992,7 +3116,8 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = FALSE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
   Vsw_freq <- pffr_compute_sandwich(
     bread,
@@ -3000,11 +3125,14 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = TRUE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   n_capped <- attr(Vsw, "n_capped_clusters") %||% 0L
   max_lev <- attr(Vsw, "max_leverage") %||% NA_real_
+  n_adjusted <- attr(Vsw, "n_adjusted") %||% 0L
+  resolved_adjustment <- attr(Vsw, "cl2_adjustment") %||% NA_character_
 
   gam_obj$pffr$Vsandwich <- Vsw
   gam_obj$pffr$Vsandwich_freq <- Vsw_freq
@@ -3014,6 +3142,10 @@ apply_sandwich_correction <- function(
     G = if (!is.null(cluster_id)) length(unique(cluster_id)) else NA_integer_,
     n_capped = n_capped,
     max_leverage = max_lev,
+    cl2_adjustment = resolved_adjustment,
+    n_adjusted = n_adjusted,
+    min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_,
+    max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_,
     dof_correction = if (type == "cluster") dof_correction else "none",
     edf_type = edf_type,
     version = as.character(utils::packageVersion("refund")),
@@ -3021,6 +3153,8 @@ apply_sandwich_correction <- function(
   )
   # Keep the legacy CL2 leverage-cap diagnostic slot populated.
   gam_obj$pffr$cl2_n_capped <- if (type == "cl2") n_capped else NULL
+  gam_obj$pffr$cl2_adjustment <- if (type == "cl2") resolved_adjustment else
+    NULL
   # Fresh cache for on-demand recomputation of other sandwich types
   # (fit$pffr$Vsandwich_cache[[type]]).
   gam_obj$pffr$Vsandwich_cache <- new.env(parent = emptyenv())
diff --git a/R/pffr-methods.R b/R/pffr-methods.R
index 05736df3..3324aecc 100755
--- a/R/pffr-methods.R
+++ b/R/pffr-methods.R
@@ -1291,6 +1291,11 @@ compute_pointwise_ci <- function(
 #' @param edf_type Which EDF the \code{"edf"} correction uses
 #'   (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 #'   Defaults to \code{NULL} (inherit from the fit).
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}.
+#'   \code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+#'   applies the documented feasibility rule, while \code{"exact"} and
+#'   \code{"shortcut"} force either CL2 variant. Supplying a value forces a
+#'   covariance recomputation.
 #' @param seWithMean logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?
 #' @param n1 see below
 #' @param n2 see below
@@ -1368,6 +1373,7 @@ coef.pffr <- function(
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -1398,6 +1404,12 @@ coef.pffr <- function(
   if (is.null(edf_type)) edf_type <- model_edf_type
   dof_correction <- match.arg(dof_correction, c("none", "edf"))
   edf_type <- match.arg(edf_type, c("trace", "edf2", "basis"))
+  if (!is.null(cl2_adjustment)) {
+    cl2_adjustment <- match.arg(
+      cl2_adjustment,
+      c("auto", "exact", "shortcut")
+    )
+  }
   # Only warn when the user *explicitly* asked for a dof correction on a
   # non-cluster sandwich; an inherited "edf" (from the fit) stays silent.
   if (dof_explicitly_set && dof_correction != "none" && sandwich != "cluster") {
@@ -1586,7 +1598,8 @@ coef.pffr <- function(
       dof_correction = dof_correction,
       edf_type = edf_type,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     )
 
     # Pointwise critical-value reference (S3). `crit` selects the pointwise
diff --git a/R/pffr.R b/R/pffr.R
index 0f6668da..a85fb987 100644
--- a/R/pffr.R
+++ b/R/pffr.R
@@ -180,7 +180,9 @@
 #'   handles both heteroskedasticity and within-curve autocorrelation — the
 #'   recommended choice for functional data when a fixed estimator is wanted.
 #'   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-#'   style CL2), mainly relevant in smaller samples.
+#'   style CL2), mainly relevant in smaller samples. Within CL2, the exact
+#'   block calculation is used by default where feasible; see
+#'   \code{cl2_adjustment} below.
 #'   \code{"hc"}: observation-level HC sandwich via
 #'   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
 #'   heteroskedasticity but ignores within-curve correlation.
@@ -220,6 +222,13 @@
 #'   the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 #'   \code{"basis"} (the basis dimension). Only relevant when
 #'   \code{dof_correction = "edf"}.
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}:
+#'   \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+#'   dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+#'   10^8}; it otherwise falls back to the historical shortcut. The exact block
+#'   is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+#'   \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+#'   variant. Ignored unless the resolved sandwich is \code{"cl2"}.
 #' @param ... additional arguments that are valid for \code{\link[mgcv]{gam}},
 #'   \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 #'   \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.
@@ -320,6 +329,7 @@ pffr <- function(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 ) {
   call <- match.call()
@@ -330,6 +340,7 @@ pffr <- function(
   sandwich <- match.arg(sandwich)
   dof_correction <- match.arg(dof_correction)
   edf_type <- match.arg(edf_type)
+  cl2_adjustment <- match.arg(cl2_adjustment)
   # "auto" is resolved after fitting (needs G / max D_g); defer the CR1 dof
   # compatibility check to then, so a legitimate auto -> "cluster" keeps its dof.
   if (dof_correction != "none" && !(sandwich %in% c("cluster", "auto"))) {
@@ -503,6 +514,7 @@ pffr <- function(
     prep$algorithm,
     type = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 }
diff --git a/hardening-exactcl2-results.csv b/hardening-exactcl2-results.csv
new file mode 100644
index 00000000..39f55aab
--- /dev/null
+++ b/hardening-exactcl2-results.csv
@@ -0,0 +1,17 @@
+"family","k","design","G","max_cluster_size","n_adjusted","min_block_eig","max_block_kappa","max_abs_exact_nocap","median_se_ratio_to_shortcut","se_finite_positive","se_ratio_sane"
+"gaussian",4,"unbalanced",5,50,0,0.453590234808712,2.20463299969789,0,1.01263947332551,TRUE,TRUE
+"poisson",4,"unbalanced",5,50,0,0.348940004869249,2.82730371704937,0,1.0079082080914,TRUE,TRUE
+"binomial",4,"unbalanced",5,50,0,0.552153568687061,1.81109035005941,0,1.0010220332576,TRUE,TRUE
+"Gamma",4,"unbalanced",5,50,0,0.332486151450086,3.00764406468859,0,1.00903670061469,TRUE,TRUE
+"gaussian",7,"unbalanced",5,50,0,0.374043657671271,2.67348471091803,0,1.0189182605098,TRUE,TRUE
+"poisson",7,"unbalanced",5,50,0,0.301435120875162,3.31746346310504,0,1.00000000257308,TRUE,TRUE
+"binomial",7,"unbalanced",5,50,0,0.478396994803052,2.09031413420915,0,1.00000000004754,TRUE,TRUE
+"Gamma",7,"unbalanced",5,50,0,0.386030967291895,2.590465752049,0,1.00648632821506,TRUE,TRUE
+"gaussian",4,"influential",5,50,0,1.39826045642e-06,709818.235108635,0,1.00791990080075,TRUE,TRUE
+"poisson",4,"influential",5,50,1,3.20272329212834e-07,2982577.58727008,7.06659823389664e-06,1.01755206252338,TRUE,TRUE
+"binomial",4,"influential",5,50,0,1.58358505475409e-06,623867.247044085,0,1.01039988193779,TRUE,TRUE
+"Gamma",4,"influential",5,50,0,2.35791791958206e-06,423665.797661914,0,1.01294453124424,TRUE,TRUE
+"gaussian",7,"influential",5,50,0,1.12859513423658e-06,882176.398770039,0,1.01074965381931,TRUE,TRUE
+"poisson",7,"influential",5,50,1,7.21418514260108e-07,1385176.29210241,8.75780657148509e-07,1.01500324062661,TRUE,TRUE
+"binomial",7,"influential",5,50,0,1.35355708907881e-06,738794.066446026,0,1.00000000006042,TRUE,TRUE
+"Gamma",7,"influential",5,50,1,7.39633434387876e-07,1348528.03666834,4.54114080495374e-07,1.00490798665056,TRUE,TRUE
diff --git a/hardening-exactcl2.R b/hardening-exactcl2.R
new file mode 100644
index 00000000..96a10d1e
--- /dev/null
+++ b/hardening-exactcl2.R
@@ -0,0 +1,141 @@
+# Scoped exact-CL2 hardening gate (DESIGN amendment A18).
+#
+# Run from the package root with:
+# OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 nice -n 15 Rscript hardening-exactcl2.R
+# The script writes hardening-exactcl2-results.csv in the package root.
+
+make_exactcl2_fixture <- function(family_name, k, design) {
+  set.seed(
+    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
+      100 +
+      k +
+      match(design, c("unbalanced", "influential"))
+  )
+  n_curve <- 12L
+  n_grid <- 10L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  eta <- outer(xlin, rep(0.35, n_grid)) +
+    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
+    0.4
+  Y <- switch(
+    family_name,
+    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
+    poisson = matrix(
+      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
+      n_curve
+    ),
+    binomial = matrix(
+      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
+      n_curve
+    ),
+    Gamma = matrix(
+      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
+      n_curve
+    )
+  )
+  # Keep the response in its ordinary range, then make one covariate curve
+  # extreme. This isolates high leverage from a separable/outcome-saturation
+  # failure and exercises the shipped eigenvalue floor (especially Gaussian,
+  # Poisson, and Gamma at k = 9).
+  if (design == "influential") xlin[1] <- 3000
+  # Five deliberately unbalanced independent clusters, with cluster 1 made
+  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
+  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
+  list(
+    data = data.frame(Y = I(Y), xlin = xlin),
+    yind = grid,
+    cluster = cluster,
+    family = switch(
+      family_name,
+      gaussian = stats::gaussian(),
+      poisson = stats::poisson(),
+      binomial = stats::binomial(),
+      Gamma = stats::Gamma(link = "log")
+    ),
+    k = k
+  )
+}
+
+run_exactcl2_hardening <- function() {
+  cases <- expand.grid(
+    family = c("gaussian", "poisson", "binomial", "Gamma"),
+    k = c(4L, 7L),
+    design = c("unbalanced", "influential"),
+    stringsAsFactors = FALSE
+  )
+  rows <- vector("list", nrow(cases))
+
+  for (i in seq_len(nrow(cases))) {
+    case <- cases[i, ]
+    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
+    fit <- suppressWarnings(suppressMessages(refund::pffr(
+      Y ~ xlin,
+      data = fixture$data,
+      yind = fixture$yind,
+      family = fixture$family,
+      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+      sandwich = "none"
+    )))
+
+    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
+    # cl2_exact_nocap comparator without duplicating the package implementation.
+    V_exact <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "exact"
+    ))
+    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
+      refund:::pffr_model_based_gam(fit),
+      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
+      cl2_adjustment = "exact",
+      leverage_cap = 1 - sqrt(1e-8)
+    ))
+    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "shortcut"
+    ))
+
+    se_exact <- sqrt(pmax(diag(V_exact), 0))
+    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
+    positive <- is.finite(se_exact) &
+      is.finite(se_shortcut) &
+      se_exact > 0 &
+      se_shortcut > 0
+    ratio <- se_exact[positive] / se_shortcut[positive]
+    rows[[i]] <- data.frame(
+      family = case$family,
+      k = case$k,
+      design = case$design,
+      G = length(unique(fixture$cluster)),
+      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
+      n_adjusted = attr(V_exact, "n_adjusted"),
+      min_block_eig = attr(V_exact, "min_block_eig"),
+      max_block_kappa = attr(V_exact, "max_block_kappa"),
+      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
+      median_se_ratio_to_shortcut = stats::median(ratio),
+      se_finite_positive = all(positive),
+      se_ratio_sane = length(ratio) > 0 &&
+        min(ratio) > 1e-3 &&
+        max(ratio) < 1e3,
+      stringsAsFactors = FALSE
+    )
+  }
+  do.call(rbind, rows)
+}
+
+if (sys.nframe() == 0L) {
+  devtools::load_all(quiet = TRUE)
+  results <- run_exactcl2_hardening()
+  utils::write.csv(results, "hardening-exactcl2-results.csv", row.names = FALSE)
+  print(results)
+  stopifnot(
+    all(results$se_finite_positive),
+    all(results$se_ratio_sane),
+    all(is.finite(results$n_adjusted)),
+    any(results$n_adjusted > 0)
+  )
+}
diff --git a/man/apply_sandwich_correction.Rd b/man/apply_sandwich_correction.Rd
index 2a497f12..904a494d 100644
--- a/man/apply_sandwich_correction.Rd
+++ b/man/apply_sandwich_correction.Rd
@@ -9,7 +9,8 @@ apply_sandwich_correction(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -26,6 +27,9 @@ warning) for `"cl2"`, which already corrects per-cluster leverage.}
 
 \item{edf_type}{Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 `"basis"`).}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"` (default),
+`"exact"`, or `"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 Model with the robust covariance stored in `$pffr$Vsandwich`.
diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd
index 11eceb3b..a6e8d4de 100644
--- a/man/coef.pffr.Rd
+++ b/man/coef.pffr.Rd
@@ -13,6 +13,7 @@
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -68,6 +69,12 @@ Ignored (with a warning) for \code{sandwich} other than \code{"cluster"}.}
 (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 Defaults to \code{NULL} (inherit from the fit).}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}.
+\code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+applies the documented feasibility rule, while \code{"exact"} and
+\code{"shortcut"} force either CL2 variant. Supplying a value forces a
+covariance recomputation.}
+
 \item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?}
 
 \item{n1}{see below}
diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd
index f5d1ac61..3bd0e45f 100644
--- a/man/gam_sandwich_cluster_cl2.Rd
+++ b/man/gam_sandwich_cluster_cl2.Rd
@@ -10,6 +10,7 @@ gam_sandwich_cluster_cl2(
   freq = FALSE,
   tol = 1e-08,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 )
@@ -27,6 +28,11 @@ If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).}
 
 \item{leverage_cap}{Cap for cluster leverage eigenvalues (< 1).}
 
+\item{cl2_adjustment}{CL2 leverage adjustment: `"auto"` (default) selects
+the exact Bell--McCaffrey block where it is relevant and affordable,
+`"exact"` forces that block, and `"shortcut"` uses the historical
+per-cluster shortcut `(I - H_gg)^(-1/2)`.}
+
 \item{b2}{Internal ablation switch (default `TRUE` = current behavior). When
 `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 V_e} (X5).}
@@ -37,8 +43,10 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 (X6).}
 }
 \value{
-A p x p covariance matrix with attributes `n_capped_clusters` and
-  `max_leverage` for the CL2 leverage diagnostic.
+A p x p covariance matrix with leverage diagnostics. Exact CL2
+  returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+  `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
+  `n_capped_clusters` and `max_leverage` attributes.
 }
 \description{
 Computes a cluster-robust covariance matrix with a Bell-McCaffrey style
diff --git a/man/pffr.Rd b/man/pffr.Rd
index 4fe62cae..769aae24 100644
--- a/man/pffr.Rd
+++ b/man/pffr.Rd
@@ -17,6 +17,7 @@ pffr(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 )
 }
@@ -81,7 +82,9 @@ penalty.}
   handles both heteroskedasticity and within-curve autocorrelation — the
   recommended choice for functional data when a fixed estimator is wanted.
   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-  style CL2), mainly relevant in smaller samples.
+  style CL2), mainly relevant in smaller samples. Within CL2, the exact
+  block calculation is used by default where feasible; see
+  \code{cl2_adjustment} below.
   \code{"hc"}: observation-level HC sandwich via
   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
   heteroskedasticity but ignores within-curve correlation.
@@ -124,6 +127,14 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 \code{"basis"} (the basis dimension). Only relevant when
 \code{dof_correction = "edf"}.}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}:
+\code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+10^8}; it otherwise falls back to the historical shortcut. The exact block
+is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+\eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+variant. Ignored unless the resolved sandwich is \code{"cl2"}.}
+
 \item{...}{additional arguments that are valid for \code{\link[mgcv]{gam}},
 \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.}
diff --git a/man/pffr_build_metadata.Rd b/man/pffr_build_metadata.Rd
index 925ec667..9d183451 100644
--- a/man/pffr_build_metadata.Rd
+++ b/man/pffr_build_metadata.Rd
@@ -24,7 +24,8 @@ pffr_build_metadata(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \description{
diff --git a/man/pffr_compute_sandwich.Rd b/man/pffr_compute_sandwich.Rd
index 5d60df7e..6069f230 100644
--- a/man/pffr_compute_sandwich.Rd
+++ b/man/pffr_compute_sandwich.Rd
@@ -12,7 +12,8 @@ pffr_compute_sandwich(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -34,6 +35,9 @@ only).}
 \item{center_scores}{Internal ablation switch (default `FALSE`); center the
 per-cluster score sums before the meat when `TRUE` (X6). Applies to
 `"cluster"`/`"cl2"` only.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"`, `"exact"`, or
+`"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 A covariance matrix (with CL2 leverage attributes for `type =
diff --git a/man/pffr_vcov.Rd b/man/pffr_vcov.Rd
index 051bc665..697cbe58 100644
--- a/man/pffr_vcov.Rd
+++ b/man/pffr_vcov.Rd
@@ -12,7 +12,8 @@ pffr_vcov(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 )
 }
 \arguments{
@@ -39,6 +40,9 @@ recomputation and are never served from or written to the cache.}
 \item{center_scores}{Internal ablation switch (default `FALSE` = current
 behavior). When `TRUE`, center the per-cluster score sums before forming
 the cluster/CL2 meat (X6). Same cache semantics as `b2`.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment. `NULL` (default) inherits
+the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.}
 }
 \value{
 A covariance matrix.
diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd
new file mode 100644
index 00000000..bdc3c084
--- /dev/null
+++ b/man/resolve_cl2_adjustment.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pffr-core.R
+\name{resolve_cl2_adjustment}
+\alias{resolve_cl2_adjustment}
+\title{Resolve the CL2 leverage adjustment}
+\usage{
+resolve_cl2_adjustment(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+)
+}
+\arguments{
+\item{cl2_adjustment}{One of `"auto"`, `"exact"`, or `"shortcut"`.}
+
+\item{G}{Number of clusters.}
+
+\item{maxDg}{Largest cluster size.}
+
+\item{p}{Coefficient-space dimension.}
+}
+\value{
+`"exact"` or `"shortcut"`.
+}
+\description{
+Exact CL2 is the default when the finite-sample correction can matter
+(`G <= 100`) and the dense block multiplication is affordable. Its
+conservative cost proxy `G * maxDg * p^2` covers the dominant
+`(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
+the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
+CL2 uses the historical shortcut because the correction tends to zero and
+the exact block has no material finite-sample benefit there.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
new file mode 100644
index 00000000..0b2db55f
--- /dev/null
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -0,0 +1,60 @@
+# A18 exact-CL2 hardening gate. The standalone script writes the same table
+# for release records; this fast test keeps its family/design coverage in CI.
+
+source(testthat::test_path("..", "..", "hardening-exactcl2.R"), local = TRUE)
+
+test_that("exact CL2 is robust on high-leverage family and basis fixtures", {
+  results <- run_exactcl2_hardening()
+
+  expect_equal(nrow(results), 16L)
+  expect_true(all(is.finite(results$n_adjusted)))
+  expect_gt(sum(results$n_adjusted), 0)
+  expect_true(all(is.finite(results$min_block_eig)))
+  expect_true(all(is.finite(results$max_block_kappa)))
+  expect_true(all(is.finite(results$max_abs_exact_nocap)))
+  expect_true(all(results$se_finite_positive))
+  expect_true(all(results$se_ratio_sane))
+})
+
+test_that("CL2 adjustment selection can be forced and is recorded", {
+  fixture <- make_exactcl2_fixture("gaussian", 4L, "influential")
+  fit_auto <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2"
+  )))
+  fit_exact <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "exact"
+  )))
+  fit_shortcut <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "shortcut"
+  )))
+
+  expect_identical(fit_exact$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_identical(fit_shortcut$pffr$sandwich_info$cl2_adjustment, "shortcut")
+  expect_identical(fit_auto$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_equal(fit_auto$pffr$Vsandwich, fit_exact$pffr$Vsandwich)
+  expect_true(is.numeric(fit_exact$pffr$sandwich_info$n_adjusted))
+  expect_true(is.numeric(fit_shortcut$pffr$sandwich_info$n_adjusted))
+})
+
+test_that("automatic exact CL2 uses the documented relevance and cost rule", {
+  resolve <- refund:::resolve_cl2_adjustment
+  expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact")
+  expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut")
+  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "shortcut")
+  expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
+  expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
+})
diff --git a/tests/testthat/test-pffr.R b/tests/testthat/test-pffr.R
index c9dca72c..37ad8d3a 100644
--- a/tests/testthat/test-pffr.R
+++ b/tests/testthat/test-pffr.R
@@ -1855,17 +1855,18 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     V_benign <- gam_sandwich_cluster_cl2(b, cluster_id, freq = FALSE),
     NA
   )
-  expect_identical(attr(V_benign, "n_capped_clusters"), 0L)
-  expect_true(is.finite(attr(V_benign, "max_leverage")))
-  expect_gt(attr(V_benign, "max_leverage"), 0)
-  expect_lt(attr(V_benign, "max_leverage"), 0.999)
+  expect_identical(attr(V_benign, "cl2_adjustment"), "exact")
+  expect_identical(attr(V_benign, "n_adjusted"), 0L)
+  expect_true(is.finite(attr(V_benign, "min_block_eig")))
+  expect_gt(attr(V_benign, "min_block_eig"), 0)
 
   expect_warning(
     V_capped <- gam_sandwich_cluster_cl2(
       b,
       cluster_id,
       freq = FALSE,
-      leverage_cap = 0.001
+      leverage_cap = 0.001,
+      cl2_adjustment = "shortcut"
     ),
     "CL2 leverage adjustment hit the leverage cap"
   )
@@ -1876,7 +1877,8 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     b,
     cluster_id,
     freq = FALSE,
-    leverage_cap = 0.001
+    leverage_cap = 0.001,
+    cl2_adjustment = "shortcut"
   ))
   expect_equal(
     as.matrix(V_capped),

From 33d8fad840fcc481347efc3b809faa126305b65c Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 11:55:34 +0200
Subject: [PATCH 2/6] auto-rule cost cap 5e8 -> 5e9, calibrated by measurement

Both adjustments pay the same per-cluster eigendecomposition, so the exact
block's marginal cost is only the extra multiplications: timed at ~0.03 s
extra at proxy 5e8 and ~0.7 s total at 5e9 (single-thread BLAS). The old
cap could deny a sub-second correctness upgrade in the saturated small-G
regime where Study EX measured the largest gain (+1.3pp at G=20). G <= 100
relevance bound unchanged (EX-validated range; gain < 0.15pp beyond).
Boundary tests updated; 20 PASS.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md    | 10 ++++++++++
 NEWS.md                             |  2 +-
 R/pffr-core.R                       | 19 ++++++++++++-------
 man/resolve_cl2_adjustment.Rd       | 17 +++++++++++------
 tests/testthat/test-pffr-exactcl2.R |  5 ++++-
 5 files changed, 38 insertions(+), 15 deletions(-)

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
index 940bdca4..0ae36851 100644
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -82,3 +82,13 @@ Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
 Read-only file system`. The branch `exact-cl2-default` exists and all changes
 remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
 change was touched.
+
+## Addendum 2026-07-21 (PI-side review): cost cap raised 5e8 -> 5e9
+
+Timing at the shipped cap showed the exact block's marginal cost over the
+shortcut is ~0.03 s (both paths pay the same per-cluster eigendecomposition;
+the marginal cost is only the extra multiplications), and ~0.7 s total at
+10x the cap. 5e8 could therefore deny a sub-second correctness upgrade in
+the saturated small-G regime where Study EX measured the largest gain
+(+1.3pp at G=20). Cap raised to 5e9; G <= 100 unchanged (EX-validated
+range; gain < 0.15pp beyond).
diff --git a/NEWS.md b/NEWS.md
index 84073aa3..a0e4f25f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -53,7 +53,7 @@
 * CL2 now uses the exact Bell--McCaffrey leverage block by default when the
   finite-sample correction is relevant and the dense-block calculation is
   affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <=
-  5e8`). At larger/prohibitive problems it uses the historical shortcut.
+  5e9`, a measured sub-second marginal cost). At larger/prohibitive problems it uses the historical shortcut.
   `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
   explicitly selects either variant. Exact-CL2 diagnostics include the number
   of eigenvalue-floored blocks (`n_adjusted`).
diff --git a/R/pffr-core.R b/R/pffr-core.R
index e225eba0..d3712dc5 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -1804,12 +1804,17 @@ gam_sandwich_cluster_cl2 <- function(
 #' Resolve the CL2 leverage adjustment
 #'
 #' Exact CL2 is the default when the finite-sample correction can matter
-#' (`G <= 100`) and the dense block multiplication is affordable. Its
-#' conservative cost proxy `G * maxDg * p^2` covers the dominant
-#' `(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
-#' the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
-#' CL2 uses the historical shortcut because the correction tends to zero and
-#' the exact block has no material finite-sample benefit there.
+#' (`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+#' decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+#' block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+#' covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+#' calibrated by measurement, not guessed: both adjustments pay the same
+#' per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+#' these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+#' total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+#' upgrade at well under a second and in particular never denies it in the
+#' saturated small-G regime where its benefit is largest. Large-G CL2 uses
+#' the historical shortcut because the correction tends to zero there.
 #'
 #' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`.
 #' @param G Number of clusters.
@@ -1828,7 +1833,7 @@ resolve_cl2_adjustment <- function(
 
   cost <- G * maxDg * p^2
   if (
-    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e8
+    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e9
   ) {
     "exact"
   } else {
diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd
index bdc3c084..009b4c2f 100644
--- a/man/resolve_cl2_adjustment.Rd
+++ b/man/resolve_cl2_adjustment.Rd
@@ -25,11 +25,16 @@ resolve_cl2_adjustment(
 }
 \description{
 Exact CL2 is the default when the finite-sample correction can matter
-(`G <= 100`) and the dense block multiplication is affordable. Its
-conservative cost proxy `G * maxDg * p^2` covers the dominant
-`(Xw_g Vp) (Xw'Xw)` multiplication; 5e8 scalar operations keeps this below
-the cost of a typical pffr fit on the supported auto-CL2 grids. Large-G
-CL2 uses the historical shortcut because the correction tends to zero and
-the exact block has no material finite-sample benefit there.
+(`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+calibrated by measurement, not guessed: both adjustments pay the same
+per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+upgrade at well under a second and in particular never denies it in the
+saturated small-G regime where its benefit is largest. Large-G CL2 uses
+the historical shortcut because the correction tends to zero there.
 }
 \keyword{internal}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
index 0b2db55f..d4faf7af 100644
--- a/tests/testthat/test-pffr-exactcl2.R
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -54,7 +54,10 @@ test_that("automatic exact CL2 uses the documented relevance and cost rule", {
   resolve <- refund:::resolve_cl2_adjustment
   expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact")
   expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut")
-  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "shortcut")
+  # 9e8 ops: below the 5e9 cost cap (raised from 5e8, see notes 2026-07-21)
+  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "exact")
+  # 6.4e10 ops: above the cap even at eligible G
+  expect_identical(resolve("auto", G = 100, maxDg = 1000, p = 800), "shortcut")
   expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
   expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
 })

From 9a4e3c8b3f250b2dd1258dd631667a98d9a62fef Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:21:05 +0200
Subject: [PATCH 3/6] review fixes: stale 5e8 docs -> 5e9; wire cl2_adjustment
 into metadata call

Independent review (+ Copilot comments 1-3, all valid): the cap raise in
33d8fad8 missed the pffr.R roxygen (and its generated Rd) and the notes
file's Default-rule body. All now 5e9. Also pass cl2_adjustment through to
pffr_build_metadata() - previously the parameter existed but the sole
caller never supplied it (the field was authoritative only after
apply_sandwich_correction overwrote it).

Copilot comment 4 (use max(abs(eigenvalues)) for kappa) NOT actioned:
B_g is a principal submatrix of the PSD matrix (I-H)^2, hence PSD
(verified numerically, min block eigenvalue +0.23), so abs(max(v)) ==
max(abs(v)) and the diagnostic is already the true condition number.

Review also verified the exact-block identity against a dense (I-H)^2
reference to 4.4e-16 incl. cross-cluster terms. Tests: 20/20.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md | 2 +-
 R/pffr.R                         | 6 ++++--
 man/pffr.Rd                      | 3 ++-
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
index 0ae36851..84710e00 100644
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ b/EXACTCL2-IMPLEMENTATION-NOTES.md
@@ -26,7 +26,7 @@ resolved adjustment plus `n_adjusted`, `min_block_eig`, and
 Within the existing CL2 path only, `"auto"` selects exact CL2 when:
 
 * `G <= 100`; and
-* `G * max(D_g) * p^2 <= 5e8`.
+* `G * max(D_g) * p^2 <= 5e9` (raised from 5e8, see Addendum below).
 
 The first bound confines exact CL2 to the small-to-moderate-cluster regime
 where its finite-sample correction is relevant. The second is a conservative
diff --git a/R/pffr.R b/R/pffr.R
index a85fb987..8f71ce73 100644
--- a/R/pffr.R
+++ b/R/pffr.R
@@ -225,7 +225,8 @@
 #' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}:
 #'   \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
 #'   dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
-#'   10^8}; it otherwise falls back to the historical shortcut. The exact block
+#'   10^9} (a measured sub-second marginal cost); it otherwise falls back to
+#'   the historical shortcut. The exact block
 #'   is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
 #'   \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
 #'   variant. Ignored unless the resolved sandwich is \code{"cl2"}.
@@ -496,7 +497,8 @@ pffr <- function(
     ydata = prep$ydata,
     sandwich = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   m <- pffr_attach_metadata(m, prep$algorithm, ret)
diff --git a/man/pffr.Rd b/man/pffr.Rd
index 769aae24..1521b53f 100644
--- a/man/pffr.Rd
+++ b/man/pffr.Rd
@@ -130,7 +130,8 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 \item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}:
 \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
 dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
-10^8}; it otherwise falls back to the historical shortcut. The exact block
+10^9} (a measured sub-second marginal cost); it otherwise falls back to
+the historical shortcut. The exact block
 is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
 \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
 variant. Ignored unless the resolved sandwich is \code{"cl2"}.}

From 39b0a356c97cb974337a5baedb9a2a65463969ca Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:30:43 +0200
Subject: [PATCH 4/6] move implementation notes + hardening results out of the
 package tree

Both are PR documentation, not package files (non-standard top-level files
would also trip R CMD check); their full content is attached to PR #121 as
a comment. The hardening *script* stays on the branch so the results
remain regenerable.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 EXACTCL2-IMPLEMENTATION-NOTES.md | 94 --------------------------------
 hardening-exactcl2-results.csv   | 17 ------
 2 files changed, 111 deletions(-)
 delete mode 100644 EXACTCL2-IMPLEMENTATION-NOTES.md
 delete mode 100644 hardening-exactcl2-results.csv

diff --git a/EXACTCL2-IMPLEMENTATION-NOTES.md b/EXACTCL2-IMPLEMENTATION-NOTES.md
deleted file mode 100644
index 84710e00..00000000
--- a/EXACTCL2-IMPLEMENTATION-NOTES.md
+++ /dev/null
@@ -1,94 +0,0 @@
-# Exact-CL2 implementation notes
-
-Date: 2026-07-21
-
-## Change made
-
-`gam_sandwich_cluster_cl2()` now implements the exact Bell--McCaffrey
-leverage block
-
-```
-B_g = I - 2 H_gg + (H_t^2)_gg
-```
-
-where `(H_t^2)_gg` is evaluated as
-`(Xw_g Vp) (Xw' Xw) (Xw_g Vp)'`. Its eigendecomposition floors eigenvalues at
-`(1 - leverage_cap)^2`; the historical shortcut continues to cap the
-eigenvalues of `H_gg` and uses `(I - H_gg)^(-1/2)`.
-
-`pffr()` and `coef.pffr()` have a new `cl2_adjustment` argument:
-`"auto"` (default), `"exact"`, or `"shortcut"`. Fit metadata records the
-resolved adjustment plus `n_adjusted`, `min_block_eig`, and
-`max_block_kappa` for exact CL2. Legacy shortcut diagnostics remain available.
-
-## Default rule
-
-Within the existing CL2 path only, `"auto"` selects exact CL2 when:
-
-* `G <= 100`; and
-* `G * max(D_g) * p^2 <= 5e9` (raised from 5e8, see Addendum below).
-
-The first bound confines exact CL2 to the small-to-moderate-cluster regime
-where its finite-sample correction is relevant. The second is a conservative
-dense-matrix proxy for the per-block multiplication
-`(Xw_g Vp) (Xw'Xw)`. Outside either bound the package visibly records and uses
-the shortcut. The package-level `sandwich = "auto"` policy is unchanged.
-
-## Tests run
-
-All R invocations used `OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1` and
-`nice -n 15`.
-
-* `Rscript hardening-exactcl2.R` — passed; wrote
-  `hardening-exactcl2-results.csv`.
-* Targeted test files run with `NOT_CRAN=true` — passed:
-  `test-pffr-exactcl2.R`, `test-pffr-s2-autopolicy.R`,
-  `test-pffr-sandwich-storage.R`, `test-pffr-sandwich-refit.R`,
-  `test-pffr-satterthwaite.R`, and `test-pffr.R`.
-* `devtools::load_all()` — passed.
-* `devtools::document()` — completed and regenerated the affected Rd files.
-  It printed pre-existing roxygen warnings about unrelated tags/method exports
-  in `pffr-core.R`, `poridge.R`, `pffr-methods.R`, and plot methods.
-
-## Hardening output summary
-
-The 16-cell table covers families Gaussian, Poisson, binomial, and Gamma;
-small/moderate response bases `k = 4, 7`; and unbalanced and
-influential-cluster designs. Each has five unbalanced clusters with curve
-counts `1, 1, 2, 3, 5`; the influential design gives one curve an extreme
-covariate value while retaining an ordinary-range outcome.
-
-Actual results:
-
-* All 16 exact covariances and coefficient SE vectors were finite and strictly
-  positive.
-* The median exact/shortcut SE ratio ranged from 1.00 to 1.02.
-* Three influential Poisson/Gamma blocks hit the shipped exact floor
-  (`n_adjusted = 1`); all other cells reported zero adjusted blocks.
-* The smallest observed pre-floor block eigenvalue was `3.20e-07`.
-* Exact versus the `cl2_exact_nocap` analogue (floor `tol = 1e-8`) differed by
-  at most `8.76e-06`; the nonzero differences occur in the floored cells.
-* The largest reported block condition proxy was `2.98e+06`.
-
-## Open questions
-
-The pointwise Satterthwaite degrees-of-freedom routine remains the established
-working-iid approximation based on the shortcut adjustment; this change
-updates the CL2 covariance block only. Exact Bell--McCaffrey df remains a
-separate methodological task, as already documented in the package.
-
-The requested local commits could not be created in this execution environment:
-Git returned `Unable to create .../.git/worktrees/refund-wt-s1/index.lock:
-Read-only file system`. The branch `exact-cl2-default` exists and all changes
-remain unstaged in this worktree; no unrelated `tests/testthat/Rplots.pdf`
-change was touched.
-
-## Addendum 2026-07-21 (PI-side review): cost cap raised 5e8 -> 5e9
-
-Timing at the shipped cap showed the exact block's marginal cost over the
-shortcut is ~0.03 s (both paths pay the same per-cluster eigendecomposition;
-the marginal cost is only the extra multiplications), and ~0.7 s total at
-10x the cap. 5e8 could therefore deny a sub-second correctness upgrade in
-the saturated small-G regime where Study EX measured the largest gain
-(+1.3pp at G=20). Cap raised to 5e9; G <= 100 unchanged (EX-validated
-range; gain < 0.15pp beyond).
diff --git a/hardening-exactcl2-results.csv b/hardening-exactcl2-results.csv
deleted file mode 100644
index 39f55aab..00000000
--- a/hardening-exactcl2-results.csv
+++ /dev/null
@@ -1,17 +0,0 @@
-"family","k","design","G","max_cluster_size","n_adjusted","min_block_eig","max_block_kappa","max_abs_exact_nocap","median_se_ratio_to_shortcut","se_finite_positive","se_ratio_sane"
-"gaussian",4,"unbalanced",5,50,0,0.453590234808712,2.20463299969789,0,1.01263947332551,TRUE,TRUE
-"poisson",4,"unbalanced",5,50,0,0.348940004869249,2.82730371704937,0,1.0079082080914,TRUE,TRUE
-"binomial",4,"unbalanced",5,50,0,0.552153568687061,1.81109035005941,0,1.0010220332576,TRUE,TRUE
-"Gamma",4,"unbalanced",5,50,0,0.332486151450086,3.00764406468859,0,1.00903670061469,TRUE,TRUE
-"gaussian",7,"unbalanced",5,50,0,0.374043657671271,2.67348471091803,0,1.0189182605098,TRUE,TRUE
-"poisson",7,"unbalanced",5,50,0,0.301435120875162,3.31746346310504,0,1.00000000257308,TRUE,TRUE
-"binomial",7,"unbalanced",5,50,0,0.478396994803052,2.09031413420915,0,1.00000000004754,TRUE,TRUE
-"Gamma",7,"unbalanced",5,50,0,0.386030967291895,2.590465752049,0,1.00648632821506,TRUE,TRUE
-"gaussian",4,"influential",5,50,0,1.39826045642e-06,709818.235108635,0,1.00791990080075,TRUE,TRUE
-"poisson",4,"influential",5,50,1,3.20272329212834e-07,2982577.58727008,7.06659823389664e-06,1.01755206252338,TRUE,TRUE
-"binomial",4,"influential",5,50,0,1.58358505475409e-06,623867.247044085,0,1.01039988193779,TRUE,TRUE
-"Gamma",4,"influential",5,50,0,2.35791791958206e-06,423665.797661914,0,1.01294453124424,TRUE,TRUE
-"gaussian",7,"influential",5,50,0,1.12859513423658e-06,882176.398770039,0,1.01074965381931,TRUE,TRUE
-"poisson",7,"influential",5,50,1,7.21418514260108e-07,1385176.29210241,8.75780657148509e-07,1.01500324062661,TRUE,TRUE
-"binomial",7,"influential",5,50,0,1.35355708907881e-06,738794.066446026,0,1.00000000006042,TRUE,TRUE
-"Gamma",7,"influential",5,50,1,7.39633434387876e-07,1348528.03666834,4.54114080495374e-07,1.00490798665056,TRUE,TRUE

From 276ba7179e5db2702cea3de9110b3d99ec0d9484 Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 21 Jul 2026 12:52:40 +0200
Subject: [PATCH 5/6] remove the hardening script from the package tree as well

Non-standard top-level file; full script preserved as PR #121
documentation alongside the notes and results it generates.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95
---
 hardening-exactcl2.R | 141 -------------------------------------------
 1 file changed, 141 deletions(-)
 delete mode 100644 hardening-exactcl2.R

diff --git a/hardening-exactcl2.R b/hardening-exactcl2.R
deleted file mode 100644
index 96a10d1e..00000000
--- a/hardening-exactcl2.R
+++ /dev/null
@@ -1,141 +0,0 @@
-# Scoped exact-CL2 hardening gate (DESIGN amendment A18).
-#
-# Run from the package root with:
-# OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 nice -n 15 Rscript hardening-exactcl2.R
-# The script writes hardening-exactcl2-results.csv in the package root.
-
-make_exactcl2_fixture <- function(family_name, k, design) {
-  set.seed(
-    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
-      100 +
-      k +
-      match(design, c("unbalanced", "influential"))
-  )
-  n_curve <- 12L
-  n_grid <- 10L
-  grid <- seq(0, 1, length.out = n_grid)
-  xlin <- stats::rnorm(n_curve)
-  eta <- outer(xlin, rep(0.35, n_grid)) +
-    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
-    0.4
-  Y <- switch(
-    family_name,
-    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
-    poisson = matrix(
-      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
-      n_curve
-    ),
-    binomial = matrix(
-      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
-      n_curve
-    ),
-    Gamma = matrix(
-      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
-      n_curve
-    )
-  )
-  # Keep the response in its ordinary range, then make one covariate curve
-  # extreme. This isolates high leverage from a separable/outcome-saturation
-  # failure and exercises the shipped eigenvalue floor (especially Gaussian,
-  # Poisson, and Gamma at k = 9).
-  if (design == "influential") xlin[1] <- 3000
-  # Five deliberately unbalanced independent clusters, with cluster 1 made
-  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
-  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
-  list(
-    data = data.frame(Y = I(Y), xlin = xlin),
-    yind = grid,
-    cluster = cluster,
-    family = switch(
-      family_name,
-      gaussian = stats::gaussian(),
-      poisson = stats::poisson(),
-      binomial = stats::binomial(),
-      Gamma = stats::Gamma(link = "log")
-    ),
-    k = k
-  )
-}
-
-run_exactcl2_hardening <- function() {
-  cases <- expand.grid(
-    family = c("gaussian", "poisson", "binomial", "Gamma"),
-    k = c(4L, 7L),
-    design = c("unbalanced", "influential"),
-    stringsAsFactors = FALSE
-  )
-  rows <- vector("list", nrow(cases))
-
-  for (i in seq_len(nrow(cases))) {
-    case <- cases[i, ]
-    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
-    fit <- suppressWarnings(suppressMessages(refund::pffr(
-      Y ~ xlin,
-      data = fixture$data,
-      yind = fixture$yind,
-      family = fixture$family,
-      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
-      sandwich = "none"
-    )))
-
-    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
-    # cl2_exact_nocap comparator without duplicating the package implementation.
-    V_exact <- suppressWarnings(refund:::pffr_vcov(
-      fit,
-      sandwich = "cl2",
-      cluster = fixture$cluster,
-      cl2_adjustment = "exact"
-    ))
-    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
-      refund:::pffr_model_based_gam(fit),
-      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
-      cl2_adjustment = "exact",
-      leverage_cap = 1 - sqrt(1e-8)
-    ))
-    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
-      fit,
-      sandwich = "cl2",
-      cluster = fixture$cluster,
-      cl2_adjustment = "shortcut"
-    ))
-
-    se_exact <- sqrt(pmax(diag(V_exact), 0))
-    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
-    positive <- is.finite(se_exact) &
-      is.finite(se_shortcut) &
-      se_exact > 0 &
-      se_shortcut > 0
-    ratio <- se_exact[positive] / se_shortcut[positive]
-    rows[[i]] <- data.frame(
-      family = case$family,
-      k = case$k,
-      design = case$design,
-      G = length(unique(fixture$cluster)),
-      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
-      n_adjusted = attr(V_exact, "n_adjusted"),
-      min_block_eig = attr(V_exact, "min_block_eig"),
-      max_block_kappa = attr(V_exact, "max_block_kappa"),
-      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
-      median_se_ratio_to_shortcut = stats::median(ratio),
-      se_finite_positive = all(positive),
-      se_ratio_sane = length(ratio) > 0 &&
-        min(ratio) > 1e-3 &&
-        max(ratio) < 1e3,
-      stringsAsFactors = FALSE
-    )
-  }
-  do.call(rbind, rows)
-}
-
-if (sys.nframe() == 0L) {
-  devtools::load_all(quiet = TRUE)
-  results <- run_exactcl2_hardening()
-  utils::write.csv(results, "hardening-exactcl2-results.csv", row.names = FALSE)
-  print(results)
-  stopifnot(
-    all(results$se_finite_positive),
-    all(results$se_ratio_sane),
-    all(is.finite(results$n_adjusted)),
-    any(results$n_adjusted > 0)
-  )
-}

From f8a8fa6f391b5ded7bfbf2b15d34d3a6f5d14867 Mon Sep 17 00:00:00 2001
From: fabian-s <fabian.scheipl@googlemail.com>
Date: Tue, 4 Aug 2026 11:24:31 +0200
Subject: [PATCH 6/6] CL2: detect a penalized hat that has broken its own
 bounds (LB P-LB5)

Study LB recorded Poisson replicates whose cluster-robust interval widths were
numerically exploded -- up to 1e133 -- while the fit was recorded as a success,
so no failure gate saw them. Forensics: those replicates are ones where the
simulated response ran many orders of magnitude above its nominal mean, leaving
a fit whose model-based bread is no longer the inverse of the same weighted
cross-product. The penalized hat then violates bounds that hold exactly
(0 <= h_ii <= 1, 0 <= eigen(H_gg) <= 1): per-cluster hat eigenvalues up to
1.8e131 were recorded. CR1 explodes there too, so this is not a defect of the
leverage adjustment.

gam_sandwich_cluster_cl2() now monitors those invariants -- both monitors are
free, h_ii falls out of T_g and the block eigenvalues are computed anyway --
and warns when they break by more than round-off, exposing max_obs_leverage,
min_block_eig_rel and hat_invariant_violation as attributes.

It warns rather than substituting: capping H_gg silently turns an impossible
leverage into the largest legitimate one, so the shortcut answers garbage with
its maximal inflation while the exact block, PSD by construction, deflates.
Both are bounded, neither is meaningful, and falling back to CR1 would only
swap in a covariance built from the same broken bread. Measured on the new
fixture at max(Y) = 8.3e17: CR1 5.2e85, shortcut 1.7e87, exact 8.7e8 -- the
exact path is far less explosive, but not trustworthy either.

Also restores tests/testthat/test-pffr-exactcl2.R, which has been unable to run
since 276ba717 removed the top-level script it source()d; its fixtures now live
in tests/testthat/helper-exactcl2.R, together with the new P-LB5 fixture.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FaZwk9uigjyiVkuTTigSaH
---
 NEWS.md                             |  14 +++
 R/pffr-core.R                       | 126 +++++++++++++++++++++-
 man/gam_sandwich_cluster_cl2.Rd     |   7 +-
 man/pffr_hat_invariant_violation.Rd |  54 ++++++++++
 tests/testthat/helper-exactcl2.R    | 158 ++++++++++++++++++++++++++++
 tests/testthat/test-pffr-exactcl2.R |  84 ++++++++++++++-
 6 files changed, 437 insertions(+), 6 deletions(-)
 create mode 100644 man/pffr_hat_invariant_violation.Rd
 create mode 100644 tests/testthat/helper-exactcl2.R

diff --git a/NEWS.md b/NEWS.md
index a0e4f25f..aaff3adf 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -57,6 +57,20 @@
   `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
   explicitly selects either variant. Exact-CL2 diagnostics include the number
   of eigenvalue-floored blocks (`n_adjusted`).
+* The CL2 sandwich now checks the penalized hat matrix against its own bounds
+  (`0 <= h_ii <= 1`, `0 <= eigen(H_gg) <= 1`, and positive semi-definiteness of
+  the exact Bell--McCaffrey block) and **warns** when they are violated by more
+  than round-off. Such a violation means the model-based bread and the weighted
+  design have become numerically inconsistent -- a barely converged or
+  extremely ill-conditioned fit -- so the cluster-robust covariance is
+  meaningless however plausible it looks. Previously this produced silently
+  exploded standard errors (interval widths up to 1e133 were observed on
+  degenerate Poisson fits, with nothing to distinguish them from a legitimately
+  wide interval). The covariance is still returned, now carrying
+  `max_obs_leverage` and `hat_invariant_violation` attributes. The CR1
+  (`sandwich = "cluster"`) covariance is built from the same bread and is
+  equally affected, so switching sandwich type is not a remedy: inspect and
+  refit the model.
 * AR(1) support improvements: `pffr()` now automatically switches to
   `algorithm = "bam"` and `method = "fREML"` when `rho` is supplied, and
   sets `discrete = TRUE` for non-Gaussian families.
diff --git a/R/pffr-core.R b/R/pffr-core.R
index d3712dc5..354ab921 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -1631,8 +1631,11 @@ gam_sandwich_cluster <- function(
 #'   (X6).
 #' @returns A p x p covariance matrix with leverage diagnostics. Exact CL2
 #'   returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
-#'   `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
-#'   `n_capped_clusters` and `max_leverage` attributes.
+#'   `min_block_eig`, `max_block_kappa` and `min_block_eig_rel`; the shortcut
+#'   returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both
+#'   return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the
+#'   penalized hat respects its bounds, otherwise a description of the
+#'   violation; see [pffr_hat_invariant_violation()]).
 #' @keywords internal
 gam_sandwich_cluster_cl2 <- function(
   b,
@@ -1706,6 +1709,10 @@ gam_sandwich_cluster_cl2 <- function(
   n_adjusted <- 0L
   min_block_eig <- Inf
   max_block_kappa <- 0
+  # Hat-invariant monitors (see pffr_hat_invariant_violation()). Both are free:
+  # h_ii is a by-product of T_g, and the block eigenvalues are computed anyway.
+  max_obs_leverage <- NA_real_
+  min_block_eig_rel <- 0
 
   for (g in groups) {
     idx <- which(cluster_id_work == g)
@@ -1716,6 +1723,14 @@ gam_sandwich_cluster_cl2 <- function(
     Hgg <- Tg %*% t(Xwg)
     Hgg <- 0.5 * (Hgg + t(Hgg))
 
+    # Per-observation leverage h_ii = diag(H)_ii, mathematically in [0, 1] for
+    # the penalized hat. Tracked in both branches; a value above 1 means the
+    # bread and the weighted design are numerically inconsistent.
+    h_ii <- rowSums(Tg * Xwg)
+    if (any(is.finite(h_ii))) {
+      max_obs_leverage <- max(max_obs_leverage, h_ii, na.rm = TRUE)
+    }
+
     if (cl2_adjustment == "shortcut") {
       ee_H <- eigen(Hgg, symmetric = TRUE)
       cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
@@ -1743,6 +1758,14 @@ gam_sandwich_cluster_cl2 <- function(
       ee_B <- eigen(0.5 * (Bg + t(Bg)), symmetric = TRUE)
       block_min <- min(ee_B$values)
       min_block_eig <- min(min_block_eig, block_min)
+      # B_g = ((I - H)^2)_gg is a principal block of a squared symmetric
+      # matrix and is therefore positive semi-definite in exact arithmetic.
+      # A materially negative eigenvalue, measured relative to the block's own
+      # scale, is an invariant violation rather than round-off.
+      min_block_eig_rel <- min(
+        min_block_eig_rel,
+        block_min / max(abs(max(ee_B$values)), 1e-300)
+      )
       max_block_kappa <- max(
         max_block_kappa,
         abs(max(ee_B$values)) / max(abs(block_min), 1e-300)
@@ -1776,6 +1799,30 @@ gam_sandwich_cluster_cl2 <- function(
   attr(V, "min_block_eig") <- if (is.finite(min_block_eig)) min_block_eig else
     NA_real_
   attr(V, "max_block_kappa") <- max_block_kappa
+  attr(V, "max_obs_leverage") <- max_obs_leverage
+  attr(V, "min_block_eig_rel") <- min_block_eig_rel
+
+  hat_violation <- pffr_hat_invariant_violation(
+    max_obs_leverage = max_obs_leverage,
+    max_leverage = if (cl2_adjustment == "shortcut") max_leverage else NA_real_,
+    min_block_eig_rel = if (cl2_adjustment == "exact") min_block_eig_rel else 0
+  )
+  attr(V, "hat_invariant_violation") <- hat_violation
+  if (!is.null(hat_violation)) {
+    warning(
+      "Cluster-robust covariance is NOT trustworthy for this fit: ",
+      hat_violation,
+      " The penalized hat matrix satisfies 0 <= h_ii <= 1 and ",
+      "0 <= eigen(H_gg) <= 1 exactly, so this indicates that the model-based ",
+      "bread and the weighted design have become numerically inconsistent -- ",
+      "typically an ill-conditioned or barely converged fit (extreme fitted ",
+      "values, a near-singular penalized Hessian, or a basis far too rich for ",
+      "the data). Both the CL2 and the CR1 cluster-robust covariances inherit ",
+      "the problem, so switching sandwich type does not repair it: inspect and ",
+      "refit the model instead.",
+      call. = FALSE
+    )
+  }
   if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) {
     max_leverage_label <- if (is.finite(max_leverage)) {
       sprintf("%.3f", max_leverage)
@@ -1801,6 +1848,81 @@ gam_sandwich_cluster_cl2 <- function(
   V
 }
 
+#' Numerical-sanity check on the per-cluster hat blocks
+#'
+#' The penalized hat \eqn{H = X_w V_p X_w'} obeys \eqn{0 \le h_{ii} \le 1} and,
+#' since any principal block of a symmetric matrix has its eigenvalues inside
+#' the parent's range, \eqn{0 \le \mathrm{eigen}(H_{gg}) \le 1}. The exact
+#' Bell--McCaffrey block \eqn{B_g = ((I - H)^2)_{gg}} is positive semi-definite
+#' for the same reason. When a fit is numerically degenerate --- an
+#' ill-conditioned or barely converged fit, e.g. a Poisson fit whose fitted
+#' means span many orders of magnitude --- the bread \eqn{V_p} stops being the
+#' inverse of the same weighted cross-product, and these invariants break by
+#' orders of magnitude rather than by round-off. This was the mechanism behind
+#' the numerically exploded interval widths observed on hard Poisson
+#' replicates in the locked benchmark (study LB, claim P-LB5).
+#'
+#' Detection is deliberately separate from the leverage cap. Capping
+#' \eqn{H_{gg}} at `leverage_cap` silently converts an impossible leverage into
+#' the *largest legitimate* one, so the shortcut adjustment responds to garbage
+#' with its maximal variance inflation; the exact block, being a squared
+#' operator, instead deflates. Both are bounded, neither is meaningful, and
+#' only a diagnostic can tell the user which situation they are in. We
+#' therefore warn and keep going rather than substituting a differently wrong
+#' number: the CR1 covariance is built from the same bread and is no more
+#' trustworthy here.
+#'
+#' @param max_obs_leverage Largest per-observation leverage \eqn{h_{ii}} seen,
+#'   or `NA`.
+#' @param max_leverage Largest eigenvalue of any \eqn{H_{gg}} (shortcut path
+#'   only; `NA` otherwise).
+#' @param min_block_eig_rel Smallest eigenvalue of any \eqn{B_g} relative to
+#'   that block's largest eigenvalue (exact path only; `0` otherwise).
+#' @param tol Relative slack allowed before an invariant counts as violated.
+#' @returns `NULL` when every invariant holds, otherwise a one-sentence
+#'   character description of the violation.
+#' @keywords internal
+pffr_hat_invariant_violation <- function(
+  max_obs_leverage = NA_real_,
+  max_leverage = NA_real_,
+  min_block_eig_rel = 0,
+  tol = 1e-6
+) {
+  msgs <- character(0)
+  if (isTRUE(is.finite(max_obs_leverage) && max_obs_leverage > 1 + tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        "the largest per-observation leverage is %.3g, above the bound 1;",
+        max_obs_leverage
+      )
+    )
+  }
+  if (isTRUE(is.finite(max_leverage) && max_leverage > 1 + tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        "the largest per-cluster hat eigenvalue is %.3g, above the bound 1;",
+        max_leverage
+      )
+    )
+  }
+  if (isTRUE(is.finite(min_block_eig_rel) && min_block_eig_rel < -tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        paste0(
+          "the exact Bell-McCaffrey block has a relative eigenvalue of %.3g, ",
+          "but it is positive semi-definite by construction;"
+        ),
+        min_block_eig_rel
+      )
+    )
+  }
+  if (length(msgs) == 0) return(NULL)
+  paste(msgs, collapse = " ")
+}
+
 #' Resolve the CL2 leverage adjustment
 #'
 #' Exact CL2 is the default when the finite-sample correction can matter
diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd
index 3bd0e45f..ee0668f0 100644
--- a/man/gam_sandwich_cluster_cl2.Rd
+++ b/man/gam_sandwich_cluster_cl2.Rd
@@ -45,8 +45,11 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 \value{
 A p x p covariance matrix with leverage diagnostics. Exact CL2
   returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
-  `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy
-  `n_capped_clusters` and `max_leverage` attributes.
+  `min_block_eig`, `max_block_kappa` and `min_block_eig_rel`; the shortcut
+  returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both
+  return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the
+  penalized hat respects its bounds, otherwise a description of the
+  violation; see [pffr_hat_invariant_violation()]).
 }
 \description{
 Computes a cluster-robust covariance matrix with a Bell-McCaffrey style
diff --git a/man/pffr_hat_invariant_violation.Rd b/man/pffr_hat_invariant_violation.Rd
new file mode 100644
index 00000000..3c0ed713
--- /dev/null
+++ b/man/pffr_hat_invariant_violation.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pffr-core.R
+\name{pffr_hat_invariant_violation}
+\alias{pffr_hat_invariant_violation}
+\title{Numerical-sanity check on the per-cluster hat blocks}
+\usage{
+pffr_hat_invariant_violation(
+  max_obs_leverage = NA_real_,
+  max_leverage = NA_real_,
+  min_block_eig_rel = 0,
+  tol = 1e-06
+)
+}
+\arguments{
+\item{max_obs_leverage}{Largest per-observation leverage \eqn{h_{ii}} seen,
+or `NA`.}
+
+\item{max_leverage}{Largest eigenvalue of any \eqn{H_{gg}} (shortcut path
+only; `NA` otherwise).}
+
+\item{min_block_eig_rel}{Smallest eigenvalue of any \eqn{B_g} relative to
+that block's largest eigenvalue (exact path only; `0` otherwise).}
+
+\item{tol}{Relative slack allowed before an invariant counts as violated.}
+}
+\value{
+`NULL` when every invariant holds, otherwise a one-sentence
+  character description of the violation.
+}
+\description{
+The penalized hat \eqn{H = X_w V_p X_w'} obeys \eqn{0 \le h_{ii} \le 1} and,
+since any principal block of a symmetric matrix has its eigenvalues inside
+the parent's range, \eqn{0 \le \mathrm{eigen}(H_{gg}) \le 1}. The exact
+Bell--McCaffrey block \eqn{B_g = ((I - H)^2)_{gg}} is positive semi-definite
+for the same reason. When a fit is numerically degenerate --- an
+ill-conditioned or barely converged fit, e.g. a Poisson fit whose fitted
+means span many orders of magnitude --- the bread \eqn{V_p} stops being the
+inverse of the same weighted cross-product, and these invariants break by
+orders of magnitude rather than by round-off. This was the mechanism behind
+the numerically exploded interval widths observed on hard Poisson
+replicates in the locked benchmark (study LB, claim P-LB5).
+}
+\details{
+Detection is deliberately separate from the leverage cap. Capping
+\eqn{H_{gg}} at `leverage_cap` silently converts an impossible leverage into
+the *largest legitimate* one, so the shortcut adjustment responds to garbage
+with its maximal variance inflation; the exact block, being a squared
+operator, instead deflates. Both are bounded, neither is meaningful, and
+only a diagnostic can tell the user which situation they are in. We
+therefore warn and keep going rather than substituting a differently wrong
+number: the CR1 covariance is built from the same bread and is no more
+trustworthy here.
+}
+\keyword{internal}
diff --git a/tests/testthat/helper-exactcl2.R b/tests/testthat/helper-exactcl2.R
new file mode 100644
index 00000000..711e7bd2
--- /dev/null
+++ b/tests/testthat/helper-exactcl2.R
@@ -0,0 +1,158 @@
+# Fixtures for the exact-CL2 hardening tests (DESIGN amendments A18 and the
+# 2026-08-04 study-LB P-LB5 follow-up). The full standalone gates, which also
+# write the archived result tables, live in the pffr-ci repo
+# (diagnostics/hardening-exactcl2.R); these helpers keep the same coverage
+# available to CI without a non-standard top-level file in the package.
+
+make_exactcl2_fixture <- function(family_name, k, design) {
+  set.seed(
+    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
+      100 +
+      k +
+      match(design, c("unbalanced", "influential"))
+  )
+  n_curve <- 12L
+  n_grid <- 10L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  eta <- outer(xlin, rep(0.35, n_grid)) +
+    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
+    0.4
+  Y <- switch(
+    family_name,
+    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
+    poisson = matrix(
+      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
+      n_curve
+    ),
+    binomial = matrix(
+      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
+      n_curve
+    ),
+    Gamma = matrix(
+      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
+      n_curve
+    )
+  )
+  # Keep the response in its ordinary range, then make one covariate curve
+  # extreme. This isolates high leverage from a separable/outcome-saturation
+  # failure and exercises the shipped eigenvalue floor.
+  if (design == "influential") xlin[1] <- 3000
+  # Five deliberately unbalanced independent clusters, with cluster 1 made
+  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
+  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
+  list(
+    data = data.frame(Y = I(Y), xlin = xlin),
+    yind = grid,
+    cluster = cluster,
+    family = switch(
+      family_name,
+      gaussian = stats::gaussian(),
+      poisson = stats::poisson(),
+      binomial = stats::binomial(),
+      Gamma = stats::Gamma(link = "log")
+    ),
+    k = k
+  )
+}
+
+run_exactcl2_hardening <- function() {
+  cases <- expand.grid(
+    family = c("gaussian", "poisson", "binomial", "Gamma"),
+    k = c(4L, 7L),
+    design = c("unbalanced", "influential"),
+    stringsAsFactors = FALSE
+  )
+  rows <- vector("list", nrow(cases))
+
+  for (i in seq_len(nrow(cases))) {
+    case <- cases[i, ]
+    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
+    fit <- suppressWarnings(suppressMessages(refund::pffr(
+      Y ~ xlin,
+      data = fixture$data,
+      yind = fixture$yind,
+      family = fixture$family,
+      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+      sandwich = "none"
+    )))
+
+    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
+    # cl2_exact_nocap comparator without duplicating the implementation.
+    V_exact <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "exact"
+    ))
+    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
+      refund:::pffr_model_based_gam(fit),
+      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
+      cl2_adjustment = "exact",
+      leverage_cap = 1 - sqrt(1e-8)
+    ))
+    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "shortcut"
+    ))
+
+    se_exact <- sqrt(pmax(diag(V_exact), 0))
+    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
+    positive <- is.finite(se_exact) &
+      is.finite(se_shortcut) &
+      se_exact > 0 &
+      se_shortcut > 0
+    ratio <- se_exact[positive] / se_shortcut[positive]
+    rows[[i]] <- data.frame(
+      family = case$family,
+      k = case$k,
+      design = case$design,
+      G = length(unique(fixture$cluster)),
+      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
+      n_adjusted = attr(V_exact, "n_adjusted"),
+      min_block_eig = attr(V_exact, "min_block_eig"),
+      max_block_kappa = attr(V_exact, "max_block_kappa"),
+      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
+      median_se_ratio_to_shortcut = stats::median(ratio),
+      se_finite_positive = all(positive),
+      se_ratio_sane = length(ratio) > 0 &&
+        min(ratio) > 1e-3 &&
+        max(ratio) < 1e3,
+      stringsAsFactors = FALSE
+    )
+  }
+  do.call(rbind, rows)
+}
+
+#' Study-LB P-LB5 fixture: a Poisson design whose latent roughness amplitude
+#' `amp` pushes individual fitted means far above the nominal marginal mean of
+#' 4, the mechanism behind the exploded interval widths in the locked
+#' benchmark. `amp = 1` is the benign control; `amp = 10` is degenerate.
+make_lb5_fixture <- function(amp, n_grid, k, seed) {
+  set.seed(seed)
+  n_curve <- 8L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  nfpc <- 6L
+  rough <- matrix(stats::rnorm(n_curve * nfpc), n_curve, nfpc) %*%
+    t(sapply(seq_len(nfpc), function(j) sin(j * pi * grid)))
+  eta <- outer(xlin, rep(0.35, n_grid)) + amp * rough + log(4)
+  Y <- matrix(
+    stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 600))),
+    n_curve
+  )
+  list(data = data.frame(Y = I(Y), xlin = xlin), yind = grid, k = k)
+}
+
+fit_lb5_fixture <- function(fixture) {
+  suppressWarnings(suppressMessages(refund::pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    family = stats::poisson(),
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "none"
+  )))
+}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
index d4faf7af..d9b88b10 100644
--- a/tests/testthat/test-pffr-exactcl2.R
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -1,7 +1,8 @@
 # A18 exact-CL2 hardening gate. The standalone script writes the same table
 # for release records; this fast test keeps its family/design coverage in CI.
-
-source(testthat::test_path("..", "..", "hardening-exactcl2.R"), local = TRUE)
+# Fixtures live in helper-exactcl2.R (they used to be sourced from a top-level
+# script that was removed from the package tree, which left this file unable to
+# run at all).
 
 test_that("exact CL2 is robust on high-leverage family and basis fixtures", {
   results <- run_exactcl2_hardening()
@@ -61,3 +62,82 @@ test_that("automatic exact CL2 uses the documented relevance and cost rule", {
   expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
   expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
 })
+
+# ---------------------------------------------------------------------------
+# Study-LB P-LB5 follow-up (2026-08-04): numerically exploded interval widths
+# on degenerate Poisson fits must be DETECTED, not returned silently.
+# ---------------------------------------------------------------------------
+
+test_that("pffr_hat_invariant_violation() flags only impossible hat values", {
+  viol <- refund:::pffr_hat_invariant_violation
+
+  # Every invariant holds -> NULL.
+  expect_null(viol(max_obs_leverage = 0.99, max_leverage = 0.999))
+  expect_null(viol(max_obs_leverage = NA_real_, max_leverage = NA_real_))
+  # Round-off-sized excess is tolerated.
+  expect_null(viol(max_obs_leverage = 1 + 1e-9))
+  # Real violations are described, with the offending number in the message.
+  expect_match(
+    viol(max_obs_leverage = 22.6),
+    "per-observation leverage is 22.6"
+  )
+  expect_match(viol(max_leverage = 8.18e84), "hat eigenvalue")
+  expect_match(viol(min_block_eig_rel = -0.3), "positive semi-definite")
+  # A non-finite monitor must not trip the check.
+  expect_null(viol(max_obs_leverage = Inf, max_leverage = NaN))
+})
+
+test_that("a benign Poisson CL2 fit reports leverages inside the bounds", {
+  skip_on_cran()
+  fit <- fit_lb5_fixture(make_lb5_fixture(amp = 1, n_grid = 20L, k = 8L, 21L))
+  V_exact <- expect_no_warning(
+    refund:::pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "exact")
+  )
+  V_short <- expect_no_warning(
+    refund:::pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "shortcut")
+  )
+  expect_lte(attr(V_exact, "max_obs_leverage"), 1)
+  expect_lte(attr(V_short, "max_leverage"), 1)
+  expect_null(attr(V_exact, "hat_invariant_violation"))
+  expect_null(attr(V_short, "hat_invariant_violation"))
+  # exact and shortcut agree closely in the well-behaved regime
+  ratio <- sqrt(max(abs(diag(V_exact)))) / sqrt(max(abs(diag(V_short))))
+  expect_gt(ratio, 0.5)
+  expect_lt(ratio, 2)
+})
+
+test_that("a degenerate Poisson fit warns instead of returning a silent 1e30", {
+  skip_on_cran()
+  # amp = 10 drives individual fitted means many orders of magnitude above the
+  # nominal marginal mean, the study-LB mechanism; the penalized hat then
+  # breaks its own bounds by dozens of orders of magnitude.
+  fit <- fit_lb5_fixture(make_lb5_fixture(amp = 10, n_grid = 30L, k = 12L, 21L))
+
+  expect_warning(
+    V_exact <- refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cl2_adjustment = "exact"
+    ),
+    "NOT trustworthy"
+  )
+  expect_warning(
+    V_short <- refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cl2_adjustment = "shortcut"
+    ),
+    "NOT trustworthy"
+  )
+  expect_gt(attr(V_exact, "max_obs_leverage"), 1)
+  expect_type(attr(V_exact, "hat_invariant_violation"), "character")
+
+  se_exact <- sqrt(max(abs(diag(V_exact))))
+  se_short <- sqrt(max(abs(diag(V_short))))
+  expect_true(is.finite(se_exact) && is.finite(se_short))
+  # The exact block is a principal block of a squared symmetric matrix and so
+  # stays positive semi-definite; the shortcut's I - H_gg turns indefinite and
+  # is floored into the largest legitimate variance inflation. The exact path
+  # is therefore far less explosive here, though neither is meaningful.
+  expect_lt(se_exact, se_short)
+})