Geometric Robustness of Sampling (Block RG)

Diagnosing posterior geometry and climbing a hierarchy of geometry-adaptive samplers

José Mauricio Gómez Julián

2026-07-06

What this vignette is. A complete, abundantly documented narrative of Block RG (Robustness of the Geometry of Sampling) of gdpar: why it exists, every design hinge and the decision taken at it, every validation, and how to use the resulting capability. It is opt-in: nothing here changes the default gdpar() fit, whose draws stay bit-identical.

Two levels of reproducibility. The chunks evaluated when this vignette is built are deliberately cheap – they inspect the ladder, the budget, the success criteria and the calibrated thresholds, and they run the geometry engine on small closure targets in pure R, with no Stan compilation. The heavy, cmdstan-backed runs (the diagnostic on a compiled model, the full orchestrator, the bridge over a real fit) are shown with eval = FALSE for reading, and reproduced end to end by the gated script inst/scripts/geometry_pilots_deep.R (set GDPAR_RUN_GEOMETRY_PILOTS=1), referenced throughout. This keeps the package build light without hiding the real computation.

Why this block exists

The block was opened (session B9.21) by a concrete failure. In the external re-validation 9.2.O, the count coordinate of a bivariate eBird model – a Tweedie outcome – would not converge under the no-U-turn Hamiltonian explorer (“NUTS”): split-\(\hat R \approx 3\), effective sample size \(\approx 2\), nonsensical predictive accuracy, more than an hour per fit at high warm-up. A three-layer forensic (B9.21) separated two things:

  1. A false culprit, corrected. The infinite \(\hat R\) was produced by per-slot additive scales sigma_a_k that were declared but unused for distributional slots without an a() term (phi ~ 1, p ~ 1): a flat direction (non-identification), a genuine geometric pathology. The fix (“Option A”, session B9.30) compacts those scales away (Section @ref(optiona)); it cleans the infinite \(\hat R\) but does not unstick the Tweedie.
  2. The real cause: the intrinsic geometry of the posterior. The eBird count is almost entirely explained by geography and climate, so the posterior is near-deterministic – the typical set collapses toward a lower-dimensional manifold (a razor-thin canyon), aggravated by the Tweedie shape parameter \(p\) pinned against its bound \((1.01, 1.99)\). The competitor mgcv wins because it does not traverse anything (REML plus a saddlepoint); a Hamiltonian sampler does traverse, and there it jams.

The user’s decision was not to patch the case but to give gdpar first-class geometric robustness: a capability that (a) diagnoses the geometry of a posterior and (b) offers geometry-adaptive sampling that climbs a hierarchy of geometries – Riemannian and beyond, not stopping at the Riemannian. Opt-in, never the default. “Let us go for the glory.”

The conceptual bridge (geometry of the user’s document ↔︎ sampling)

The block is grounded in the user’s manuscript ON_RIEMMANIAN_STATISTICS.Rmd, which supplies the ontological-geometric frame and a hierarchy of geometries (its Part IV, §12). Hamiltonian sampling simulates a conservative dynamics over the relief of the log-density; its efficiency depends on how well the metric (“mass” of the system) fits the local shape of that relief. This is applied Riemannian geometry:

An organic critical reading of the source documents

Following the package’s standing posture (the user’s manuscripts are input to improve in our context, not holy writ), three points are worth stating plainly.

  • §14 of the manuscript (“On the use of statistical models”) is programmatic: it declares that “the systematic development of the statistical implications of this framework constitutes a future line of research.” Block RG is the concrete operationalization of that line: it turns §12’s hierarchy of geometries into a working diagnostic and a working sampler ladder, with the Fisher information as the §5 local inner product and the Hamiltonian trajectories as the §6 geodesics.
  • We realize sub-Riemannian geometry faithfully but distinctly. The manuscript motivates it through non-holonomic reachability (the Chow–Rashevsky theorem, §12.4: restricted instantaneous directions, yet total reachability). In sampling we read the accessible distribution \(D_\theta\) from the near-null space of the expected Fisher and integrate the stiff “walls” by their exact Gaussian flow; the manuscript’s reachability intuition maps to our “typical set on a thinner manifold”, but the construction is our own (and Metropolis-exact, so reachability of the full space is preserved by the correction, not assumed).
  • We deliberately stop where the manuscript keeps going. Its §12.5 (Alexandrov / RCD metric spaces) and §12.6 (fractal geometry, no tangent space) describe regimes our differentiable machinery does not cover. We treat multimodality and non-smoothness as out-of-ladder remedies (tempering, a general-metric fallback) rather than pretend the geometry ladder solves them (Section @ref(orchestrator)). Honesty about the boundary is part of the design.

The honesty convention itself we borrow from the companion manuscript A-ORPHEUS-PIMC.Rmd (§16.3): keep demonstrated (correctness) strictly apart from conjectured (efficiency). Every sampler level below is Metropolis-exact, so which geometry is chosen governs only efficiency, never the validity of the draws. The speed-ups are measured, never asserted.

The taxonomy: pathologies as a hierarchy of geometries

“Hard geometry” is not a scalar but a taxonomy, each entry tied to a geometric concept and to a remedy (a level of the geometry hierarchy):

Pathology Geometric concept Level / remedy
Isotropy (easy) identity metric Euclidean diagonal (default)
Anisotropy / elongation (straight canyon) constant non-identity metric Euclidean dense
Funnel / variable curvature position-dependent metric Riemannian (Fisher / SoftAbs)
Heavy tails / directional asymmetry non-inner-product norm Finsler / relativistic
Near-determinism → lower-dim manifold distribution \(D_p \subsetneq T_pM\) sub-Riemannian
Non-smoothness / multimodality general metric space out of ladder (tempering / fallback)
Flat direction (non-identification) zero Hessian eigenvalue reparametrize / eliminate (Option A)

The calibration of the diagnostic that maps a posterior to one of these classes is validated against a synthetic suite of geometries of known difficulty, the eight targets G0–G7 (RG.1). Building the suite is free (no Stan):

suite <- gdpar_geometry_suite()
data.frame(
  target  = names(suite),
  pathology = vapply(suite, `[[`, character(1), "pathology"),
  remedy    = vapply(suite, `[[`, character(1), "geometry_remedy")
)
#>                                        target           pathology
#> G0_isotropic                     G0_isotropic           isotropic
#> G1_anisotropic                 G1_anisotropic         anisotropic
#> G2_funnel                           G2_funnel              funnel
#> G3_heavy_tails                 G3_heavy_tails         heavy_tails
#> G4_quasi_deterministic G4_quasi_deterministic quasi_deterministic
#> G5_multimodal                   G5_multimodal          multimodal
#> G6_boundary                       G6_boundary            boundary
#> G7_flat_direction           G7_flat_direction      flat_direction
#>                                      remedy
#> G0_isotropic             euclidean_diagonal
#> G1_anisotropic              euclidean_dense
#> G2_funnel                        riemannian
#> G3_heavy_tails         finsler_relativistic
#> G4_quasi_deterministic       sub_riemannian
#> G5_multimodal                     tempering
#> G6_boundary                boundary_reparam
#> G7_flat_direction         reparam_eliminate

Each target ships in dual form – a Stan program (for cmdstan-backed pilots) and a pure-R closure (log-density and gradient, cross-checked against each other to \(< 10^{-9}\)). The closure side is what makes the engine demonstrations below run with no compilation.

The diagnostic (RG.1): size-invariant signals

The diagnostic runs cheap pilots and reads signals that are invariant to the sample size – divergences, the energy fraction of missing information (“E-BFMI”), tree depth, the condition number, the step-to-scale ratio – and a difficulty-versus-\(n\) curve that separates near-determinism (which grows with \(n\)) from a structural pathology (constant). It localizes the culprit parameter and classifies the pathology with a transparent, calibrated rule-based classifier.

The decisive lesson (B9.20/B9.21). The diagnostic never uses \(\hat R\) or the effective sample size on short runs. The infamous infinite \(\hat R\) at warm-up 100 was a false positive; the size-invariant signals are the honest ones. This lesson is enforced everywhere downstream, including the success gate of the orchestrator.

The classifier’s thresholds were calibrated in RG.1.c over a 540-cell grid (eight targets \(\times\) three difficulties \(\times\) two pilot budgets \(\times\) replicates) with a minimax adaptive allocation (a Wilson-interval top-up of the unresolved cells). Reporting honestly: the macro accuracy rose from \(0.78\) to \(0.91\) and the held-out balanced accuracy from \(0.60\) to \(0.89\); six of eight classes land at \(0.93\)–\(1.00\), while the funnel (\(0.63\)) and the heavy tail (\(0.71\)) remain mutually confusable – a real limit that the orchestrator’s level selection is designed to be robust to. The recalibrated thresholds are plain, inspectable data:

str(gdpar_geometry_thresholds())
#> List of 13
#>  $ divergent_rate_high : num 0.01
#>  $ funnel_ebfmi_low    : num 0.35
#>  $ heavy_cond_max      : num 25
#>  $ treedepth_sat_high  : num 0.2
#>  $ condition_high      : num 12
#>  $ step_scale_ratio_low: num 0.1
#>  $ nslope_grows        : num 0.8
#>  $ flat_var_high       : num 600
#>  $ boundary_prox_high  : num 0.02
#>  $ boundary_eps        : num 0.01
#>  $ multimodal_high     : num 2.5
#>  $ heavy_kurtosis_high : num 1.8
#>  $ target_ess          : num 400

Running the diagnostic on a real model is heavier (it compiles and samples), so it is shown but not evaluated here:

# Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
diag <- gdpar_geometry_diagnostic(suite$G4_quasi_deterministic, n_grid = 3)
diag$pathology      # "quasi_deterministic"
diag$culprit        # the localized direction(s)

The engine and the ladder (RG.2–RG.4)

The motor decision (RG.2, decision A) is an integrator written in R over a Stan backend: the generalized leapfrog, SoftAbs, relativistic and sub-Riemannian trajectories live in pure R, delegating the log-density, gradient and Hessian to the compiled cmdstan model via $log_prob / $grad_log_prob / $hessian. Under the cornerstone rule (maximum multidimensional robustness, cost irrelevant) this wins: “slow” is only constant-factor R orchestration, opt-in and reversible to compiled code, with zero new external dependency. The same engine runs on the suite’s closures, which is what the demonstrations below exploit.

The ladder the orchestrator can climb is {0, 1, 3, 4, 5}:

Level 0 – Euclidean diagonal (the default)

The identity-mass HMC. On a closure target it is pure R. An anisotropic Gaussian (precision \(\mathrm{diag}(1, 100)\), scales \(1\) and \(0.1\)) is recovered well when the elongation is mild:

P   <- diag(c(1, 100))
tgt <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * drop(t %*% P %*% t),
  grad_log_prob = function(t) -drop(P %*% t), dim = 2L)

fit0 <- gdpar_geom_hmc(tgt, gdpar_geom_metric_euclidean(dim = 2L),
                       n_iter = 400L, n_warmup = 200L, epsilon = 0.12,
                       L = 20L, seed = 1L)
c(accept = round(fit0$accept_rate, 3),
  sd1 = round(sd(fit0$draws[, 1]), 3),   # truth 1
  sd2 = round(sd(fit0$draws[, 2]), 3))   # truth 0.1
#> accept    sd1    sd2 
#>  0.912  0.943  0.108

Level 1 – Euclidean dense

A constant non-identity mass (a linear preconditioner) straightens a tilted canyon. It is the right remedy when the elongation is severe but the curvature does not vary across the space. Construct it from a mass matrix:

# A dense mass equal to the posterior precision whitens a straight canyon.
metric_dense <- gdpar_geom_metric_euclidean(M = P)

Level 3 – Riemannian (Fisher / SoftAbs, and the learned GP-Fisher)

When curvature varies with position (a funnel) a constant metric cannot keep up; the metric must depend on \(\theta\). Block RG builds this in three complementary ways (RG.3):

The integrator is the generalized implicit leapfrog of Girolami–Calderhead; the metric is a preconditioner, not part of the target, so the Metropolis correction with the exact density is the intrinsic repair – no delayed acceptance is needed (a deliberate improvement over the ORPHEUS surrogate, where the surrogate enters the acceptance). The funnel is the canonical heavy demonstration (Section @ref(appendix)).

Level 4 – Finsler / relativistic (heavy tails)

A bounded, non-Gaussian kinetic energy coupled to the position-dependent Riemannian metric, \[K(\theta, p) = c\sqrt{p^\top M(\theta)^{-1}p + m^2c^2} + \tfrac12\log\det M(\theta),\] whose velocity has \(M\)-norm strictly below the speed \(c\) – this is what tames the overshoot of heavy tails. The \(\tfrac12\log\det M\) normalizer keeps the marginal exact for every \(c\) and \(m\) (they govern only efficiency), and the non-relativistic limit \(c \to \infty\) recovers the Riemannian kinetic of Level 3. Because the resulting Hamiltonian is non-separable, a dedicated generalized implicit integrator handles it, routed in opt-in so the default leapfrog stays bit-identical. The honest Finsler reading: this kinetic is the Legendre dual of a Finsler norm; the asymmetric Randers extension \(F=\sqrt{g(v,v)}+\beta(v)\) is deliberately deferred because it is odd in \(p\) and would break the reversibility Metropolis depends on (it models irreversible dynamics, a different goal than sampling a fixed target).

Level 5 – sub-Riemannian (the count; “the glory”)

The remedy for the near-deterministic posterior. The accessible distribution \(D_\theta\) is read from the near-null space of the expected Fisher; a continuous spectral filter (no hard cut-off) splits each direction into a soft “floor” and a stiff “wall”; a Strang symplectic splitting integrates the stiff wall by its exact harmonic flow (so the wall does not penalize the step at any \(\varepsilon\)) while the soft residual follows free drift, with the Metropolis correction on the true density keeping the sampling exact (the reference Gaussian is an internal device of the integrator). The default threshold is the floor scale \(\tau = \lambda_{\min}\), so the step is governed by the floor, never the walls. On a closure – a mild canyon with Fisher \(\mathrm{diag}(1, 50)\) – it runs in pure R and recovers both scales with no divergences at a step the Euclidean sampler could not take:

canyon <- gdpar_geom_target(
  log_prob      = function(th) -0.5 * (th[1]^2 + 50 * th[2]^2),
  grad_log_prob = function(th) -c(th[1], 50 * th[2]), dim = 2L)

metric_sr <- gdpar_geom_metric_subriemannian(
  canyon, fisher = function(th) diag(c(1, 50)))

fit_sr <- gdpar_geom_hmc(canyon, metric = metric_sr, epsilon = 0.5, L = 12L,
                         n_iter = 800L, n_warmup = 300L, seed = 3L)
c(accept = round(fit_sr$accept_rate, 3),
  divergent = fit_sr$n_divergent,
  floor_sd = round(sd(fit_sr$draws[, 1]), 3),   # truth 1
  wall_sd  = round(sd(fit_sr$draws[, 2]), 3))   # truth 1/sqrt(50) ~ 0.141
#>    accept divergent  floor_sd   wall_sd 
#>     0.953     0.000     1.059     0.136

This is the level RG.7 will point at the real Tweedie count, supplying the expected Fisher through gdpar_geom_fisher_simulator() over the compiled model.

The orchestrator (RG.5): diagnose → select → sample → re-diagnose → escalate or certify

gdpar_geom_orchestrate() closes the loop. It diagnoses the geometry, selects an entry level, samples, re-diagnoses, and either resolves or emits a certified limit. The budget and the multi-signal success gate are tunable plain data:

str(gdpar_geom_orchestrate_budget())
#> List of 16
#>  $ max_rounds         : int 8
#>  $ max_levels         : int 5
#>  $ probe_warmup       : int 150
#>  $ probe_iter         : int 150
#>  $ full_warmup        : int 500
#>  $ full_iter          : int 500
#>  $ epsilon            : num 0.25
#>  $ L                  : int 25
#>  $ tune_epsilon       : logi TRUE
#>  $ tune_iter          : int 60
#>  $ max_seconds        : num Inf
#>  $ max_seconds_per_fit: num Inf
#>  $ max_fits           : int 40
#>  $ n_rediagnose       : int 1
#>  $ stall_limit        : int 2
#>  $ hysteresis         : num 0.1
str(gdpar_geom_orchestrate_criteria())
#> List of 4
#>  $ accept_low         : num 0.5
#>  $ accept_high        : num 0.999
#>  $ divergent_rate_high: num 0.02
#>  $ ebfmi_low          : num 0.3

The four design hinges (decided by the user, B9.29) and how they are honoured:

Crucially, the out-of-ladder remedies – multimodality (tempering), boundary (reparametrization), flat direction (Option A) – are not sampled (no overreach): a diagnosis pointing at them short-circuits to a certificate that names the proper remedy. The sub-Riemannian level requires an expected Fisher; without one the certificate prescribes supplying it – the showcase of an actionable prescription. The full orchestration over a compiled model is heavy:

# Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
res <- gdpar_geom_orchestrate(suite$G0_isotropic, n_grid = 1)
res$status     # "resolved" at euclidean_diagonal for the isotropic control

Option A: the flat direction (RG.6 part i, D96)

The forensic’s flat direction is the diana G7 (remedy reparam_eliminate), an out-of-ladder pathology. Option A removes it at the source: the per-slot additive scale sigma_a_k is compacted to the slots that actually carry free a coefficients (\(J_a^{\text{free}} > 0\)), via a sigma_a_idx map computed in transformed data. The criterion is \(J_a^{\text{free}} > 0\) (not “an a() was declared”): it captures both sub-cases of the flat direction – slots with no a() (the Tweedie’s phi ~ 1 / p ~ 1) and an intercept-only a() (a = ~ 1). For every model that carries an a() with a covariate in each slot, the index is the identity and the model is mathematically identical – the draws are bit-exact, which is how the K-path goldens stay green. The single golden with an intercept-only slot was re-bootstrapped and validated as a faithful marginalization (one parameter fewer, no surviving marginal moving beyond Monte Carlo error). Option A cleans the infinite \(\hat R\); it does not unstick the Tweedie – the geometry does (RG.7).

Integration: the bridge and the one-call fit (RG.6 part ii)

Two opt-in, net-new layers connect the orchestrator to the package fit engine, without touching gdpar():

The bridge return is plain data plus closures; you can inspect the engine target on a closure with no Stan:

geom_target <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * sum(t^2),
  grad_log_prob = function(t) -t, dim = 3L)
geom_target$dim
#> [1] 3
geom_target$grad_log_prob(c(1, 2, 3))
#> [1] -1 -2 -3

The cmdstan-backed bridge and the one-call fit are heavy (they compile and sample) and are reproduced by the gated script:

# Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
fit    <- gdpar(gdpar_bf(y ~ a(x), sigma ~ a(z)), data = d,
                family = gdpar_family("gaussian"), skip_id_check = TRUE)
bridge <- gdpar_geom_bridge(fit)
res    <- gdpar_geom_orchestrate(bridge$target, bridge$geom_target,
                                 reference = bridge$reference)

# Or, in one call, building the model through the shared .gdpar_K_build() seam:
res2 <- gdpar_geom_fit(gdpar_bf(y ~ a(x), sigma ~ a(z)), data = d,
                       family = gdpar_family("gaussian"), skip_id_check = TRUE)
res2$status

RG.7: applying the capability to the real count, and the certified limit

With the capability in place, RG.7 revisited the real Tweedie count of 9.2.O (four NE-USA sub-regions, unconstrained dimension d = 14). The case turned out to be much harder than the ladder anticipated, and the honest verdict is worth stating plainly because it sharpens one of this block’s own assumptions.

The canyon is genuinely non-Gaussian, and that is certified, not assumed. The count posterior is a stiff canyon: the observed-information condition number at the mode is 9e7–4e8. The sub-Riemannian level (the count’s intended remedy) and a constant-mass Laplace preconditioner can be made to accept on a single chain, but acceptance plus zero divergences on one chain do not prove mixing. The decisive test is a multi-chain gate (R-hat / ESS across independent chains) – a diagnostic added to the orchestrator after a careful external review, and now a permanent part of it. It caught a single-chain false positive: what looked “resolved” was not.

The decisive experiment was a whitened NUTS – Stan’s adaptive sampler with a dense metric equal to the Laplace precision (the -Hessian at the climbed mode), every chain initialised at the mode, 1000 warmup + 4 x 1000 sampling. If the posterior were locally Gaussian, whitening would make it isotropic and NUTS would mix trivially. It does not: R-hat = 1.117, ESS-bulk = 22.5, 5% divergences. Divergences in already-whitened coordinates can only mean the geometry is genuinely non-Gaussian – the canyon curves; a frozen Gaussian metric cannot track it, and a per-step Riemannian metric is infeasible here (minutes per exact Hessian). This revises the rigidity assumption used to motivate the sub-Riemannian level (Section Level 5): rigidity does not imply Gaussianity once the canyon curves – Bernstein–von Mises holds only in a vanishingly small neighbourhood of the mode. We keep the assumption where it was validated (RG.4, condition ~90) and record, honestly, where the real count breaks it.

The honest endpoint is the Laplace / plug-in predictive – competitor parity. When a posterior is a certified non-sampleable non-Gaussian canyon, the first-class outcome (charter section 2.4) is a certified, reproducible limit plus the Laplace / plug-in predictive. This is not a consolation: it is exactly the regime of the competitors on this coordinate. mgcv/REML reports a coefficient-Gaussian + dispersion/power plug-in; INLA reports a Laplace approximation; both are O(d/n)-exact Laplace / plug-in predictives, not exact MCMC. gdpar matches that regime: the mode plus the exact-Hessian Gaussian posterior N(mode, M^{-1}), pushed through the model’s own constrain_variables and predict() so the Tweedie density and change-of-variables Jacobian are scored identically to every competitor.

Crucially, the Laplace approximation is diagnosed, never trusted blindly. The extractor reports its fidelity against the true posterior over the same draws: the importance-sampling ESS, the PSIS Pareto-k of the weights p/q, and the mean log-density drop log p(mode) - log p(theta_s) against its Gaussian expectation d/2. On the real count these say loudly that the Laplace Gaussian is a wide, crude proposal for the curved canyon (Pareto-k 0.77–1.05, ESS < 5%, drop >> d/2) – which is the scientific finding, a measurement of how non-Gaussian the canyon is, not a defect to hide.

The result (full table in inst/benchmarks/results/block9_revalidation.md). gdpar’s mode (plug-in) Tweedie ELPD matches mgcv’s REML plug-in to < 0.4 ELPD units in all four cells – the same predictive surface, reached by a different optimiser. gdpar’s full Laplace average is 7–36 lower because it propagates the dispersion and power uncertainty that mgcv/REML plug in (a conservative, not a favourable, choice – gdpar gets no credit for plugging in); the gap is within ~1–1.5 standard errors, so the two are statistically indistinguishable. 9.2.O closes at 80/80. The geometry capability earned its keep here not by sampling the canyon but by diagnosing it (cond ~1e8, genuine non-Gaussianity), certifying the limit with a falsifiable multi-chain experiment, and delivering a competitor-parity predictive where naive NUTS produced nothing usable.

The automated fallback (gdpar_geom_laplace(), laplace_fallback)

Because this regime (a near-deterministic, genuinely non-Gaussian canyon) is now characterised, the Laplace machinery is promoted into R/ as a first-class, exported capability, gdpar_geom_laplace(), and wired into the orchestrator behind an opt-in flag. Given any engine target, it climbs to the mode (reading the same target and gradient the sampler uses), forms the precision M = -Hessian, and returns the mode + precision Gaussian N(mode, M^{-1}), optional draws, and the fidelity diagnostics above distilled into a single scalar label. On a clean Gaussian the label is "good":

A   <- matrix(c(2, 0.8, 0.8, 1), 2, 2)
mu  <- c(1, -0.5)
tgt <- gdpar_geom_target(
  log_prob      = function(t) -0.5 * as.numeric(t(t - mu) %*% A %*% (t - mu)),
  grad_log_prob = function(t) -as.numeric(A %*% (t - mu)),
  hessian       = function(t) -A, dim = 2L)
lap <- gdpar_geom_laplace(tgt, draws = 500L, seed = 1L)
c(label = lap$fit_quality_label, max_mode_err = max(abs(lap$mode - mu)))
#>        label max_mode_err 
#>       "good"          "0"

The orchestrator gains laplace_fallback = FALSE (the bit-identical opt-in default) and laplace_draws. When opted in and a run ends in a certified limit, it attaches the Laplace ($laplace) and relabels the status "certified_limit_laplace" – the sampling limit still stands and a labelled, fidelity-diagnosed Laplace posterior is provided, never advertised as exact:

# Reproduced by inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1)
res <- gdpar_geom_orchestrate(bridge$target, bridge$geom_target,
                              reference = bridge$reference,
                              laplace_fallback = TRUE, laplace_draws = 2000L)
res$status                 # "certified_limit_laplace"
res$laplace$fit_quality_label   # on the real 9.2.O canyon: "very_poor"

On the out-of-scope path (multimodality, a flat direction, a boundary) the Gaussian premise is violated, so the fallback is attached only when the curvature at the mode is genuinely positive-definite – otherwise the certificate is left untouched (no overreach). The full RG.7 result table is in inst/benchmarks/scripts/rg7_laplace_elpd.R and inst/benchmarks/results/block9_revalidation.md.

Honesty: demonstrated vs. conjectured, and no overreach

Following ORPHEUS §16.3: what is demonstrated here is correctness – every sampler level is Metropolis-exact with respect to the true posterior, so the geometry only ever changes efficiency, never the validity of the draws, and a mis-diagnosis wastes a bounded amount of budget, never corrupts the answer. What is conjectured / measured is efficiency: which level wins, the speed-ups, the prescriptions. gdpar integrates and orchestrates established methods; it does not claim to invent them. The anchors are Girolami–Calderhead 2011 (RMHMC, generalized implicit leapfrog), Betancourt (SoftAbs, E-BFMI, sampling pathologies), Neal (the funnel), Lu et al. 2017 (relativistic Monte Carlo), Livingstone–Faulkner–Roberts 2019 (kinetic energy under heavy tails), Montgomery 2002 (sub-Riemannian geometry, Chow–Rashevsky), Shahbaba–Lan–Johnson–Neal 2014 (split HMC), Roberts–Rosenthal 2007 (diminishing adaptation / containment), Li et al. 2018 (Hyperband / successive halving), Xu et al. 2008 (SATzilla, runtime-predicted algorithm selection), and, for the statistical-manifold frame, Rao–Amari information geometry, Pennec–Sommer–Fletcher and Bhattacharya–Patrangenaru (Riemannian statistics) as cited in the user’s manuscript.

Reproducing the heavy runs

Every eval = FALSE chunk above is reproduced end to end by the gated script

inst/scripts/geometry_pilots_deep.R

Run it with the gate set (it compiles cmdstan models and samples, so it takes several minutes and is never run during a normal package check):

Sys.setenv(GDPAR_RUN_GEOMETRY_PILOTS = "1")
source(system.file("scripts", "geometry_pilots_deep.R", package = "gdpar"))

It covers: the diagnostic on a compiled diana, each ladder level on a real cmdstan model, the full orchestrator (resolution and certified limit), and the bridge plus the one-call gdpar_geom_fit. This keeps the package build light while preserving full, falsifiable reproducibility of the geometry capability.