Skip to content

Part IV Function Reference 5

Jump to bottom
José Mauricio Gómez Julián edited this page Jun 29, 2026 · 5 revisions

← Part IV — Exhaustive Function Reference (4/7) · gdpar Wiki Home · Part IV — Exhaustive Function Reference (6/7) →

R/geometry_engine.R

gdpar_geom_target(object = NULL, log_prob = NULL, grad_log_prob = NULL, dim = NULL, hessian = NULL, param_names = NULL, data = NULL, simulate = NULL, score = NULL)

Purpose Constructs a target object that wraps a posterior's log-density and its gradient so the Block RG geometric sampling engine (gdpar_geom_hmc) can evaluate them on the unconstrained scale. Acts as a three-way adapter (decision A): density work is delegated to a compiled cmdstan model, an R closure, or a gdpar_geometry_suite() instance, while the integrator itself runs in pure R.

Arguments

Argument Type Meaning
object NULL, a list (suite instance), or a cmdstanr fit/model (CmdStanFit, CmdStanMCMC, CmdStanModel) Optional. If supplied the closure arguments are ignored. A suite instance must carry $log_prob, $grad_log_prob, and optionally $dim. A cmdstan object must have been compiled with compile_model_methods = TRUE exposing $grad_log_prob.
log_prob function(theta) or NULL Function of an unconstrained parameter vector returning the log-posterior density. Used only when object is NULL.
grad_log_prob function(theta) or NULL Function returning the gradient of log_prob at theta (unconstrained scale). Used only when object is NULL.
dim positive integer scalar or NULL Dimension of the unconstrained parameter vector. Required for the closure form and for a cmdstan object whose unconstrained dimension cannot always be inferred automatically.
hessian function(theta) or NULL Optional function returning the Hessian matrix of log_prob. Used by the Riemannian level (RG.3). For a cmdstan object the built-in $hessian method is used when available.
param_names character vector or NULL Optional names for the parameters. Falls back to paste0("theta[", seq_len(dim), "]").
data list or NULL Optional data list forwarded to a one-iteration sampling run when object is an uninstantiated CmdStanModel (needed only to expose $log_prob / $grad_log_prob).
simulate function(theta) or NULL Optional generative function returning one synthetic data set drawn from the model at theta. Supplying both simulate and score makes the target a generative target consumed by gdpar_geom_fisher_simulator.
score function(theta, y) or NULL Optional function returning the gradient of the model's log-likelihood of a data set y at theta on the unconstrained scale (the per-draw score whose outer product is the expected Fisher information).

Logic / Algorithm

  1. Input validation: simulate and score, if non-NULL, must be functions; otherwise gdpar_abort with class gdpar_input_error.
  2. Suite-instance path (object is a list with $log_prob and $grad_log_prob functions): Constructs the target via .gdpar_geom_target_obj, merging dim, param_names, simulate, score from object where the user-supplied arguments are NULL. Backend is "closure".
  3. cmdstan path (object is a CmdStanFit/CmdStanMCMC/CmdStanModel or has $grad_log_prob): Delegates to .gdpar_geom_target_cmdstan.
  4. Unrecognised object: gdpar_abort.
  5. Closure path (object is NULL): Both log_prob and grad_log_prob must be functions and dim must be supplied; dim is asserted to be a count. Constructs via .gdpar_geom_target_obj with backend "closure".

Returns A list of class c("gdpar_geom_target", "list") with elements:

Element Type Description
log_prob function(theta) Log-posterior on the unconstrained scale.
grad_log_prob function(theta) Gradient of log_prob.
hessian function(theta) or NULL Hessian of log_prob when available.
dim integer Dimension of the unconstrained parameter space.
param_names character vector Parameter names.
backend "closure" or "cmdstan" Which backend performs density evaluations.
simulate function(theta) or NULL Generative simulator if provided.
score function(theta, y) or NULL Score function if provided.

Notes

  • The function operates on the unconstrained scale (the scale on which Hamiltonian dynamics act and cmdstan's $grad_log_prob is defined).
  • When object is a CmdStanModel (not yet sampled), a throwaway one-chain, one-warmup-iteration, one-sampling-iteration fit is performed to expose the methods; the data are forwarded via data.
  • The %||% null-coalescing operator is used throughout to prefer user-supplied values over object-carried defaults.
  • Errors are raised via gdpar_abort with class "gdpar_input_error".

.gdpar_geom_target_obj(log_prob, grad_log_prob, hessian, dim, param_names, backend, simulate = NULL, score = NULL)

Purpose Internal constructor that assembles the raw list with class c("gdpar_geom_target", "list") from already-validated components. Called by gdpar_geom_target for both the closure and suite-instance paths and by .gdpar_geom_target_cmdstan for the cmdstan path.

Arguments

Argument Type Meaning
log_prob function(theta) Log-posterior (unconstrained scale).
grad_log_prob function(theta) Gradient of the log-posterior.
hessian function(theta) or NULL Hessian of the log-posterior (optional).
dim integer Unconstrained parameter dimension. Coerced to integer.
param_names character vector or NULL Names; defaults to paste0("theta[", seq_len(dim), "]") when NULL.
backend "closure" or "cmdstan" Label indicating the density backend.
simulate function(theta) or NULL Optional generative simulator.
score function(theta, y) or NULL Optional score function.

Returns A list of class c("gdpar_geom_target", "list") whose elements are the arguments stored directly (log_prob, grad_log_prob, hessian, dim, param_names, backend, simulate, score).

Notes

  • No input validation is performed; callers are responsible for passing valid arguments.
  • param_names falls back via %||% to sequential "theta[i]" labels.
  • dim is coerced with as.integer.

.gdpar_geom_target_cmdstan(object, dim, data, param_names, simulate = NULL, score = NULL)

Purpose Internal helper that wraps a cmdstanr fit or model so that its $log_prob, $grad_log_prob, and optionally $hessian methods are exposed as plain R closures taking an unconstrained parameter vector.

Arguments

Argument Type Meaning
object cmdstanr object (CmdStanModel, CmdStanFit, CmdStanMCMC) The compiled cmdstan object whose methods are to be wrapped.
dim integer or NULL Unconstrained parameter dimension. Required; the code aborts if NULL.
data list or NULL Data list forwarded to a one-iteration sampling run when object is a CmdStanModel.
param_names character vector or NULL Optional parameter names.
simulate function(theta) or NULL Optional generative simulator.
score function(theta, y) or NULL Optional score function.

Logic / Algorithm

  1. If object is a CmdStanModel, a minimal fit is obtained: object$sample(data = data, chains = 1, iter_warmup = 1, iter_sampling = 1, refresh = 0, show_messages = FALSE, show_exceptions = FALSE). This requires cmdstanr (suggested dependency check via require_suggested).
  2. If the resulting fit does not expose $grad_log_prob, an error is raised instructing the user to compile with compile_model_methods = TRUE.
  3. If dim is NULL, an error is raised.
  4. Three closures are built:
    • lp(theta) calls fit$log_prob(unconstrained_variables = theta).
    • gl(theta) calls fit$grad_log_prob(unconstrained_variables = theta), coerced with as.numeric.
    • he(theta) (when fit$hessian exists) calls fit$hessian(unconstrained_variables = theta)$hessian; otherwise NULL.
  5. Delegates to .gdpar_geom_target_obj with backend "cmdstan".

Returns A gdpar_geom_target list (same structure as .gdpar_geom_target_obj).

Notes

  • The throwaway fit is never used for inference; it exists solely to expose the model methods.
  • The as.numeric coercion on the gradient ensures a plain numeric vector regardless of cmdstanr's internal return type.
  • Errors raised: "gdpar_input_error" if $grad_log_prob is missing, or if dim is NULL.

print.gdpar_geom_target(x, ...)

Purpose S3 print method for gdpar_geom_target objects. Provides a concise summary of the target's backend, dimension, Hessian availability, and generative status.

Arguments

Argument Type Meaning
x gdpar_geom_target The target object to print.
... — Unused; present for S3 method compatibility.

Returns Invisibly returns x (standard print-method convention).

Side Effects Prints two lines to the console:

<gdpar_geom_target> backend: <backend> | dim: <dim>
  hessian available: <TRUE/FALSE> | generative: <TRUE/FALSE>

The generative flag is TRUE when both x$simulate and x$score are non-NULL.

Notes

  • No validation is performed on x.

gdpar_geom_metric_euclidean(dim = NULL, M = NULL)

Purpose Constructs the level-0/1 Euclidean (position-independent) metric for the Block RG geometry hierarchy. With the identity matrix (default) it is the standard diagonal Euclidean metric; with a user-supplied symmetric positive-definite matrix or diagonal variance vector it becomes a constant dense preconditioner — the remedy for anisotropic posteriors (straight canyons).

Arguments

Argument Type Meaning
dim positive integer or NULL Dimension. Required when M is NULL.
M symmetric positive-definite dim × dim matrix, length-dim positive numeric vector, or NULL Mass matrix. A vector is interpreted as diagonal entries and expanded to diag(M). Defaults to the dim × dim identity.

Logic / Algorithm

  1. If M is NULL, dim must be supplied and is asserted to be a count. M is set to diag(dim).
  2. If M is a vector (not a matrix), its entries are validated as finite and positive, then M <- diag(M).
  3. M is coerced to a matrix and checked for squareness.
  4. Cholesky decomposition ch <- chol(M) is attempted; failure raises an error (not symmetric positive-definite).
  5. Derived quantities are computed once:
    • Minv <- chol2inv(ch) — the inverse $M^{-1}$.
    • L <- t(ch) — the lower Cholesky factor $\mathbf{L}$ such that $M = \mathbf{L}\mathbf{L}^\top$; momenta are drawn as $\mathbf{p} = \mathbf{L}\mathbf{z}$ with $\mathbf{z} \sim \mathcal{N}(0, I)$.
    • ld <- 2 * sum(log(diag(ch))) — $\log\det M$.

Mathematics

$$M = \mathbf{L}\mathbf{L}^\top, \qquad \log\det M = 2\sum_{i=1}^d \log L_{ii}$$

$$M^{-1} = (\mathbf{L}\mathbf{L}^\top)^{-1} = \mathbf{L}^{-\top}\mathbf{L}^{-1}$$

Returns A list of class c("gdpar_geom_metric", "list") with:

Element Type Description
position_dependent FALSE Euclidean metric is constant.
dim integer Dimension $d$.
mass function(theta) → matrix Returns $M$ (ignores theta).
inv_mass function(theta) → matrix Returns $M^{-1}$ (ignores theta).
chol_mass function(theta) → matrix Returns lower Cholesky factor $\mathbf{L}$ (ignores theta).
logdet function(theta) → scalar Returns $\log\det M$ (ignores theta).

Notes

  • All four accessor functions accept theta for interface compatibility with the Riemannian metric (which is position-dependent) but ignore it entirely.
  • The identity metric (M = NULL, no dim with a vector) is the default for Euclidean HMC.
  • Errors: "gdpar_input_error" if neither dim nor M is supplied; if M is not square; if M has non-finite or non-positive diagonal entries (vector form); if M is not symmetric positive-definite (Cholesky fails).

gdpar_geom_metric_riemannian(target, curvature, fisher = NULL, dfisher = NULL, alpha = 1e6, floor = 1e-8, fd_step = 1e-5)

Purpose Constructs the level-3 Riemannian (position-dependent) metric of the Block RG geometry hierarchy: a mass matrix $M(\theta)$ that adapts to the local curvature of the log-posterior. This is the remedy for funnels and other variable-curvature geometries. Two curvature sources are supported: the expected Fisher information (model-specific, robust) and the SoftAbs regularisation of the observed Hessian (fully general, model-agnostic).

Arguments

Argument Type Meaning
target gdpar_geom_target or object accepted by gdpar_geom_target Provides the Hessian/gradient for SoftAbs and the parameter dimension.
curvature "fisher" or "softabs" Curvature source selector.
fisher function(theta) → matrix or NULL For curvature = "fisher": returns the $d \times d$ expected Fisher information matrix $I(\theta)$ at theta.
dfisher function(theta) → list or NULL For curvature = "fisher": returns a length-$d$ list whose $k$-th element is $\partial I(\theta)/\partial\theta_k$. When NULL, derivatives are finite-differenced from fisher.
alpha positive scalar (default 1e6) SoftAbs softening parameter $\alpha$. Larger values track the true eigenvalues more faithfully ($\lambda\coth(\alpha\lambda) \to
floor positive scalar (default 1e-8) Minimum eigenvalue imposed on $M(\theta)$ to ensure strict positive-definiteness at flat/non-identified directions.
fd_step positive scalar (default 1e-5) Finite-difference step size for: (a) the Hessian in the SoftAbs path when target$hessian is NULL; (b) the metric derivative $\partial M/\partial\theta_k$ when an analytic derivative is unavailable.

Mathematics

SoftAbs regularisation. Given the Hessian $H = -\nabla^2 \log\pi(\theta)$ with eigendecomposition $H = Q,\mathrm{diag}(\lambda_1,\dots,\lambda_d),Q^\top$, the SoftAbs metric is:

$$M(\theta) = Q,\mathrm{diag}!\bigl(\lambda_i\coth(\alpha,\lambda_i)\bigr),Q^\top$$

where $\coth(x) = \frac{e^x + e^{-x}}{e^x - e^{-x}}$. Near zero, $\lambda\coth(\alpha\lambda) \approx \frac{1}{\alpha}$, so flat directions are floored at $1/\alpha$. An additional hard floor of floor is applied to the eigenvalues.

Fisher metric. $M(\theta) = I(\theta)$, the expected Fisher information. Derivatives $\partial M/\partial\theta_k$ are supplied analytically via dfisher or computed by central finite differences:

$$\frac{\partial M}{\partial\theta_k} \approx \frac{M(\theta + \varepsilon, e_k) - M(\theta - \varepsilon, e_k)}{2\varepsilon}$$

Returns A list of class c("gdpar_geom_metric", "list") with:

Element Type Description
position_dependent TRUE Riemannian metric varies with $\theta$.
dim integer Parameter dimension $d$.
mass function(theta) → matrix $M(\theta)$.
inv_mass function(theta) → matrix $M(\theta)^{-1}$.
chol_mass function(theta) → matrix Lower Cholesky factor of $M(\theta)$.
logdet function(theta) → scalar $\log\det M(\theta)$.
dmass function(theta) → list Length-$d$ list of $\partial M(\theta)/\partial\theta_k$, consumed by the generalised implicit leapfrog.

Notes

  • The metric is a preconditioner, not part of the target density; the Metropolis correction with the exact log-density ensures correctness regardless of metric quality.
  • The Riemannian interface is a strict superset of the Euclidean interface (it adds dmass), so gdpar_geom_hmc can dispatch uniformly.
  • SoftAbs requires the Hessian: it uses target$hessian if available, otherwise finite-differences it from target$grad_log_prob with step fd_step.
  • The fisher source is the primary, maximally robust choice; its learned amortisation across model families is a separate sub-phase.
  • Errors follow gdpar_abort with class "gdpar_input_error".

gdpar_geom_metric_riemannian(target, curvature = c("fisher", "softabs"), fisher = NULL, dfisher = NULL, alpha = 1e6, floor = 1e-8, fd_step = 1e-4)

Purpose
Constructs a position-dependent Riemannian metric (mass matrix) for a target distribution using either the exact expected Fisher information matrix (if available and closed-form) or the SoftAbs transform of the Hessian of the negative log-probability. This is the primary factory for Capa 1 metrics in the Block RG hierarchy.

Arguments

  • target: A gdpar_geom_target object or an object coercible to it via gdpar_geom_target(). Must provide at least dim, grad_log_prob, and optionally hessian.
  • curvature: Character string, either "fisher" or "softabs", selecting the curvature estimator.
  • fisher: A function theta -> matrix returning the expected Fisher information matrix at theta. Required when curvature = "fisher".
  • dfisher: Optional function theta -> list of matrices returning the partial derivatives of the Fisher matrix with respect to each component of theta. If NULL and curvature = "fisher", derivatives are computed via finite differences.
  • alpha: Positive scalar controlling the SoftAbs softening parameter (larger values yield a sharper approximation of the absolute eigenvalue). Used only for curvature = "softabs".
  • floor: Non-negative scalar; eigenvalues of the mass matrix are floored to this value to ensure positive definiteness.
  • fd_step: Positive scalar step size for central finite-difference approximations of gradients and Hessians.

Mathematics

  • Fisher case: $M(\theta) = \text{floor\_spd}(\text{fisher}(\theta))$, where $\text{floor\_spd}$ projects onto the SPD cone.
  • SoftAbs case: Let $H_U(\theta)$ be the Hessian of $-\log p(\theta)$. The mass matrix is $M(\theta) = Q \cdot \text{diag}\bigl(\lambda_i \coth(\alpha \lambda_i)\bigr) \cdot Q^\top$, where $H_U = Q \Lambda Q^\top$ is the eigendecomposition, with eigenvalues $\lambda_i$ floored to floor. Derivatives use the Daleckii–Krein formula.

Returns
A list of class gdpar_geom_metric (with position_dependent = TRUE) containing functions:

  • mass(theta): the mass matrix $M(\theta)$.
  • inv_mass(theta): $M(\theta)^{-1}$.
  • chol_mass(theta): the upper Cholesky factor of $M(\theta)$.
  • logdet(theta): $\log\det M(\theta)$.
  • dmass(theta): a list of matrices $\partial M / \partial \theta_k$.

Also includes metadata fields: dim, metric_kind ("fisher" or "softabs"), alpha.

Notes

  • Raises an error of class "gdpar_input_error" if curvature = "fisher" and fisher is not a function.
  • For the SoftAbs case, if the target does not expose a hessian method, the Hessian is approximated by central finite differences of the gradient (symmetrized).
  • The returned object uses internal caching to reuse the Cholesky factorization across calls for the same theta (optimization for the implicit leapfrog integrator).

.gdpar_geom_metric_from_mass(mass_fn, dmass_fn, d, kind, alpha)

Purpose
Internal helper that assembles a gdpar_geom_metric object from supplied mass and derivative functions, incorporating a caching mechanism for the mass matrix, its Cholesky factor, and inverse at the last queried position. This avoids redundant factorizations during fixed-point iterations of the implicit leapfrog.

Arguments

  • mass_fn: A function theta -> matrix returning the mass matrix $M(\theta)$.
  • dmass_fn: A function theta -> list of matrices returning $\partial M(\theta)/\partial \theta_k$.
  • d: Integer, the dimension of the parameter space.
  • kind: Character string indicating the metric type (e.g., "fisher", "softabs", "gp_fisher").
  • alpha: Numeric scalar; the SoftAbs smoothing parameter (used for metadata).

Returns
A list of class gdpar_geom_metric with the standard interface (mass, inv_mass, chol_mass, logdet, dmass), plus fields position_dependent = TRUE, dim, metric_kind, alpha. The functions mass, inv_mass, chol_mass, and logdet are backed by a cache that stores the last evaluated theta, its mass matrix, Cholesky factor, inverse, and log-determinant.

Notes

  • The cache is implemented via an environment (new.env(parent = emptyenv())).
  • The ensure function updates the cache only if the input theta differs from the cached one (using identical).

.gdpar_geom_chol_spd(M)

Purpose
Computes the lower Cholesky factor of a symmetric matrix, with symmetrization and fallback jitter to handle numerically near-singular but valid metric matrices.

Arguments

  • M: A square numeric matrix (assumed symmetric or nearly symmetric).

Mathematics

  1. Symmetrize: $M \leftarrow (M + M^\top)/2$.
  2. Attempt Cholesky decomposition $M = L L^\top$.
  3. If it fails, add a jitter $jI$ to the diagonal, where $j = 10^{-8} \cdot \max(|\mathrm{diag}(M)|, 1)$, and recompute.

Returns
The lower Cholesky factor $L$ (a lower-triangular matrix).

Notes

  • The jitter ensures the matrix becomes strictly positive definite, preserving SPD structure for metrics that are theoretically SPD but numerically borderline.

.gdpar_geom_floor_spd(M, floor)

Purpose
Projects a symmetric matrix onto the symmetric positive-definite (SPD) cone by flooring its eigenvalues.

Arguments

  • M: A square numeric matrix (symmetric).
  • floor: Non-negative scalar; the minimum allowed eigenvalue.

Mathematics

  1. Symmetrize: $M \leftarrow (M + M^\top)/2$.
  2. Eigendecomposition: $M = Q \Lambda Q^\top$.
  3. Floor eigenvalues: $\Lambda' = \max(\Lambda, \text{floor})$.
  4. Reconstruct: $M' = Q \Lambda' Q^\top$.

Returns
The projected matrix $M'$, which is symmetric and has all eigenvalues at least floor.

Notes

  • This operation is idempotent if M is already SPD with eigenvalues ≥ floor.

.gdpar_geom_softabs_vals(lambda, alpha)

Purpose
Evaluates the SoftAbs transform $f(\lambda) = \lambda \coth(\alpha \lambda)$ for a vector of eigenvalues, with a stable series expansion near zero.

Arguments

  • lambda: Numeric vector of eigenvalues.
  • alpha: Positive scalar.

Mathematics
Let $u = \alpha \lambda$. Then
$$ f(\lambda) = \frac{1}{\alpha} \times \begin{cases} \frac{u}{\tanh u} & \text{if } |u| \ge 10^{-4}, \ 1 + \frac{u^2}{3} & \text{if } |u| < 10^{-4} \end{cases} $$
The near-zero expansion uses the series of $u \coth u \approx 1 + u^2/3 + O(u^4)$.

Returns
A numeric vector of the same length as lambda containing $f(\lambda)$.

Notes

  • The threshold $10^{-4}$ is chosen to avoid loss of precision in the direct formula.

.gdpar_geom_softabs_deriv(lambda, alpha)

Purpose
Computes the derivative of the SoftAbs transform $f'(\lambda) = \frac{d}{d\lambda} [\lambda \coth(\alpha \lambda)]$.

Arguments

  • lambda: Numeric vector of eigenvalues.
  • alpha: Positive scalar.

Mathematics
Let $u = \alpha \lambda$. Then
$$ f'(\lambda) = \frac{1}{\alpha} \times \begin{cases} \frac{1}{\tanh u} - \frac{u}{\sinh^2 u} & \text{if } |u| \ge 10^{-3}, \ \frac{2}{3} u & \text{if } |u| < 10^{-3} \end{cases} $$
The near-zero expansion uses the Taylor series $f'(\lambda) \approx \frac{2}{3} \alpha \lambda + O(\lambda^3)$.

Returns
A numeric vector of derivatives, same length as lambda.

Notes

  • The threshold $10^{-3}$ balances accuracy and stability; the linear approximation is valid for small $u$.

.gdpar_geom_softabs_mass(HU, alpha, floor)

Purpose
Constructs the SoftAbs mass matrix from the Hessian of the negative log-probability.

Arguments

  • HU: Square matrix, the Hessian of $-\log p(\theta)$ (symmetric).
  • alpha: Positive scalar for SoftAbs.
  • floor: Non-negative scalar for eigenvalue flooring.

Mathematics

  1. Symmetrize $HU$.
  2. Eigendecomposition: $HU = Q \Lambda Q^\top$.
  3. Transform eigenvalues: $\tilde{\lambda} = \max(f(\lambda), \text{floor})$, where $f$ is the SoftAbs function.
  4. Reconstruct: $M = Q \mathrm{diag}(\tilde{\lambda}) Q^\top$.

Returns
The SPD mass matrix $M$.

Notes

  • The flooring is applied after the SoftAbs transform to ensure positive definiteness even for zero or negative eigenvalues of $HU$.

.gdpar_geom_softabs_dmass(theta, hess_lp_fn, alpha, h)

Purpose
Computes the derivative of the SoftAbs mass matrix with respect to each component of $\theta$ using the Daleckii–Krein formula and central finite differences for the third-order derivative tensor.

Arguments

  • theta: Numeric vector, current parameter position.
  • hess_lp_fn: Function returning the Hessian of the log-probability (i.e., $-\text{Hessian of } \log p$ is hess_lp_fn(theta)).
  • alpha: Positive scalar for SoftAbs.
  • h: Positive scalar finite-difference step.

Mathematics

  1. Compute $HU = -\text{hess\_lp\_fn}(\theta)$ and its eigendecomposition $HU = Q \Lambda Q^\top$.
  2. Compute the matrix of divided differences $R$:
    $R_{ij} = \begin{cases} \frac{f(\lambda_i) - f(\lambda_j)}{\lambda_i - \lambda_j} & \text{if } |\lambda_i - \lambda_j| > 10^{-8}, \ f'((\lambda_i + \lambda_j)/2) & \text{otherwise}. \end{cases}$
  3. For each component $k = 1,\dots,d$:
    • Finite-difference the Hessian: $dHU_k = \frac{HU(\theta + h e_k) - HU(\theta - h e_k)}{4h} + \text{symmetric part}$? Wait, the code does: 
      HUp <- -hess_lp_fn(theta + e); HUm <- -hess_lp_fn(theta - e)
dHU <- ((HUp + t(HUp)) - (HUm + t(HUm))) / (4 * h)
      This is a central difference of the symmetrized Hessian.
    • Compute $A = Q^\top dHU_k Q$.
    • Then $\frac{\partial M}{\partial \theta_k} = Q (R \odot A) Q^\top$, where $\odot$ is the Hadamard (element-wise) product.

Returns
A list of length $d$ (the dimension of $\theta$), where each element is a square matrix $\frac{\partial M}{\partial \theta_k}$.

Notes

  • The divided-difference matrix $R$ accounts for repeated eigenvalues; the derivative at a repeated eigenvalue is taken as $f'(\lambda)$ (the derivative, which is well-defined).
  • The finite-difference step h is applied symmetrically to reduce bias.

.gdpar_geom_hessian_fn(target, h)

Purpose
Returns a function that computes the Hessian of the log-probability, either analytically (if the target exposes a hessian method) or via central finite differences of the gradient.

Arguments

  • target: A gdpar_geom_target object (or similar) that must have a grad_log_prob method and optionally a hessian method.
  • h: Positive scalar finite-difference step.

Mathematics
If no analytic Hessian is available, the Hessian $H$ is approximated columnwise:
$$ H_{:,k} \approx \frac{\nabla \log p(\theta + h e_k) - \nabla \log p(\theta - h e_k)}{2h}, $$
and then symmetrized: $H \leftarrow (H + H^\top)/2$.

Returns
A function theta -> matrix returning the Hessian of the log-probability at theta.

Notes

  • The returned function is used by the SoftAbs machinery to obtain $HU = -\text{Hessian}$.

.gdpar_geom_fd_dmass(mass_fn, theta, h)

Purpose
Approximates the derivative of a matrix-valued mass function by central finite differences.

Arguments

  • mass_fn: A function theta -> matrix returning the mass matrix.
  • theta: Numeric vector, current position.
  • h: Positive scalar finite-difference step.

Mathematics
For each component $k$:
$$ \frac{\partial M}{\partial \theta_k} \approx \frac{\text{mass_fn}(\theta + h e_k) - \text{mass_fn}(\theta - h e_k)}{2h}. $$

Returns
A list of matrices, one per component of theta, representing the partial derivatives of the mass matrix.

Notes

  • Used when an analytic derivative is not provided (e.g., for the Fisher case when dfisher is NULL).

.gdpar_geom_tri_index(d)

Purpose
Generates the row and column indices of the lower-triangular part of a $d \times d$ matrix (including the diagonal) in column-major order. This is used to map between a symmetric matrix and a vector of unique entries.

Arguments

  • d: Integer, the dimension of the matrix.

Returns
A matrix with two columns and $d(d+1)/2$ rows. Each row gives the (row, col) index of an element in the lower triangle (including diagonal), ordered column by column (i.e., all indices for column 1, then column 2, etc.).

Notes

  • The order is consistent with R's lower.tri(..., diag = TRUE) and the column-major storage order.

.gdpar_geom_logchol(M, idx)

Purpose
Computes the log-Cholesky parametrisation of an SPD matrix $M$. The lower Cholesky factor $L$ is represented by a vector $\psi$ where the diagonal entries are $\log(L_{ii})$ and the off-diagonal entries are $L_{ij}$ (for $i &gt; j$), in column-major lower-triangular order.

Arguments

  • M: A symmetric positive-definite matrix.
  • idx: A matrix of lower-triangular indices as returned by .gdpar_geom_tri_index(d).

Mathematics

  1. Symmetrize $M$ and compute its lower Cholesky factor: $M = L L^\top$.
  2. Define $\psi$ such that for each index $(i,j)$ in idx:
    • If $i = j$ (diagonal), $\psi_{\text{idx}} = \log(L_{ii})$.
    • If $i &gt; j$ (off-diagonal), $\psi_{\text{idx}} = L_{ij}$.

Returns
A list with components:

  • L: The lower Cholesky factor (a lower-triangular matrix).
  • psi: Numeric vector of length $d(d+1)/2$ containing the log-Cholesky parameters.

Notes

  • This parametrisation ensures that any $\psi \in \mathbb{R}^{d(d+1)/2}$ maps to a unique SPD matrix $M = L L^\top$ with $L_{ii} = \exp(\psi_{ii}) &gt; 0$, making the map smooth and surjective onto the SPD cone.

.gdpar_geom_L_from_psi(psi, idx, d)

Purpose
Reconstructs the lower Cholesky factor $L$ from the log-Cholesky parameters $\psi$.

Arguments

  • psi: Numeric vector of length $d(d+1)/2$, the log-Cholesky parameters.
  • idx: Matrix of lower-triangular indices (from .gdpar_geom_tri_index(d)).
  • d: Integer, the dimension.

Mathematics
For each index $(i,j)$ in idx:

  • If $i = j$ (diagonal), $L_{ij} = \exp(\psi_{\text{idx}})$.
  • If $i &gt; j$ (off-diagonal), $L_{ij} = \psi_{\text{idx}}$.
    All other entries of $L$ are zero.

Returns
A lower-triangular matrix $L$ of dimension $d \times d$.

.gdpar_geom_dM_from_dpsi(L, dpsi, idx, d)

Purpose
Computes the derivative of the mass matrix $M = L L^\top$ with respect to the log-Cholesky parameters $\psi$, given the current $L$ and a vector of derivatives $d\psi$.

Arguments

  • L: The current lower Cholesky factor (lower-triangular matrix).
  • dpsi: Numeric vector of the same length as psi, representing $d\psi$.
  • idx: Matrix of lower-triangular indices.
  • d: Integer, the dimension.

Mathematics
First, construct $dL$ from $d\psi$:

  • For diagonal entries ($i = j$): $dL_{ii} = L_{ii} \cdot d\psi_{ii}$ (since $\psi_{ii} = \log L_{ii}$, so $dL_{ii}/L_{ii} = d\psi_{ii}$).
  • For off-diagonal entries ($i &gt; j$): $dL_{ij} = d\psi_{ij}$.

Then, the derivative of $M$ is:
$$ dM = dL \cdot L^\top + L \cdot dL^\top. $$

Returns
A symmetric matrix $dM$ of dimension $d \times d$.

.gdpar_geom_dpsi_from_dM(L, dM, idx, d)

Purpose
Computes the derivative of the log-Cholesky parameters $\psi$ with respect to the mass matrix $M$, given the current $L$ and a matrix derivative $dM$. This is the inverse map of .gdpar_geom_dM_from_dpsi.

Arguments

  • L: The current lower Cholesky factor (lower-triangular matrix).
  • dM: A symmetric matrix representing $dM$.
  • idx: Matrix of lower-triangular indices.
  • d: Integer, the dimension.

Mathematics

  1. Compute $P = L^{-1} dM L^{-\top}$ (symmetric).
  2. Define $\Phi$ as the lower-triangular part of $P$ with the diagonal halved:
    • $\Phi_{ij} = P_{ij}$ for $i &gt; j$,
    • $\Phi_{ii} = P_{ii}/2$,
    • $\Phi_{ij} = 0$ for $i &lt; j$.
  3. Compute $dL = L \Phi$.
  4. Extract $d\psi$ from $dL$:
    • For diagonal entries: $d\psi_{ii} = dL_{ii} / L_{ii}$.
    • For off-diagonal entries: $d\psi_{ij} = dL_{ij}$.

Returns
A numeric vector $d\psi$ of length $d(d+1)/2$.

Notes

  • This function is the differential of the inverse map from $M$ to $\psi$, useful for backpropagation or sensitivity analysis in the log-Cholesky space.

gdpar_geom_metric_gp_fisher(target, fisher, sites, weights = NULL, lengthscale = NULL, nugget = 1e-6, alpha = 1e6, floor = 1e-8, fd_step = 1e-4)

Purpose
Constructs a learned Gaussian-process (GP) Riemannian metric that serves as a surrogate for the expected Fisher information matrix. The GP models the residual between the SoftAbs curvature (the always-available Capa 1 metric) and the true Fisher information, in log-Cholesky parametrisation. This is the Capa 2 metric for the general case where the Fisher has no closed form.

Arguments

  • target: A gdpar_geom_target object (or coercible) providing the SoftAbs mean function (via its Hessian/gradient) and the dimension.
  • fisher: A function theta -> matrix that returns the expected Fisher information matrix at theta (the ground truth for training).
  • sites: An $m \times d$ matrix of reservoir positions (one per row) where the Fisher information is evaluated to train the GP.
  • weights: Optional positive numeric vector of length $m$; importance weights for the sites (e.g., to correct for biased reservoir sampling). Defaults to equal weights.
  • lengthscale: Optional positive scalar; the RBF kernel lengthscale on standardised inputs. If NULL, defaults to the median pairwise Euclidean distance among the standardised sites.
  • nugget: Non-negative scalar; the GP observation noise variance. Small values interpolate the Fisher at sites; larger values smooth it.
  • alpha, floor, fd_step: Parameters passed to the SoftAbs machinery for the mean function (same as in gdpar_geom_metric_riemannian).

Mathematics
The surrogate models the log-Cholesky parameters $\psi$ of the Fisher information as:
$$ \psi(\theta) = \psi_{\text{SoftAbs}}(\theta) + \delta(\theta), $$
where $\psi_{\text{SoftAbs}}$ are the log-Cholesky parameters of the SoftAbs mass matrix, and $\delta(\theta)$ is a GP with zero mean and a radial basis function kernel:
$$ k(\theta, \theta') = \sigma^2 \exp\left( -\frac{|\theta - \theta'|^2}{2 \ell^2} \right) + \text{nugget} \cdot \delta_{\theta,\theta'}. $$
The GP is trained on the residuals $\delta_i = \psi_{\text{Fisher}}(\theta_i) - \psi_{\text{SoftAbs}}(\theta_i)$ at the reservoir sites $\theta_i$. Predictions at a new $\theta$ yield the posterior mean $\psi_{\text{GP}}(\theta)$, from which the mass matrix is reconstructed via $M(\theta) = L(\theta) L(\theta)^\top$ with $L$ from $\psi_{\text{GP}}(\theta)$.

Returns
A list of class gdpar_geom_metric (with position_dependent = TRUE, metric_kind = "gp_fisher") containing the standard metric interface:

  • mass(theta), inv_mass(theta), chol_mass(theta), logdet(theta), dmass(theta).
    Additionally:
  • novelty(theta): A function returning the predictive standard deviation of $\psi$ at theta (a measure of epistemic uncertainty; higher values indicate out-of-distribution regions).
  • n_sites: Integer, the number of training sites.
  • lengthscale: The lengthscale used (may be the computed median if not provided).

Notes

  • The metric degrades gracefully to the SoftAbs mean when far from the reservoir because the kernel decays to zero.
  • The sampler remains exact because the metric is only a preconditioner; the Metropolis correction uses the exact target density.
  • The derivative dmass(theta) is computed analytically: the GP derivative plus the Daleckii–Krein derivative of the SoftAbs mean, all pushed through the log-Cholesky map.
  • The GP is implemented using exact (non-sparse) inference; for large reservoirs, this may be computationally expensive.
  • The function performs input standardisation (zero mean, unit variance) internally before applying the kernel.

gdpar_geom_metric_gp_fisher(target, fisher, sites, weights = NULL, lengthscale = NULL, nugget = 1e-6, alpha = 1e6, floor = 1e-8, fd_step = 1e-4)

Purpose Constructs a learned, position-dependent Riemannian metric whose mass matrix is a Gaussian-process surrogate of the expected Fisher information (in log-Cholesky coordinates). The GP's mean function is the SoftAbs approximation to the negative log-posterior Hessian (the deterministic "cold-start" curvature), and the process learns the smooth residual between the true expected Fisher and that SoftAbs mean evaluated at a set of reservoir sites. This implements the highest fidelity level of the Block RG hierarchy, where the Fisher information is approximated point-wise via an RBF interpolant over the log-Cholesky parameterisation of SPD matrices.

Arguments

  • target — A gdpar_geom_target object (or any object accepted by its constructor). Provides log_prob, grad_log_prob, and dim. Converted internally if needed.
  • fisher — A function fisher(theta) returning the dim × dim expected Fisher information matrix at position theta. Must be supplied by the caller (e.g. the output of gdpar_geom_fisher_simulator or a closed-form evaluator).
  • sites — A numeric matrix with d columns, each row a reservoir site where the Fisher and SoftAbs are evaluated to train the GP. If a bare numeric vector of length d is passed it is reshaped into a single-row matrix. Must have at least two rows.
  • weights — Optional positive numeric vector of length m = nrow(sites). Importance weights for the reservoir sites; when non-NULL, the nugget noise at each site is scaled as nugget / (w / mean(w)), implementing importance-weighted GP regression (inverse-variance weighting where higher-weight sites receive less noise).
  • lengthscale — Positive scalar length-scale of the RBF kernel. If NULL, set automatically by the median-distance heuristic: compute pairwise Euclidean distances in the standardised site matrix, take the median of strictly positive distances; if that fails, default to sqrt(d).
  • nugget — Non-negative scalar; the base diagonal noise added to the kernel matrix. Ensures numerical positive-definiteness of K + noise.
  • alpha — Non-negative scalar SoftAbs hyperparameter controlling the smoothness of the absolute-value approximation to eigenvalues (larger = sharper). Passed to .gdpar_geom_softabs_mass.
  • floor — Non-negative scalar; minimum eigenvalue floor applied inside the SoftAbs mass computation.
  • fd_step — Positive scalar finite-difference step for the Hessian approximation and for the derivative of the SoftAbs mass matrix. Passed to .gdpar_geom_hessian_fn and .gdpar_geom_softabs_dmass.

Mathematics

The function maps each site $\mathbf{s}_r$ and the query point $\theta$ to a standardised coordinate $\mathbf{z} = (\theta - \mathbf{c})/\mathbf{s}$ where $\mathbf{c}$ and $\mathbf{s}$ are the column means and standard deviations of the site matrix.

The RBF (squared-exponential) kernel is

$$K(\mathbf{s}_i, \mathbf{s}_j) = \exp!\left(-\frac{|\mathbf{z}_i - \mathbf{z}_j|^2}{2\ell^2}\right)$$

with $\ell$ the length-scale. The kernel matrix with importance noise is

$$\mathbf{K}_n = \mathbf{K} + \mathrm{diag}(\boldsymbol{\nu}), \qquad \nu_r = \begin{cases}\texttt{nugget} & \text{if weights is NULL},\ \texttt{nugget},\big/,(w_r / \bar{w}) & \text{otherwise}.\end{cases}$$

At every site $r$, two quantities are computed and mapped through the log-Cholesky parametrisation $\psi = \text{logchol}(M)$ of SPD matrices (indexed by the upper-triangular index set $\mathrm{idx}$ with $q$ entries):

$$\boldsymbol{\psi}^{(F)}_r = \text{logchol}\big(I(\mathbf{s}_r)\big), \qquad \boldsymbol{\psi}^{(0)}_r = \text{logchol}\big(M_{\text{SoftAbs}}(\mathbf{s}_r)\big).$$

The GP target is the residual $\mathbf{Y}_{\mathrm{res},r} = \boldsymbol{\psi}^{(F)}_r - \boldsymbol{\psi}^{(0)}_r$. The dual coefficients are

$$\mathbf{C} = \mathbf{K}_n^{-1},\mathbf{Y}_{\mathrm{res}} \in \mathbb{R}^{m \times q}.$$

At a new position $\theta$, the predictive mean in log-Cholesky space is

$$\boldsymbol{\psi}(\theta) = \boldsymbol{\psi}^{(0)}(\theta) + \mathbf{k}_\star^\top \mathbf{C}$$

where $\mathbf{k}_\star$ is the kernel vector between $\theta$ and the sites. The mass matrix is recovered via the inverse log-Cholesky map: $M(\theta) = L(\theta),L(\theta)^\top$ with $L = \texttt{L\_from\_psi}(\boldsymbol{\psi})$.

For derivatives, the GP correction to the log-Cholesky differential is

$$\frac{\partial \boldsymbol{\psi}}{\partial \theta_j} = \frac{\partial \boldsymbol{\psi}^{(0)}}{\partial \theta_j} + \left(\frac{\partial \mathbf{k}_\star}{\partial \theta_j}\right)^\top \mathbf{C}$$

where $\frac{\partial k_{\star,r}}{\partial \theta_j} = k_{\star,r},\frac{z_j - z_{r,j}}{\ell^2 s_j}$.

The novelty (predictive standard deviation at a test point, without the noise term) is

$$\text{novelty}(\theta) = \sqrt{\max!\left(0,; 1 - \mathbf{k}_\star^\top \mathbf{K}_n^{-1},\mathbf{k}_\star\right)}.$$

Returns A list of class c("gdpar_geom_metric", "list") with position_dependent = TRUE and metric_kind = "gp_fisher". Contains closures mass(theta), inv_mass(theta), chol_mass(theta), logdet(theta), dmass(theta) (returning a length-d list of SPD matrices), and novelty(theta). Also exposes dim, alpha, n_sites, lengthscale.

Notes

  • All mutable per-position state is cached in an internal environment st; repeated calls at the same theta skip recomputation (guarded by identical).
  • The SoftAbs mean serves as a structural control variate: it is deterministic and carries the bulk of the curvature; the GP learns only the smooth residual, so few reservoir sites suffice.
  • dmass computes the derivative of the SoftAbs mean analytically via .gdpar_geom_softabs_dmass (cached on first derivative call) and adds the GP correction analytically.
  • Input validation raises errors of class "gdpar_input_error" via gdpar_abort for malformed inputs.
  • The log-Cholesky index set .gdpar_geom_tri_index(d) produces $q$ upper-triangular entries of the $d \times d$ factor.

gdpar_geom_reservoir(target, n_sites = 50L, metric = NULL, epsilon = 0.25, L = 15L, n_warmup = 200L, init = NULL, seed = NULL)

Purpose Runs a short HMC warmup to collect a reservoir of positions that are later used as the sites argument of gdpar_geom_metric_gp_fisher. This is the first phase of the two-phase, decoupled-archivist design (ORPHEUS-PIMC-A section 16): a cheap Euclidean (or SoftAbs) warmup provides geometrically representative sites at which the expensive expected Fisher is evaluated.

Arguments

  • target — A gdpar_geom_target (or any object accepted by its constructor). Converted internally.
  • n_sites — Integer: the number of retained post-warmup draws to return (the reservoir size). Must be a positive integer (validated by assert_count).
  • metric — Optional metric object for the warmup HMC run; defaults to the Euclidean identity inside gdpar_geom_hmc. A gdpar_geom_metric_riemannian SoftAbs metric is recommended for curved targets.
  • epsilon — Positive scalar: leapfrog step size for the warmup run.
  • L — Positive integer: trajectory length (number of leapfrog steps) for the warmup run.
  • n_warmup — Non-negative integer: number of discarded warmup iterations before retaining draws.
  • init — Optional numeric vector of length dim giving the initial position; defaults to zeros inside gdpar_geom_hmc.
  • seed — Optional integer seed for reproducibility of the HMC run.

Mathematics None beyond the underlying HMC sampler.

Returns A numeric matrix of dimension n_sites × dim, the retained reservoir positions.

Notes

  • The function is a thin wrapper: it calls gdpar_geom_hmc with n_iter = n_sites and returns fit$draws.
  • The reservoir design decouples the (potentially expensive) Fisher evaluation from the sampling run, allowing the surrogate to be trained offline.

.gdpar_geom_seed_from_theta(theta)

Purpose Derives a deterministic, non-negative integer key from a position vector so that simulation-based evaluations at $\theta$ are reproducible regardless of the order of calls. Used by gdpar_geom_fisher_simulator to reseed the RNG before each expected-Fisher evaluation.

Arguments

  • theta — A numeric vector (the position). Non-finite entries are replaced with zero.

Mathematics

$$\text{key} = \left(\sum_{i=1}^{d} \lfloor|\theta_i|\cdot 1009 + 1\rfloor \cdot (7919,i)\right) \bmod 2147483646 + 1$$

The primes 1009 and 7919 are chosen to decorrelate nearby positions while keeping the computation cheap. The modulus 2147483646 = $2^{31} - 2$ ensures the result fits in a positive 32-bit signed integer (1 through $2^{31}-1$).

Returns A single integer in ${1, 2, \ldots, 2147483647}$, suitable for passing to set.seed.

Notes

  • Non-finite entries of theta are silently zeroed before the computation.
  • The key is designed to be a fast hash, not a cryptographic hash; collisions are possible but benign (they merely produce correlated simulations, which the GP surrogate absorbs).

gdpar_geom_fisher_simulator(target, n_sim = 64L, seed = 1L, floor = 1e-8)

Purpose Constructs a deterministic, simulation-based estimator of the expected Fisher information matrix for generative models where no closed-form expression exists. The returned function is suitable as the fisher argument of gdpar_geom_metric_gp_fisher or gdpar_geom_metric_subriemannian.

Arguments

  • target — A gdpar_geom_target (or any object accepted by its constructor) that must carry both simulate(theta) (draws one synthetic data set from the model at $\theta$) and score(theta, y) (returns the gradient of the log-likelihood of $y$ at $\theta$). If either is missing, an error is raised.
  • n_sim — Positive integer: the number of simulated data sets averaged per evaluation. Should comfortably exceed dim to ensure the raw estimate is non-singular before the eigenvalue floor is applied.
  • seed — Positive integer: the base seed combined with the position key to produce the sub-seed for each evaluation.
  • floor — Non-negative scalar: minimum eigenvalue imposed on the averaged matrix to guarantee strict positive-definiteness (required downstream by the log-Cholesky map).

Mathematics

The expected Fisher information is

$$I(\theta) = \mathbb{E}_{y \sim p(\cdot \mid \theta)}!\left[,s(\theta, y),s(\theta, y)^\top\right]$$

where $s(\theta, y) = \nabla_\theta \log p(y \mid \theta)$ is the score. The estimator is the sample average of $n_{\text{sim}}$ rank-one outer products:

$$\hat{I}(\theta) = \frac{1}{n_{\text{sim}}} \sum_{i=1}^{n_{\text{sim}}} s(\theta, y_i),s(\theta, y_i)^\top, \qquad y_i \overset{\text{iid}}{\sim} p(\cdot \mid \theta).$$

This is positive semi-definite by construction (sum of PSD rank-one terms) and unbiased. The eigenvalue floor is applied as

$$\hat{I}_{\text{floor}} = U,\max(\Lambda, \texttt{floor}\cdot I),U^\top$$

via .gdpar_geom_floor_spd.

The sub-seed for reproducibility is

$$\text{sub} = \bigl(\texttt{base_seed} + \texttt{.gdpar_geom_seed_from_theta}(\theta)\bigr) \bmod 2147483646 + 1.$$

Returns A closure fisher(theta) returning a dim × dim symmetric strictly-positive-definite matrix (the estimated expected Fisher), with attribute "n_sim" recording the number of simulations.

Notes

  • Each call saves and restores the global RNG state (.Random.seed) so that evaluations are independent of call order — the position key makes the function deterministic in $\theta$.
  • The docstring explicitly notes that antithetic simulation would give no variance reduction for location families because the score is odd in the centred data, hence $(-s)(-s)^\top = ss^\top$.
  • Input validation raises "gdpar_input_error" if simulate or score are missing from the target, or if the score output is not a finite vector of length dim.

.gdpar_geom_subriemann_flow(theta, p, ref, U, omega, t)

Purpose Performs the exact, closed-form harmonic (Gaussian) flow of the reference quadratic used as the middle step of the sub-Riemannian Strang integrator. This is the exact symplectic rotation in each eigenmode, replacing the leapfrog step for the stiff wall modes.

Arguments

  • theta — Numeric vector of length $d$: the current position.
  • p — Numeric vector of length $d$: the current momentum.
  • ref — Numeric vector of length $d$: the reference position $\theta^\ast$.
  • U — Orthogonal $d \times d$ matrix: the eigenvectors of the reference Fisher (columns are eigenvectors), defining the modal coordinates $z = U^\top(\theta - \text{ref})$, $q = U^\top p$.
  • omega — Numeric vector of length $d$: the modal frequencies $\omega_i = \sqrt{w_i \lambda_i}$ where $w_i$ is the verticality weight and $\lambda_i$ the Fisher eigenvalue.
  • t — Scalar: the flow time (half the leapfrog step, or the full step for the exact flow).

Mathematics

Each mode $i$ evolves under the Hamiltonian $H_i = \frac{1}{2}q_i^2 + \frac{1}{2}\omega_i^2 z_i^2$ with exact flow

$$\begin{pmatrix} z_i(t) \ q_i(t) \end{pmatrix} = \begin{pmatrix} \cos(\omega_i t) & \sin(\omega_i t)/\omega_i \ -\omega_i\sin(\omega_i t) & \cos(\omega_i t) \end{pmatrix} \begin{pmatrix} z_i(0) \ q_i(0) \end{pmatrix}.$$

For $\omega_i \leq 10^{-8}$ (degenerate "canyon floor" mode), the flow continuously degrades to the free drift

$$\begin{pmatrix} z_i(t) \ q_i(t) \end{pmatrix} = \begin{pmatrix} 1 & t \ 0 & 1 \end{pmatrix} \begin{pmatrix} z_i(0) \ q_i(0) \end{pmatrix}.$$

The implementation uses ifelse(drift, ...) element-wise to blend between the two regimes without branching.

In global coordinates:

$$\theta_{\text{new}} = \text{ref} + U,\mathbf{z}_{\text{new}}, \qquad p_{\text{new}} = U,\mathbf{q}_{\text{new}}.$$

Returns A named list with components theta (numeric vector, the advanced position) and p (numeric vector, the advanced momentum).

Notes

  • The flow is exact and reversible for any time step $t$; there is no step-size penalty from stiff (high-$\omega$) modes, which is the core advantage of the sub-Riemannian integrator.
  • The ifelse on the drift mask ensures numerical stability as $\omega_i \to 0$: $\sin(\omega t)/\omega \to t$ and $\omega\sin(\omega t) \to 0$.

gdpar_geom_metric_subriemannian(target, fisher, reference = NULL, tau = NULL, softness = 1, floor = 1e-8)

Purpose Constructs the level-5 geometry of the Block RG hierarchy: a sub-Riemannian metric for quasi-deterministic posteriors where the typical set contracts onto a lower-dimensional manifold. The expected Fisher information grows without bound along stiff "wall" directions while soft "floor" directions carry genuine variation. The metric equips only a distribution $D_\theta \subsetneq T_\theta M$ (the near-null space of the Fisher) with an inner product, so the sampler glides along the floor instead of fighting the walls.

Arguments

  • target — A gdpar_geom_target (or any object accepted by its constructor). Converted internally.
  • fisher — A function fisher(theta) returning the dim × dim expected Fisher information matrix at position theta. For models without a closed form, use gdpar_geom_fisher_simulator. Evaluated once at the reference position.
  • reference — Optional numeric vector of length dim at which the Fisher and the reference quadratic are evaluated. Defaults to the zero vector.
  • tau — Positive scalar: the eigenvalue threshold of the verticality filter (directions with eigenvalue $\lambda \gg \tau$ are walls). Defaults to the smallest eigenvalue of the (floored) Fisher at the reference, so that the leapfrog residual curvature is capped at the floor scale and the step is never limited by stiff walls.
  • softness — Positive scalar: the logistic width of the verticality filter in log-eigenvalue space. Small values approximate a hard floor/wall split; the default 1 keeps the blend smooth.
  • floor — Non-negative scalar: minimum eigenvalue imposed on the Fisher before the filter is applied, so genuinely flat directions remain well-defined.

Mathematics

At the reference position $\theta^\ast$, the expected Fisher is eigendecomposed:

$$I(\theta^*) = U,\Lambda,U^\top, \qquad \Lambda = \mathrm{diag}(\lambda_1, \ldots, \lambda_d).$$

A continuous verticality filter assigns each eigendirection a weight in $(0,1)$:

$$w_i = \sigma!\left(\frac{\log\lambda_i - \log\tau}{s}\right) = \frac{1}{1 + \exp!\left(-\frac{\log\lambda_i - \log\tau}{s}\right)}$$

where $s$ is the softness parameter. Directions with $\lambda_i \gg \tau$ get $w_i \approx 1$ (walls); directions with $\lambda_i \ll \tau$ get $w_i \approx 0$ (floor).

The wall curvature matrix is

$$A = U,\mathrm{diag}(w_i,\lambda_i),U^\top$$

with modal frequencies $\omega_i = \sqrt{w_i,\lambda_i}$.

The integrator uses a Strang splitting with a Euclidean kinetic energy $K(p) = \frac{1}{2}|p|^2$:

  1. Half momentum kick: $p \leftarrow p + \frac{\epsilon}{2},\nabla U_{\text{rest}}(\theta)$ where $U_{\text{rest}} = -\log p(\theta) - \frac{1}{2}(\theta - \theta^\ast)^\top A,(\theta - \theta^\ast)$.
  2. Exact harmonic flow of the reference quadratic for time $\epsilon$: .gdpar_geom_subriemann_flow.
  3. Second half momentum kick.

The residual gradient is

$$\nabla U_{\text{rest}}(\theta) = -\nabla\log p(\theta) - A(\theta - \theta^*).$$

Returns A list of class c("gdpar_geom_metric", "list") with position_dependent = FALSE and metric_kind = "sub_riemannian". The mass functions (mass, inv_mass, chol_mass, logdet) all return the identity (the kinetic energy is Euclidean). An integrator closure is provided for consumption by gdpar_geom_hmc. Diagnostic fields include: reference, eigenvalues, verticality (the $w_i$), frequencies (the $\omega_i$), n_walls (count of directions with $w_i &gt; 0.5$), tau, softness, and suggested_epsilon (a leapfrog step matched to the floor scale).

Notes

  • The Metropolis correction inside gdpar_geom_hmc uses the exact target log-density, so the sampler is exact regardless of how coarse the Gaussian wall approximation is — only efficiency depends on it.
  • The suggested_epsilon exploits the fact that the exact wall integration imposes no step-size limit; the step is set solely by the gentle floor curvature, which is the source of the speed-up over Euclidean HMC (whose step is bottlenecked by the stiffest wall eigenvalue).
  • The scheme is symplectic and time-reversible by construction (split HMC with a Gaussian part; Shahbaba et al. 2014).

gdpar_geom_metric_subriemannian(target, fisher, reference = NULL, tau = NULL, softness = 1, floor = 1e-8)

Purpose Builds the level-3 "sub-Riemannian" metric of the Block RG hierarchy. It partitions the eigen-spectrum of a reference Fisher information matrix into "wall" directions (curvature above threshold $\tau$, handled exactly) and residual directions (curvature below $\tau$, handled by the leapfrog), producing a geometry whose integrator can take step sizes governed by the soft floor rather than the stiffest eigenvalue. This is the remedy for extreme eigenvalue heterogeneity (the count/regime-switching target class).

Arguments

Argument Type Meaning
target gdpar_geom_target or coercible The target distribution. If not already a gdpar_geom_target, it is coerced via gdpar_geom_target(). Provides $dim, $grad_log_prob.
fisher function A function theta ↦ I(θ) returning the expected Fisher information matrix at theta. Must return a d × d symmetric positive-(semi)definite matrix.
reference numeric vector or NULL The reference parameter point at which fisher is evaluated to determine the spectral decomposition. Defaults to the zero vector of length d.
tau numeric scalar or NULL Spectral threshold separating walls from residuals. Eigenvalues $\lambda_k \ge \tau$ become exact wall potentials; eigenvalues $\lambda_k &lt; \tau$ are residual and integrated by the leapfrog. Defaults to min(λ) (the softest direction), meaning everything is a residual (no walls).
softness numeric scalar ≥ 0 Controls the sharpness of the sigmoid that assigns each eigenvalue to the wall or residual partition. Smaller values approach a hard threshold; larger values make a smoother blend.
floor numeric scalar ≥ 0 Eigenvalue floor; any eigenvalue of the reference Fisher below this is clamped. Also used by .gdpar_geom_floor_spd when flooring the Fisher matrix itself.

Mathematics

The reference Fisher matrix at the reference point is $$I_{\text{ref}} = \mathrm{\texttt{fisher}}(\text{ref}),$$ floored to be strictly positive-definite. Its eigendecomposition is $$I_{\text{ref}} = U ,\mathrm{diag}(\lambda_1,\dots,\lambda_d), U^\top,$$ with $\lambda_k \leftarrow \max(\lambda_k, \texttt{floor})$.

Verticality weights (logistic sigmoid): $$w_k = \frac{1}{1 + \exp!\bigl(-(\log \lambda_k - \log \tau),/,\texttt{softness}\bigr)}.$$

Wall curvature and frequencies: $$a_k = w_k ,\lambda_k, \qquad \omega_k = \sqrt{\max(a_k, 0)}.$$

Wall Hessian matrix: $$A = U,\mathrm{diag}(a_1,\dots,a_d),U^\top.$$

Residual potential gradient. The function returns a gradient closure grad_U_rest(θ) that computes $$\nabla U_{\text{rest}}(\theta) = -\nabla \log p(\theta) - A,(\theta - \text{ref}),$$ i.e. the full negative log-probability gradient minus the reference-quadratic wall term that the exact flow already integrates.

Strang (leapfrog + exact flow) integrator. For $L$ steps of size $\varepsilon$:

For each step $\ell = 1,\dots,L$:

  1. Half kick: $;p \leftarrow p - \frac{\varepsilon}{2},\nabla U_{\text{rest}}(\theta)$
  2. Exact sub-Riemannian flow via .gdpar_geom_subriemann_flow(θ, p, ref, U, ω, ε) for time $\varepsilon$ in the wall directions.
  3. Half kick: $;p \leftarrow p - \frac{\varepsilon}{2},\nabla U_{\text{rest}}(\theta)$

If any iterate is non-finite the flag converged is set to FALSE (the proposal will be rejected).

Suggested step size. The residual maximum curvature is $$\lambda_{\text{resid}}^{\max} = \max_k \bigl((1 - w_k),\lambda_k,, 0\bigr).$$ The suggested epsilon keeping the residual leapfrog stable is $$\varepsilon_{\text{suggested}} = \begin{cases} 2 & \text{if } \lambda_{\text{resid}}^{\max} \le 0,\[4pt] \min!\bigl(2,; 0.7 / \sqrt{\lambda_{\text{resid}}^{\max}}\bigr) & \text{otherwise.}\end{cases}$$

Returns A list of class c("gdpar_geom_metric", "list") with:

Field Value
position_dependent FALSE
dim integer, parameter dimension $d$
metric_kind "sub_riemannian"
mass(θ), inv_mass(θ), chol_mass(θ) Identity matrix $I_d$ (mass matrix is trivial; all action is in the potential)
logdet(θ) 0
integrator Closure (theta, p, target, eps, L) → list(theta, p, converged) implementing the Strang scheme above
reference numeric vector, the reference point
eigenvalues numeric vector $\lambda$
verticality numeric vector $w$
frequencies numeric vector $\omega$
n_walls integer, sum(w > 0.5)
tau numeric, the spectral threshold
softness numeric, the sigmoid softness
suggested_epsilon numeric, the suggested step size

Notes

  • Raises gdpar_input_error if fisher is not a function, if reference has wrong length, or if fisher does not return a d × d matrix.
  • softness = 0 would yield a degenerate sigmoid; the code does not guard against this (the logistic denominator becomes 1 + exp(±∞)).
  • When tau = min(λ) (default), $w_k \approx 0.5$ for the softest direction and $w_k \approx 1$ for stiffer ones; the exact assignment depends on softness.
  • The mass matrix is identity because the sub-Riemannian geometry encodes stiffness entirely through the potential, not through a position-dependent kinetic energy.

.gdpar_geom_relativistic_radial(d, speed, mass, n_grid = 8192L)

Purpose Internal factory that builds a pre-computed inverse-CDF sampler for the radial component of relativistic momentum draws. Under the change of variables $p = L q$ (where $M = LL^\top$), the radius $r = |q|$ has a density proportional to $r^{d-1}\exp!\bigl(-c\sqrt{r^2 + m^2 c^2}\bigr)$ independent of $\theta$. The sampler is constructed once and reused for every momentum refresh.

Arguments

Argument Type Meaning
d integer Dimension of the parameter space.
speed numeric > 0 Speed of light $c$.
mass numeric > 0 Rest mass $m$.
n_grid integer Number of grid points for the numerical CDF approximation (default 8192).

Mathematics

The unnormalised log-density of the radius $r$ is $$\log g(r) = (d-1)\log!\max(r,,10^{-300}) - c\sqrt{r^2 + m^2 c^2}.$$

The algorithm:

  1. Maximise $\log g$ on $[10^{-8}, 10^6]$ to find the peak value.
  2. Extend the upper grid bound $r_{\text{hi}}$ until $\log g(r_{\text{hi}}) &lt; \text{peak} - 60$ (or $r_{\text{hi}} = 10^9$).
  3. Build a grid $r_0 = 0, r_1, \dots, r_{n-1} = r_{\text{hi}}$.
  4. Evaluate $\log g$ on the grid; set $\log g(0) = -\infty$ when $d &gt; 1$ (since $r^{d-1} = 0$).
  5. Compute the CDF via the trapezoidal rule, normalise to 1.
  6. Build a monotone linear interpolation inv mapping $[0, 1] \to [0, r_{\text{hi}}]$ (only on strictly-increasing CDF segments).

Returns A function function(n) that draws n samples from the radial distribution via inverse-CDF transform of runif(n).

Notes

  • When $d = 1$, the first grid point lg[1] is kept (not set to $-\infty$), because $r^{d-1} = r^0 = 1$ at $r = 0$.
  • The pmax(r, 1e-300) guards against log(0) on the grid.
  • The CDF interpolation uses rule = 2 (extrapolation by boundary values) for numerical safety at the extremes.
  • The keep mask ensures only strictly-increasing CDF entries feed into approxfun, preventing non-invertibility.
  • The result is independent of $\theta$ and is built once per call to .gdpar_geom_kinetic_relativistic.

.gdpar_geom_kinetic_relativistic(metric, speed, mass)

Purpose Internal factory that constructs the relativistic kinetic energy object: its value, its partial gradients with respect to $p$ and $\theta$, and a momentum sampler. This is the non-separable kinetic energy that replaces the Gaussian kinetic of Riemannian HMC.

Arguments

Argument Type Meaning
metric gdpar_geom_metric A Riemannian metric object providing $dim, $inv_mass(θ), $logdet(θ), $chol_mass(θ), $dmass(θ), $position_dependent.
speed numeric > 0 Speed of light $c$.
mass numeric > 0 Rest mass $m$.

Mathematics

Define $s(\theta, p) = \sqrt{p^\top M(\theta)^{-1} p + m^2 c^2}$. The kinetic energy is $$K(\theta, p) = c,s(\theta, p) + \tfrac{1}{2}\log\det M(\theta).$$

Gradient w.r.t. $p$ (the relativistic velocity): $$\nabla_p K = \frac{c,M(\theta)^{-1} p}{s(\theta, p)}.$$

Its $M$-norm is $c\sqrt{p^\top M^{-1}p},/,s &lt; c$, which is the bounded-velocity property.

Gradient w.r.t. $\theta$: $$\frac{\partial K}{\partial \theta_i} = -\frac{c,p^\top M^{-1},\frac{\partial M}{\partial \theta_i},M^{-1} p}{2,s(\theta, p)} + \frac{1}{2}\mathrm{tr}!\left(M^{-1}\frac{\partial M}{\partial \theta_i}\right).$$

When metric$position_dependent is FALSE, this gradient is identically zero.

Momentum sampler. Under $p = Lq$ with $q = r,u$, $u$ uniform on $S^{d-1}$ and $r$ from .gdpar_geom_relativistic_radial: $$p = L,(r \cdot u), \quad u = z/|z|, \quad z \sim \mathcal{N}(0, I_d).$$ Since $p^\top M^{-1} p = q^\top q = r^2$, the joint density of $p$ is exactly $\exp(-K(\theta, p))$.

Returns A named list with four closures:

Field Signature Returns
value (theta, p) → numeric $K(\theta, p)$
grad_p (theta, p) → numeric vector $\nabla_p K$
grad_theta (theta, p) → numeric vector $\nabla_\theta K$
draw_momentum (theta) → numeric vector One sample from $e^{-K(\theta,\cdot)}$

Notes

  • The precomputed radial sampler (radial) is captured in the closure environment; it is built exactly once.
  • metric$dmass(theta) must return a length-$d$ list of $d\times d$ matrices $\partial M/\partial\theta_i$. This is only called when metric$position_dependent is TRUE.
  • The 0.5 * metric$logdet(theta) term in value is the same normaliser as in Gaussian RMHMC; it cancels in the Metropolis acceptance ratio, ensuring exactness.

.gdpar_geom_relativistic_integrator(kinetic)

Purpose Internal factory that builds the dedicated generalised implicit leapfrog integrator for the non-separable relativistic Hamiltonian $H(\theta, p) = -\log p(\theta) + K(\theta, p)$. This integrator is carried in the metric's integrator slot and dispatched by gdpar_geom_hmc when a relativistic geometry is present.

Arguments

Argument Type Meaning
kinetic list of closures The kinetic energy object returned by .gdpar_geom_kinetic_relativistic, providing grad_p, grad_theta.

Mathematics

The integrator follows Girolami & Calderhead (2011), performing three sub-steps per leapfrog iteration (each iteration is time-reversible and volume-preserving to the fixed-point tolerance). Let $\nabla_\theta H = -\nabla\log p(\theta) + \nabla_\theta K(\theta, p)$.

For each of $L$ iterations:

Sub-step 1 — Implicit half kick in $p$. Solve by fixed-point iteration: $$p_{1/2} = p - \frac{\varepsilon}{2},\nabla_\theta H(\theta, p_{1/2}).$$ Iterate $p^{(k+1)} = p - \frac{\varepsilon}{2},\nabla_\theta H(\theta, p^{(k)})$ until $\|p^{(k+1)} - p^{(k)}\|_\infty &lt; \texttt{fp\_tol}$ or fp_max iterations.

Sub-step 2 — Implicit drift in $\theta$. The velocity is $v = \nabla_p K$. Solve by fixed-point iteration: $$\theta' = \theta + \frac{\varepsilon}{2}\bigl(v(\theta, p_{1/2}) + v(\theta', p_{1/2})\bigr).$$ Iterate until convergence or fp_max iterations.

Sub-step 3 — Explicit half kick in $p$: $$p' = p_{1/2} - \frac{\varepsilon}{2},\nabla_\theta H(\theta', p_{1/2}).$$

Returns A closure function(theta, p, target, eps, L, fp_tol = 1e-9, fp_max = 100L) that returns list(theta, p, converged).

Notes

  • If either fixed-point iteration (sub-steps 1 or 2) does not converge, or any iterate is non-finite, converged is set to FALSE and the loop breaks early. The caller (the Metropolis step in gdpar_geom_hmc) will reject the proposal.
  • The entire integration is wrapped in tryCatch; any runtime error returns converged = FALSE with the last valid (theta, p).
  • fp_tol = 1e-9 and fp_max = 100 are defaults; the caller may override them.
  • The implicit scheme is necessary because $H$ is non-separable: the velocity $\nabla_p K$ depends on $\theta$ through $M(\theta)$, so the standard explicit leapfrog is not symplectic.
  • Bit-identical behaviour with the default (Gaussian) leapfrog is achieved by having gdpar_geom_hmc dispatch to this integrator only when the metric provides one, leaving the default branch unchanged.

gdpar_geom_metric_relativistic(target, curvature = c("fisher", "softabs"), fisher = NULL, dfisher = NULL, speed = 10, rest_mass = 1, alpha = 1e6, floor = 1e-8, fd_step = 1e-4)

Purpose Constructs a relativistic metric object that augments a Riemannian metric (Euclidean, Fisher, or SoftAbs) with a relativistic kinetic energy. This corresponds to the Finsler/relativistic level of the Block RG hierarchy. The metric delegates the position-dependent mass matrix and its derivatives to the underlying Riemannian layer, but wraps a relativistic kinetic energy and a corresponding implicit integrator on top.

Arguments

  • target — An object of class gdpar_geom_target or a list/closure coercible to one via gdpar_geom_target(target).
  • curvature — Character string, one of "fisher" or "softabs". Selects the curvature scheme for the underlying Riemannian metric. Passed through to gdpar_geom_metric_riemannian.
  • fisher — An optional user-supplied Fisher information function fisher(theta) returning the $d \times d$ Fisher matrix at position $\theta$. If NULL, the Riemannian metric constructor supplies its own.
  • dfisher — An optional user-supplied derivative-of-Fisher function dfisher(theta) returning a list of $d$ matrices $\partial G / \partial \theta_i$. If NULL, finite-differenced.
  • speed — Numeric scalar, $c &gt; 0$. The relativistic "speed of light" constant controlling the saturation of momentum magnitude.
  • rest_mass — Numeric scalar, $m_0 &gt; 0$. The relativistic rest mass constant.
  • alpha — Numeric scalar. SoftAbs regularization parameter forwarded to the Riemannian metric.
  • floor — Numeric scalar. Eigenvalue floor for the Riemannian metric. Forwarded to the Riemannian metric.
  • fd_step — Numeric scalar. Finite-difference step size for numerically differentiating the Fisher matrix when dfisher is NULL. Forwarded to the Riemannian metric.

Mathematics

The relativistic kinetic energy is:

$$K(\theta, p) = m_0 c^2 \left(\sqrt{1 + \frac{p^\top M(\theta)^{-1} p}{m_0^2 c^2}} - 1\right)$$

where $M(\theta)$ is the position-dependent mass matrix from the underlying Riemannian metric. This saturates the effective velocity $v = M^{-1} p$ at the speed of light $c$, providing natural gradient clipping that prevents the explosive trajectories that plague HMC on poorly conditioned geometries.

The Hamiltonian is $H(\theta, p) = U(\theta) + K(\theta, p)$ where $U(\theta) = -\log p(\theta)$.

Returns A list of class c("gdpar_geom_metric", "list") with components:

Component Description
position_dependent TRUE (relativistic metrics are always position-dependent via $M(\theta)$).
dim Integer, the parameter dimension $d$.
metric_kind "relativistic".
speed The speed constant $c$.
rest_mass The rest mass $m_0$.
curvature The curvature scheme string.
mass Function mass(theta) returning $M(\theta)$; inherited from the Riemannian base.
inv_mass Function inv_mass(theta) returning $M(\theta)^{-1}$.
chol_mass Function chol_mass(theta) returning the Cholesky factor $L$ of $M(\theta)$ (so $M = L L^\top$).
logdet Function logdet(theta) returning $\log \det M(\theta)$.
dmass Function dmass(theta) returning the list $[\partial M / \partial \theta_1, \ldots, \partial M / \partial \theta_d]$.
kinetic A list of closures for the relativistic kinetic energy (see .gdpar_geom_kinetic_relativistic).
integrator A closure for the relativistic generalised leapfrog integrator (see .gdpar_geom_relativistic_integrator).

Notes

  • speed and rest_mass are validated to be strictly positive; a zero or negative value raises a gdpar_input_error via gdpar_abort.
  • The underlying Riemannian metric is constructed via gdpar_geom_metric_riemannian, passing all curvature-related arguments through.
  • The kinetic and integrator fields are populated by two internal helpers (.gdpar_geom_kinetic_relativistic and .gdpar_geom_relativistic_integrator), which are not defined in this section — they are defined elsewhere in the package.
  • All mass-matrix accessor functions (mass, inv_mass, chol_mass, logdet, dmass) are references to the Riemannian base metric's implementations.

.gdpar_geom_kinetic_gaussian(metric)

Purpose Constructs the standard Gaussian kinetic energy closure list used by the Euclidean and Riemannian HMC levels. The Gaussian kinetic energy $K(\theta, p) = \tfrac{1}{2} p^\top M(\theta)^{-1} p + \tfrac{1}{2} \log \det M(\theta)$ is the canonical momentum distribution when $p \mid \theta \sim \mathcal{N}(0, M(\theta))$. This function is internal (leading dot) and is called by gdpar_geom_hmc when the metric does not supply its own kinetic.

Arguments

  • metric — A metric object (typically of class gdpar_geom_metric) exposing inv_mass, logdet, dmass, chol_mass, and the boolean position_dependent.

Mathematics

The kinetic energy and its derivatives are:

$$K(\theta, p) = \frac{1}{2} p^\top M(\theta)^{-1} p + \frac{1}{2} \log \det M(\theta)$$

$$\frac{\partial K}{\partial p} = M(\theta)^{-1} p$$

$$\frac{\partial K}{\partial \theta_i} = \frac{1}{2} \mathrm{tr}!\bigl(M^{-1} \tfrac{\partial M}{\partial \theta_i}\bigr) - \frac{1}{2} p^\top M^{-1} \tfrac{\partial M}{\partial \theta_i} M^{-1} p$$

The gradient with respect to $p$ is the velocity $v = M^{-1} p$. The gradient with respect to $\theta$ vanishes for a constant (Euclidean) metric ($\partial M / \partial \theta_i = 0$).

Momentum is drawn as $p = L , z$ where $L$ is the Cholesky factor of $M(\theta)$ and $z \sim \mathcal{N}(0, I)$, ensuring $p \sim \mathcal{N}(0, M(\theta))$.

Returns A list with four named function components:

Component Signature Returns
value (theta, p) Numeric scalar, the kinetic energy $K(\theta, p)$.
grad_p (theta, p) Numeric vector of length $d$, $\partial K / \partial p = M(\theta)^{-1} p$.
grad_theta (theta, p) Numeric vector of length $d$, $\partial K / \partial \theta$.
draw_momentum (theta) Numeric vector of length $d$, a draw from $\mathcal{N}(0, M(\theta))$.

Notes

  • When metric$position_dependent is FALSE, grad_theta short-circuits to rep(0, length(theta)), avoiding unnecessary computation.
  • The grad_theta computation uses vapply over the dimension index $i$, indexing into the list returned by metric$dmass(theta). Each element $dMi$ is a $d \times d$ matrix.
  • The expression sum(Minv * dMi) computes $\mathrm{tr}(M^{-1} , dM_i)$ efficiently via elementwise multiplication and summation, exploiting the fact that both matrices are the same size.
  • The expression as.numeric(crossprod(Minv_p, dMi %*% Minv_p)) computes $p^\top M^{-1} dM_i M^{-1} p$ via $\texttt{Minv\_p}^\top \, dM_i \, \texttt{Minv\_p}$.

.gdpar_geom_dH_dtheta(target, kinetic, theta, p)

Purpose Computes the gradient of the Hamiltonian $H$ with respect to the position $\theta$, used inside the leapfrog integrator. Internal helper.

Arguments

  • target — A gdpar_geom_target providing grad_log_prob(theta).
  • kinetic — A kinetic energy list (as returned by .gdpar_geom_kinetic_gaussian or .gdpar_geom_kinetic_relativistic) providing grad_theta(theta, p).
  • theta — Numeric vector, the current position.
  • p — Numeric vector, the current momentum.

Mathematics

$$\frac{\partial H}{\partial \theta} = -\nabla_\theta \log p(\theta) + \frac{\partial K}{\partial \theta}$$

Since $H(\theta, p) = U(\theta) + K(\theta, p)$ with $U = -\log p(\theta)$, the gradient is the negative target gradient plus the kinetic position gradient. For a constant Euclidean metric, $\partial K / \partial \theta = 0$ and the expression reduces to $-\nabla_\theta \log p(\theta)$.

Returns A numeric vector of length $d$ equal to $\nabla_\theta H(\theta, p)$.

Notes

  • This is a thin composition; all heavy lifting is delegated to target$grad_log_prob and kinetic$grad_theta.

.gdpar_geom_leapfrog_step(theta, p, target, metric, kinetic, eps, fp_tol = 1e-9, fp_max = 100L)

Purpose Performs a single leapfrog integration step. Dispatches to an explicit Stoermer–Verlet step for constant metrics or to the implicit generalised leapfrog of Girolami & Calderhead (2011) for position-dependent (Riemannian/relativistic) metrics. Internal helper.

Arguments

  • theta — Numeric vector of length $d$, current position.
  • p — Numeric vector of length $d$, current momentum.
  • target — A gdpar_geom_target providing grad_log_prob.
  • metric — A metric object providing inv_mass, position_dependent, and optionally integrator.
  • kinetic — A kinetic energy list providing grad_p and grad_theta.
  • eps — Numeric scalar, the leapfrog step size $\varepsilon$.
  • fp_tol — Numeric scalar, convergence tolerance $\delta$ for the fixed-point iterations. Only used for position-dependent metrics.
  • fp_max — Integer, maximum iterations per fixed-point solve. Only used for position-dependent metrics.

Mathematics

Explicit (Euclidean/constant metric) step — the standard Stoermer–Verlet integrator:

$$p_{1/2} = p_0 + \frac{\varepsilon}{2} \nabla_\theta \log p(\theta_0)$$

$$\theta_1 = \theta_0 + \varepsilon , M^{-1} p_{1/2}$$

$$p_1 = p_{1/2} + \frac{\varepsilon}{2} \nabla_\theta \log p(\theta_1)$$

Implicit (position-dependent metric) step — three sub-steps with two fixed-point solves:

Sub-step 1 (implicit half-kick in $p$): Iterate to convergence:

$$p^{(k+1)} = p_0 - \frac{\varepsilon}{2} \nabla_\theta H(\theta_0, p^{(k)})$$

converging to $p_{1/2}$ such that $p_{1/2} = p_0 - \tfrac{\varepsilon}{2} \nabla_\theta H(\theta_0, p_{1/2})$.

Sub-step 2 (implicit full drift in $\theta$): Iterate to convergence:

$$\theta^{(k+1)} = \theta_0 + \frac{\varepsilon}{2}\bigl(M(\theta_0)^{-1} + M(\theta^{(k)})^{-1}\bigr) p_{1/2}$$

converging to $\theta_1$.

Sub-step 3 (explicit half-kick in $p$):

$$p_1 = p_{1/2} - \frac{\varepsilon}{2} \nabla_\theta H(\theta_1, p_{1/2})$$

Returns A list with:

Component Type Description
theta Numeric vector Updated position after one step.
p Numeric vector Updated momentum after one step.
converged Logical TRUE if both fixed-point solves converged (always TRUE for the explicit branch); FALSE otherwise. Present only for the implicit branch.

Notes

  • The constant-metric branch returns a two-element list (theta, p) with no converged field.
  • The implicit branch is wrapped in tryCatch; any error (e.g., non-finite iterates, matrix overflow) is caught, and the function returns the original (theta, p) with converged = FALSE. This causes the Metropolis correction to reject the proposal rather than crashing the sampler.
  • Non-finite delta in any fixed-point loop immediately breaks out of that loop, leaving converged = FALSE.
  • The integrator is exactly time-reversible and volume-preserving (symplectic) up to the fixed-point tolerance.
  • When metric$integrator is non-NULL (as in the relativistic case), gdpar_geom_hmc bypasses this function entirely and calls the metric's own integrator.

.gdpar_geom_leapfrog_traj(theta, p, target, metric, kinetic, eps, L, fp_tol = 1e-9, fp_max = 100L)

Purpose Chains $L$ consecutive leapfrog steps to form a full trajectory. Returns the endpoint state and a global convergence flag. Internal helper.

Arguments

  • theta — Numeric vector of length $d$, initial position.
  • p — Numeric vector of length $d$, initial momentum.
  • target — A gdpar_geom_target.
  • metric — A metric object.
  • kinetic — A kinetic energy list.
  • eps — Numeric scalar, leapfrog step size $\varepsilon$.
  • L — Integer, number of leapfrog steps $L$.
  • fp_tol — Numeric scalar, fixed-point convergence tolerance.
  • fp_max — Integer, maximum fixed-point iterations.

Mathematics

Iterates the leapfrog map $\phi_\varepsilon$ for $L$ steps:

$$(\theta_L, p_L) = \phi_\varepsilon^L(\theta_0, p_0)$$

Returns A list with:

Component Type Description
theta Numeric vector Position after $L$ steps.
p Numeric vector Momentum after $L$ steps.
converged Logical TRUE if every step converged; FALSE if any step failed (in which case iteration halts early at the failed step's output).

Notes

  • If any step reports converged = FALSE, the loop breaks immediately and the remaining steps are not attempted. The returned state is that of the last (failed) step.
  • For a constant metric (explicit leapfrog), every step's converged is effectively always TRUE.

.gdpar_geom_hamiltonian(target, kinetic, theta, p)

Purpose Evaluates the Hamiltonian $H(\theta, p) = U(\theta) + K(\theta, p)$ used for the Metropolis accept/reject step. Internal helper.

Arguments

  • target — A gdpar_geom_target providing log_prob(theta).
  • kinetic — A kinetic energy list providing value(theta, p).
  • theta — Numeric vector, position.
  • p — Numeric vector, momentum.

Mathematics

$$H(\theta, p) = -\log p(\theta) + K(\theta, p)$$

where $K$ is either the Gaussian kinetic $K = \tfrac{1}{2} p^\top M^{-1} p + \tfrac{1}{2} \log \det M$ or the relativistic kinetic $K = m_0 c^2 (\sqrt{1 + p^\top M^{-1} p / (m_0^2 c^2)} - 1)$, depending on which kinetic object is passed in.

Returns A numeric scalar, the Hamiltonian energy at $(\theta, p)$.

Notes

  • This function does not guard against non-finite return values; callers (gdpar_geom_hmc) wrap the call in tryCatch and treat Inf as a divergent proposal.

gdpar_geom_hmc(target, metric = NULL, epsilon = 0.1, L = 20L, n_iter = 1000L, n_warmup = 500L, init = NULL, seed = NULL, fp_tol = 1e-9, fp_max = 100L)

Purpose The main exported HMC sampler implementing the Block RG engine. Runs a fixed step-size, fixed trajectory-length Hamiltonian Monte Carlo with Metropolis correction over a pluggable metric geometry. Supports Euclidean, Riemannian (position-dependent), and (via metric$integrator) relativistic/sub-Riemannian levels of the hierarchy. No adaptation is performed; warmup iterations are simply discarded.

Arguments

  • target — A gdpar_geom_target or a list/closure coercible to one.
  • metric — A metric object (e.g., from gdpar_geom_metric_euclidean, gdpar_geom_metric_riemannian, gdpar_geom_metric_relativistic). If NULL, defaults to an identity Euclidean metric of dimension target$dim.
  • epsilon — Numeric scalar $&gt; 0$. Leapfrog step size $\varepsilon$.
  • L — Positive integer. Number of leapfrog steps per proposal.
  • n_iter — Positive integer. Number of retained (post-warmup) iterations.
  • n_warmup — Non-negative integer. Number of warmup iterations (burn-in, no adaptation).
  • init — Optional numeric vector of length $d$ giving the initial position $\theta_0$. Defaults to rep(0, d).
  • seed — Optional integer. If supplied, the RNG state is saved, the seed is set, and the state is restored on exit via on.exit.
  • fp_tol — Numeric scalar. Fixed-point convergence tolerance for implicit leapfrog steps. Only used when the metric is position-dependent.
  • fp_max — Integer. Maximum fixed-point iterations per implicit sub-step.

Mathematics

The sampler implements the standard HMC algorithm:

  1. Draw momentum $p_0 \sim \mathcal{N}(0, M(\theta))$.
  2. Compute $H_0 = H(\theta, p_0)$.
  3. Integrate $L$ leapfrog steps to obtain proposal $(\theta^\ast, p^\ast)$.
  4. Compute $H_1 = H(\theta^\ast, p^\ast)$.
  5. Accept with probability $\min(1, \exp(-(H_1 - H_0)))$.

The energy difference $\Delta H = H_1 - H_0$ is the Metropolis–Hastings log acceptance ratio. A proposal is flagged as divergent if:

  • the implicit solve did not converge (converged == FALSE), or
  • $\Delta H$ is non-finite, or
  • $|\Delta H| &gt; 1000$.

The energy Bayesian fraction of missing information (E-BFMI) is computed as:

$$\text{E-BFMI} = \frac{\sum_{i} (\Delta E_i)^2}{\sum_i (E_i - \bar{E})^2}$$

where $\Delta E_i = E_{i+1} - E_i$ and $E_i$ are the Hamiltonian energies of the accepted states.

Returns A list of class c("gdpar_geom_hmc", "list") with:

Component Type Description
draws $n_{\text{iter}} \times d$ numeric matrix Post-warmup position draws. Column names from target$param_names.
accept_rate Numeric scalar Acceptance rate $= n_{\text{accept}} / (n_{\text{warmup}} + n_{\text{iter}})$.
n_divergent Integer Count of divergent proposals across all iterations.
energy Numeric vector of length $n_{\text{iter}}$ Hamiltonian trace (pre-acceptance $H_0$ for each post-warmup iteration).
ebfmi Numeric scalar Energy Bayesian fraction of missing information. NA_real_ if fewer than 2 finite energy values.
epsilon Numeric scalar The step size used.
L Integer The trajectory length used.
metric_type Character "euclidean_constant" if metric$position_dependent is FALSE, "position_dependent" otherwise.

Notes

  • The target is coerced to gdpar_geom_target if not already one.
  • epsilon must be strictly positive; validated via assert_numeric_scalar with lower = 0 (and then effectively checked by the leapfrog).
  • L, n_iter are validated via assert_count (must be positive integers).
  • n_warmup is validated manually to be a non-negative integer scalar.
  • If metric$integrator is a function (e.g., provided by a relativistic or sub-Riemannian metric), it is called in place of the default .gdpar_geom_leapfrog_traj. The Metropolis correction still uses the exact Hamiltonian from the kinetic energy object.
  • Seed handling: if seed is supplied and .Random.seed exists in .GlobalEnv, it is saved and restored via on.exit; otherwise, set.seed(seed) is called directly.
  • If metric$kinetic is non-NULL (e.g., the relativistic metric populates this), that kinetic is used; otherwise .gdpar_geom_kinetic_gaussian(metric) is the default.
  • The warmup phase simply runs iterations without recording draws or tracking energy; no adaptation of epsilon or L occurs.
  • A proposal that raises any error (e.g., target log_prob or grad_log_prob throwing on a non-finite unconstrained value, an overflow in the Hessian, a non-PD metric) is caught by the outer tryCatch in the main loop and treated as a divergent rejection.

.gdpar_geom_ebfmi(energy)

Purpose Computes the energy Bayesian fraction of missing information (E-BFMI), a diagnostic for how well the metric resolves the target geometry (Betancourt 2016). Internal helper.

Arguments

  • energy — Numeric vector of per-iteration Hamiltonian energies (may contain non-finite values).

Mathematics

$$\text{E-BFMI} = \frac{\sum_{i=1}^{n-1} (E_{i+1} - E_i)^2}{\sum_{i=1}^{n} (E_i - \bar{E})^2}$$

where the sum is over finite energies only. Low values (well below 1) indicate that the metric does not adequately capture the target's geometry, leading to slow exploration of the energy landscape.

Returns A numeric scalar. Returns NA_real_ if fewer than 2 finite energy values are available, or if the denominator is non-finite or non-positive.

Notes

  • Non-finite entries in energy are filtered out before computation via energy[is.finite(energy)].
  • A denominator of zero or negative (degenerate case: all energies identical) yields NA_real_.

print.gdpar_geom_hmc(x, ...)

Purpose S3 print method for objects of class gdpar_geom_hmc. Produces a concise summary to the console.

Arguments

  • x — An object of class gdpar_geom_hmc.
  • ... — Unused; present for S3 method compatibility.

Mathematics None.

Returns Invisibly returns x (via invisible(x)).

Notes

  • Prints the draw matrix dimensions, metric type, epsilon, L, acceptance rate, number of divergent transitions, and E-BFMI to stdout via cat.
  • The acceptance rate and E-BFMI are formatted with 3 significant digits.
  • This is an exported S3 method registered via @export.

gdpar_geom_rmhmc_adaptive(target, fisher, n_sim, n_sites_init, max_rounds, batch, n_add, decay, novelty_tol, warmup_metric, epsilon, L, n_iter, n_warmup, init, seed, alpha, nugget, lengthscale)

Purpose
Implements an adaptive Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) sampler that uses a Gaussian process (GP) surrogate to approximate the local Fisher information matrix. The algorithm sequentially builds and refines the GP metric by collecting novel parameter sites in a reservoir.

Arguments

  • target : an object of class "gdpar_geom_target" or any object coercible to one via gdpar_geom_target(). Defines the target distribution and its geometric properties.
  • fisher : a function of theta returning the Fisher information matrix at theta, or NULL. If NULL, a simulator (gdpar_geom_fisher_simulator) is constructed internally.
  • n_sim : integer. Number of draws used by the Fisher simulator when fisher = NULL.
  • n_sites_init : integer. Initial number of parameter sites collected in the reservoir.
  • max_rounds : integer. Maximum number of novelty‑driven adaptation rounds.
  • batch : integer. Number of HMC draws generated per adaptation round.
  • n_add : integer. Base number of novel sites to admit per round before decay.
  • decay : numeric (≥ 0). Multiplicative decay factor for the number of admitted sites per round.
  • novelty_tol : numeric (≥ 0). Novelty threshold; only draws with novelty above this value are eligible for admission.
  • warmup_metric : a metric object (e.g., from gdpar_geom_metric_riemannian) used for the initial reservoir generation. If NULL, a Riemannian metric with softabs curvature is constructed.
  • epsilon : numeric. Leapfrog step size for all HMC runs.
  • L : integer. Number of leapfrog steps per HMC proposal.
  • n_iter : integer. Number of post‑warmup draws in the final sampling phase.
  • n_warmup : integer. Number of warmup draws in the final sampling phase.
  • init : numeric vector of length target$dim giving the initial parameter vector. If NULL, defaults to a zero vector.
  • seed : integer or NULL. If non‑NULL, the random seed is set and restored on exit.
  • alpha : numeric. Tuning parameter for the softabs curvature and the GP metric.
  • nugget : numeric. Jitter added to the diagonal of the GP kernel matrix to ensure positive‑definiteness.
  • lengthscale : numeric or NULL. Length‑scale hyper‑parameter for the GP kernel. If NULL, it is estimated from the reservoir sites.

Mathematics
The algorithm proceeds in three phases:

  1. Cold‑start reservoir
    Collect $M_0 = \texttt{n\_sites\_init}$ parameter sites $\theta^{(1)},\dots,\theta^{(M_0)}$ by running HMC with the fixed warmup_metric. The Fisher information at each site is either provided by fisher(\theta) or approximated via the simulator. These sites form the initial reservoir $S_0$.

  2. Novelty‑driven adaptation
    For round $r = 1,\dots,R$ (where $R \le \texttt{max\_rounds}$):

    • Draw a batch of $B = \texttt{batch}$ samples using HMC with the current GP metric.
    • Compute the novelty $\nu(\theta)$ for each draw via the GP metric's novelty method.
    • Let $\nu_{\max} = \max{\nu(\theta): \theta \text{ in batch}}$. If $\nu_{\max} &lt; \texttt{novelty\_tol}$, stop.
    • Otherwise, select the $n_{\text{admit}} = \lceil \texttt{n\_add} \cdot \texttt{decay}^{\,r-1} \rceil$ draws with highest novelty that exceed the threshold, and append them to the reservoir: $S_r = S_{r-1} \cup {\text{selected draws}}$.
    • Rebuild the GP metric on the enlarged reservoir $S_r$.

  3. Final sampling
    Using the final GP metric, draw $\texttt{n\_iter}$ samples after $\texttt{n\_warmup}$ warmup iterations via HMC.

Returns
An object of class c("gdpar_geom_rmhmc_adaptive", "list") with the following components:

  • draws : numeric matrix of posterior draws (rows = iterations, columns = parameters).
  • metric : the final GP metric object.
  • reservoir : numeric matrix of all parameter sites collected in the reservoir.
  • n_sites_trace : integer vector recording the reservoir size after each adaptation round.
  • novelty_trace : numeric vector of the maximum novelty observed in each adaptation round.
  • accept_rate : acceptance rate of the final HMC run.
  • ebfmi : empirical Bayes fraction of missing information (E‑BFMI) of the final run.
  • n_divergent : number of divergent transitions in the final run.
  • n_rounds : number of adaptation rounds actually performed.
  • epsilon : step size used.
  • L : integer, number of leapfrog steps used.

Notes

  • If target is not already a "gdpar_geom_target", it is coerced via gdpar_geom_target().
  • The function aborts with class "gdpar_input_error" if n_sites_init < 2 or if fisher is neither a function nor NULL.
  • When seed is provided, the function temporarily modifies .Random.seed and restores it on exit.
  • If warmup_metric is NULL, a Riemannian metric with softabs curvature is built using alpha.
  • The GP metric is built by gdpar_geom_metric_gp_fisher; its novelty method is used to assess new points.
  • The decay formula for the number of admitted sites is $n_{\text{admit}} = \lceil \texttt{n\_add} \cdot \texttt{decay}^{\,r-1} \rceil$.

print.gdpar_geom_rmhmc_adaptive(x, ...)

Purpose
S3 print method for objects of class "gdpar_geom_rmhmc_adaptive". Prints a concise summary of the sampling results.

Arguments

  • x : an object of class "gdpar_geom_rmhmc_adaptive".
  • ... : additional arguments (currently unused).

Returns
Invisibly returns x.

Notes

  • The output includes: number of draws, dimensionality, number of adaptation rounds, final reservoir size, tuning parameters (epsilon, L), acceptance rate, number of divergences, E‑BFMI, and the final novelty value (if available).

R/geometry_laplace.R

Internal helpers (promoted extractor core)

.gdpar_geom_newton_climb(geom_target, ref, n_steps = 300L, ridge = 1e-6, tol = 1e-3, max_step_norm = 5)

Purpose

Modified-Newton ascent to the posterior mode, used as the polish stage of .gdpar_geom_laplace_climb. It is the in-package promotion of the benchmark extractor routine rg7_climb_to_mode. It reads the same log_prob, grad_log_prob, and hessian that the sampler uses, so the mode it finds is the mode of the same target the orchestrator certified.

Arguments

Argument Type Meaning
geom_target list (a gdpar_geom_target) Must expose $log_prob, $grad_log_prob, $hessian (a function; the caller checks first).
ref numeric vector Starting position for the climb.
n_steps integer (scalar) Maximum number of Newton iterations (default 300L).
ridge numeric (scalar) Relative eigenvalue floor: each eigenvalue of $A=-H$ is raised to at least `ridge * max(
tol numeric (scalar) Gradient-norm convergence tolerance (default 1e-3).
max_step_norm numeric (scalar) Trust-region radius; any step with $\lVert\text{step}\rVert_2 &gt;$ max_step_norm is rescaled (default 5).

Mathematics

At each iteration $k$ with current point $\theta_k$ and $f_0 = \log p(\theta_k)$:

  1. Evaluate $g = \nabla\log p(\theta_k)$ and $H = \nabla^2\log p(\theta_k)$.
  2. Symmetrise: $H \leftarrow \tfrac12(H + H^\top)$.
  3. Eigen-decompose $A = -H = V\Lambda V^\top$ (symmetric).
  4. Floor the spectrum: $$\tilde\lambda_i = \max!\bigl(\lambda_i,;\texttt{ridge}\cdot\max(|\lambda_i|,1)\bigr).$$
  5. Form the ascent direction $$\text{step} = V,\bigl(\tilde\Lambda^{-1}(V^\top g)\bigr) = \tilde A^{-1}g,$$ which is $A^{-1}g$ on the floored spectrum.
  6. Trust-region cap: if $\lVert\text{step}\rVert_2 &gt; \texttt{max\_step\_norm}$, rescale $\text{step}\leftarrow\text{step}\cdot\texttt{max\_step\_norm}/\lVert\text{step}\rVert_2$.
  7. Backtracking Armijo line search over $t\in{1,\tfrac12,\tfrac14,\dots}$ (up to 20 halvings): accept $\theta_k + t,\text{step}$ when $$\log p(\theta_k + t,\text{step}) ;\ge; f_0 + 10^{-4},t,g^\top\text{step}.$$
  8. Convergence test: $\lVert g\rVert_2 &lt; \texttt{tol}$ or improvement $f_{\text{new}}-f_0 &lt; \texttt{tol}\cdot 10^{-2}$.

When $H$ contains non-finite entries (a stiff off-mode region), the code falls back to a steepest-ascent step $\text{step}=g$ so the climb can traverse the region; eigen() is never fed a non-finite matrix.

Returns

A list:

Element Type Meaning
mode numeric vector Best point found.
converged logical TRUE if the gradient/improvement test fired.
iters integer Iteration count at convergence (present only when converged = TRUE).
grad_norm numeric $\lVert\nabla\log p(\text{mode})\rVert_2$.
logp numeric $\log p(\text{mode})$.
reason character "no exact Hessian" (present only on the early-return path).

Notes

  • Early return with converged = FALSE and reason = "no exact Hessian" when geom_target$hessian is not a function.
  • The loop breaks (returning converged = FALSE) when the gradient is non-finite, the step is non-finite, or the line search exhausts all 20 backtracks without acceptance.
  • No side effects; no RNG usage; no S3 dispatch.

.gdpar_geom_laplace_climb(geom_target, start, climb_steps)

Purpose

Two-stage mode-finder: an L-BFGS-B warm start on $-\log p$ (robust when start is far from the mode) followed by a modified-Newton polish when an exact Hessian is available; L-BFGS-B alone otherwise. Both stages read the same geom_target target and gradient the sampler uses.

Arguments

Argument Type Meaning
geom_target list (gdpar_geom_target) Exposes $log_prob, $grad_log_prob, optionally $hessian.
start numeric vector Warm-start position (the orchestrator passes the best on-ridge position found during sampling).
climb_steps integer (scalar) Maximum Newton steps passed as n_steps to .gdpar_geom_newton_climb.

Mathematics

Stage 1 — L-BFGS-B on the negated log-density:

$$\tilde f(\theta) = -\log p(\theta),\qquad \nabla\tilde f(\theta) = -\nabla\log p(\theta),$$

with non-finite values clamped: $\tilde f\to 10^{18}$ on error, and non-finite gradient entries set to $0$. stats::optim is called with method = "L-BFGS-B", maxit = 5000L, factr = 1e1.

Stage 2 — if geom_target$hessian is a function, the L-BFGS-B solution warm is passed to .gdpar_geom_newton_climb(geom_target, warm, n_steps = climb_steps) for Newton polishing. Otherwise the warm start is returned directly.

Returns

A list:

Element Type Meaning
mode numeric vector Best point found.
converged logical From Newton polish (converged field) when Hessian available; otherwise isTRUE(opt$convergence == 0L).
grad_norm numeric $\lVert\nabla\log p(\text{mode})\rVert_2$.
logp numeric $\log p(\text{mode})$.
iters integer Present only when the Newton polish converged (inherited from .gdpar_geom_newton_climb).

Notes

  • If stats::optim errors, opt is NULL and warm falls back to start.
  • All tryCatch blocks swallow errors silently, returning sentinel values (Inf, NULL, zero vectors).
  • No side effects beyond the optimisation; no S3 dispatch.

.gdpar_geom_observed_information(geom_target, reference, h = 1e-4)

Purpose

Compute the observed information matrix $M = -\nabla^2\log p$ at reference, symmetrised. Promotion of rg7_observed_information. Prefers the exact Hessian; falls back to central finite differences of the gradient when no Hessian is exposed or it is non-finite. The route taken is recorded in attr(., "method").

Arguments

Argument Type Meaning
geom_target list (gdpar_geom_target) Exposes $dim, $grad_log_prob, optionally $hessian.
reference numeric vector Point at which to evaluate the curvature.
h numeric (scalar) Central finite-difference step (default 1e-4).

Mathematics

Exact-Hessian route (when geom_target$hessian is a function and returns a finite matrix):

$$M = -\tfrac12\bigl(H + H^\top\bigr),\qquad H = \texttt{geom_targethessian}(\texttt{reference}).$$

Finite-difference route (fallback): the Jacobian of the gradient is approximated column-wise by central differences,

$$J_{ij} = \frac{\partial g_i}{\partial\theta_j}\Big|_{\text{ref}} \approx \frac{g_i(\text{ref}+h,e_j) - g_i(\text{ref}-h,e_j)}{2h},$$

and then

$$M = -\tfrac12\bigl(J + J^\top\bigr).$$

Returns

A $d\times d$ numeric matrix with attribute "method" set to one of:

method value Condition
"exact_hessian" Hessian function present and all entries finite.
"finite_difference" No usable Hessian; finite-difference matrix all finite.
"finite_difference_nonfinite" No usable Hessian; finite-difference matrix contains non-finite entries.

Notes

  • The exact-Hessian route is tried inside a tryCatch; any error yields NULL and triggers the fallback.
  • Finite differences loop over seq_len(d) columns, evaluating the gradient at reference ± h*e_j.
  • No side effects; no S3 dispatch.

.gdpar_geom_laplace_draws_unconstrained(mode, Lhalf, S, seed)

Purpose

Generate $S$ iid Laplace draws on the unconstrained scale: $\theta_s = \hat\theta + M^{-1/2}z_s$ with $\mathrm{Cov}(\theta_s)=M^{-1}$. Uses a local RNG that saves and restores the global .Random.seed, so the caller's RNG stream is untouched.

Arguments

Argument Type Meaning
mode numeric vector ($\hat\theta$, length $d$) The posterior mode.
Lhalf $d\times d$ numeric matrix The symmetric matrix square root $M^{-1/2}$ (i.e. $\texttt{Lhalf}\cdot\texttt{Lhalf}=M^{-1}$).
S integer (scalar) Number of draws.
seed integer (scalar) Seed passed to set.seed().

Mathematics

Draw $Z\in\mathbb{R}^{S\times d}$ with i.i.d. $N(0,1)$ entries, then

$$U = Z,M^{-1/2},\qquad \theta_s = \hat\theta + U_{s,\cdot}.$$

Because $M^{-1/2}$ is symmetric, $\mathrm{Cov}(U_{s,\cdot}) = M^{-1/2}(M^{-1/2})^\top = M^{-1}$, as required. The addition of mode is performed via sweep(U, 2L, mode, "+").

Returns

An $S\times d$ numeric matrix of unconstrained Laplace draws.

Notes

  • RNG isolation: old_seed is captured from globalenv() if .Random.seed exists; on.exit restores it (or does nothing if it did not exist). set.seed(seed) is called after the save.
  • No S3 dispatch; no error handling on Lhalf dimensions — the caller is responsible for conformability.

.gdpar_geom_laplace_fit_quality(geom_target, mode, M, logdetM, U)

Purpose

Assess the fidelity of the Laplace Gaussian $q$ against the true posterior $p$ over the same draws U. Promotion of rg7_laplace_fit_quality. Computes the self-normalised importance-sampling ESS, the PSIS Pareto-$k$ of the weights $p/q$, and the mean/max log-density drop $\log p(\hat\theta)-\log p(\theta)$ against its Gaussian expectation $d/2$.

Arguments

Argument Type Meaning
geom_target list (gdpar_geom_target) Exposes $log_prob.
mode numeric vector ($\hat\theta$) The mode.
M $d\times d$ matrix The (floored, positive-definite) precision matrix.
logdetM numeric (scalar) $\log\det M$.
U $S\times d$ matrix Draws at which to evaluate (the same draws returned to the caller or internally generated for diagnostics).

Mathematics

The Laplace Gaussian log-density at each draw $\theta_s = U[s,]$:

$$\log q(\theta_s) = -\tfrac{d}{2}\log(2\pi) + \tfrac12\log\det M - \tfrac12(\theta_s-\hat\theta)^\top M,(\theta_s-\hat\theta).$$

The quadratic form is computed as rowSums((cen %*% M) * cen) where cen = sweep(U, 2, mode, "-").

The true log-density $\log p(\theta_s)$ is evaluated one draw at a time via geom_target$log_prob, with errors caught and mapped to NA_real_.

Importance weights (log-scale, max-stabilised):

$$\ell_s = \log p(\theta_s) - \log q(\theta_s),\qquad \ell_s \leftarrow \ell_s - \max_s \ell_s,\qquad w_s = e^{\ell_s}.$$

Self-normalised IS effective sample size:

$$\text{ESS}_{\text{IS}} = \frac{(\sum_s w_s)^2}{\sum_s w_s^2}.$$

PSIS Pareto-$k$: when the loo package is available and there are more than 10 finite draws, loo::psis() is called on the log-weights $\log p - \log q$ (as a one-column matrix, r_eff = NA), and $diagnostics$pareto_k is extracted. Otherwise pareto_k = NA_real_.

Log-density drop:

$$\Delta_s = \log p(\hat\theta) - \log p(\theta_s),\qquad \text{logdrop_mean} = \overline{\Delta},\quad \text{logdrop_max} = \max_s \Delta_s.$$

Under a perfect Gaussian, $\Delta_s \sim \chi^2_d/2$ with mean $d/2$.

Returns

A list:

Element Type Meaning
ess_is numeric Self-normalised IS ESS (raw count).
ess_is_frac numeric ess_is / sum(fin) (fraction of finite draws).
pareto_k numeric PSIS Pareto-$k$, or NA if unavailable.
logdrop_mean numeric Mean of $\Delta_s$ over finite draws.
logdrop_max numeric Max of $\Delta_s$ over finite draws.
logdrop_expected numeric Always $d/2$.
n_finite integer Count of draws with both $\log p$ and $\log q$ finite.
n_total integer S (total draws).

When sum(fin) < 2L, all diagnostic fields except logdrop_expected, n_finite, and n_total are NA_real_.

Notes

  • logp_mode is evaluated separately (with error handling) for the drop calculation.
  • The loo::psis call is wrapped in suppressWarnings() and a tryCatch; any error yields NA_real_.
  • No side effects; no S3 dispatch.

.gdpar_geom_laplace_label(fq, all_pos, d, mode_offset_sd = 0)

Purpose

Distil the raw fidelity diagnostics from .gdpar_geom_laplace_fit_quality into a single, unambiguous scalar label ("good" / "poor" / "very_poor") so a caller can never mis-read raw numbers. Implements the data-driven thresholds adjudicated in design session B9.34.

Arguments

Argument Type Meaning
fq list The fidelity-diagnostics list returned by .gdpar_geom_laplace_fit_quality.
all_pos logical Was the raw (un-floored) curvature positive-definite?
d integer (scalar) Unconstrained dimension.
mode_offset_sd numeric (scalar) Newton-decrement bound on the mode offset in posterior SDs (default 0).

Mathematics / decision logic

Let $\text{ess} = \texttt{fq\$ess\_is\_frac}$, $\text{ld} = \texttt{fq\$logdrop\_mean}$, $\text{pk} = \texttt{fq\$pareto\_k}$, and $\text{off} = \texttt{mode\_offset\_sd}$ (or $\infty$ if non-finite).

  1. Non-PD curvature (saddle, not a maximum): if !all_pos $\Rightarrow$ "very_poor".
  2. Non-finite diagnostics: if ess or ld is non-finite $\Rightarrow$ "very_poor".
  3. Mode-location gate: if $\text{off} &gt; 2$ $\Rightarrow$ "very_poor". Define well_centred $= (\text{off} \le 0.5)$.
  4. Excellent fit override: if well_centred and $\text{ess}\ge 0.9$ and $\text{ld}\le d$ $\Rightarrow$ "good" (Pareto-$k$ on numerically constant weights is uninformative and must not override).
  5. Severe degradation: define pk_bad $= (\text{pk}\ge 1.0)$. If $\text{ess}&lt;0.1$ or $\text{ld}\ge 2d$ or (pk_bad and $\text{ess}&lt;0.5$) $\Rightarrow$ "very_poor".
  6. Middle ground: define pk_ok $= (\text{is.na(pk)} \text{ or } \text{pk}&lt;0.7)$. If well_centred and $\text{ess}\ge 0.5$ and $\text{ld}\le d$ and pk_ok $\Rightarrow$ "good".
  7. Otherwise: "poor".

Returns

A character scalar: "good", "poor", or "very_poor".

Notes

  • The ESS fraction is the stable primary signal (no tail fit); Pareto-$k$ only refines the middle ground and is damning only when ESS is also weak.
  • The mean log-density drop threshold uses $d$ (benign $\le d$) and $2d$ (severe $\ge 2d$) relative to the Gaussian expectation $d/2$.
  • A non-positive-definite curvature at the "mode" returns "very_poor" outright because the Laplace Gaussian is ill-defined there.
  • No side effects; no S3 dispatch.

gdpar_geom_laplace(geom_target, reference = ..., draws = 0L, climb = TRUE, seed = ..., fit_quality_draws = 256L, eigen_floor_rel = 1e-10, climb_steps = 300L, cond_warn = 1e12)

Note: The roxygen documentation for this exported function appears in this section, but the function body (gdpar_geom_laplace <- function(...) { ... }) is in section 2 of 2. The subsection below documents the interface and contract as declared by the roxygen; the implementation details will be covered when section 2 is processed.

Purpose

The exported, first-class Laplace approximation of a posterior at its mode. Compute the mode + curvature Gaussian $N(\hat\theta, M^{-1})$ on the unconstrained scale, report its covariance, optional draws, and a fidelity diagnostic of the Gaussian against the true posterior. This is the honest endpoint for a posterior the geometry-adaptive orchestrator certifies as a non-sampleable, genuinely non-Gaussian canyon (Block RG, RG.7) — the regime of mgcv/REML and INLA/Laplace, accompanied by a measurement of how good the Gaussian is.

Arguments

Argument Type Meaning
geom_target gdpar_geom_target (or coercible) Exposes log_prob, grad_log_prob, dim, and ideally hessian.
reference numeric vector Optional unconstrained warm-start position for the mode climb; defaults to the origin. When called by the orchestrator this is the best on-ridge position found during sampling.
draws integer (scalar) Number of iid Laplace draws $\hat\theta + M^{-1/2}z$ to return (default 0L: the mode plus precision Gaussian is the approximation).
climb logical Whether to climb to the mode from reference (default TRUE). FALSE treats reference as the mode and only reads the curvature there.
seed integer Seed for the local, stream-preserving draw RNG; defaults to a fixed value for reproducibility.
fit_quality_draws integer (scalar) Number of internal iid draws used to assess fidelity when draws is small (default 256L); the fidelity label is always computed independently of user-facing draws.
eigen_floor_rel numeric (scalar) Relative eigenvalue floor applied to $M$ for positive-definiteness (default 1e-10 of the largest eigenvalue).
climb_steps integer (scalar) Maximum modified-Newton steps in the mode climb (default 300L).
cond_warn numeric (scalar) Condition-number threshold above which the un-floored curvature triggers an ill-conditioning warning (default 1e12).

Mathematics

The Laplace approximation is

$$q(\theta) = N!\bigl(\hat\theta,; M^{-1}\bigr),\qquad M = -\nabla^2\log p(\hat\theta),$$

where $\hat\theta$ is the posterior mode. The mode is reached by an L-BFGS-B warm start on $-\log p$ followed, when the target exposes an exact Hessian, by a modified-Newton polish (.gdpar_geom_laplace_climb → .gdpar_geom_newton_climb). The precision $M$ is the symmetrised $-$Hessian (.gdpar_geom_observed_information), eigen-floored to be positive-definite. Draws are $\theta_s = \hat\theta + M^{-1/2}z_s$ (.gdpar_geom_laplace_draws_unconstrained). Fidelity is scored via IS-ESS, PSIS Pareto-$k$, and the log-density drop against $d/2$ (.gdpar_geom_laplace_fit_quality), then distilled into a label (.gdpar_geom_laplace_label).

Returns

An object of class gdpar_geom_laplace: a list with:

Element Type Meaning
mode numeric vector $\hat\theta$.
M $d\times d$ matrix Floored precision.
cov $d\times d$ matrix $M^{-1}$.
Lhalf $d\times d$ matrix $M^{-1/2}$ (symmetric square root).
logdet numeric $\log\det M$.
eig numeric vector Eigenvalues of $M$.
cond numeric Floored condition number.
cond_unfloored numeric Un-floored condition number.
n_floored integer Count of eigen-floored directions.
floor_value numeric The floor value used.
all_pos logical Was the raw curvature positive-definite?
mode_offset_sd numeric Newton-decrement bound on the mode offset in posterior SDs.
grad_norm numeric $\lVert\nabla\log p(\hat\theta)\rVert_2$.
logp numeric $\log p(\hat\theta)$.
converged logical Did the mode climb converge?
method character Hessian route ("exact_hessian", "finite_difference", etc.).
dim integer $d$.
draws $S\times d$ matrix Laplace draws (possibly zero rows), with attr(., "approximation") = "laplace".
fit_quality list Fidelity diagnostics from .gdpar_geom_laplace_fit_quality.
fit_quality_label character Scalar label: "good" / "poor" / "very_poor".

Notes

  • The Laplace approximation is never advertised as exact; the draws matrix carries attr(., "approximation") = "laplace" so it is never mistaken for exact MCMC downstream.
  • A warning is raised when the un-floored condition number exceeds cond_warn, indicating the eigen-floor could be masking ill-conditioning.
  • A non-positive-definite raw curvature (a saddle) is flagged loudly (the label is "very_poor").
  • The default fit path of the package is byte-for-byte unchanged when laplace_fallback = FALSE (the default); this function is only wired in behind the opt-in flag.
  • S3 class: the return value has class gdpar_geom_laplace.

gdpar_geom_laplace(geom_target, reference = NULL, draws = 0L, climb = TRUE, seed = NULL, fit_quality_draws = 256L, eigen_floor_rel = 1e-10, climb_steps = 300L, cond_warn = 1e12)

Purpose Constructor for the Laplace (Gaussian) approximation of a target density exposed through a gdpar_geom_target object. It locates the posterior mode, builds the precision matrix from the observed information, eigen-floors it to guarantee positive-definiteness, draws an internal sample for fidelity diagnostics, and returns a tagged gdpar_geom_laplace object consumed by downstream gdpar reference-anchoring routines.

Arguments

  • geom_target: object of class "gdpar_geom_target" (or coercible via gdpar_geom_target()). Must expose $dim, $log_prob(x), $grad_log_prob(x), and whatever .gdpar_geom_observed_information() requires.
  • reference: numeric vector of length d (the unconstrained dimension). Starting point for mode climbing when climb = TRUE, or used directly as the mode when climb = FALSE. Defaults to rep(0, d).
  • draws: non-negative integer. Number of Laplace draws to attach to the returned object. Coerced via as.integer.
  • climb: logical. If TRUE, the mode is sought with .gdpar_geom_laplace_climb(geom_target, reference, climb_steps); if FALSE, reference is treated as the mode and its gradient norm / log-density are evaluated directly.
  • seed: integer count seed for the internal RNG used by .gdpar_geom_laplace_draws_unconstrained(). Defaults to 20260625L when NULL; validated by assert_count(seed, "seed").
  • fit_quality_draws: integer lower-bounding the internal sample size used to compute fit quality. Coerced to max(as.integer(.), 2L). Defaults to 256L.
  • eigen_floor_rel: positive numeric. Relative eigenvalue floor: $\lambda_{\text{floor}} = \lambda_{\max}\cdot\text{eigen\_floor\_rel}$. Defaults to 1e-10.
  • climb_steps: integer passed to .gdpar_geom_laplace_climb(). Defaults to 300L.
  • cond_warn: numeric threshold on the un-floored curvature condition number above which a warning is emitted. Set to a non-finite value to suppress. Defaults to 1e12.

Mathematics

Mode: $$ \hat\mu = \begin{cases} \text{climb}(\text{geom_target}, \text{reference}, \text{climb_steps}) & \text{if } \text{climb} = \text{TRUE},\ \text{reference} & \text{otherwise.} \end{cases} $$

Precision (observed information), symmetrized: $$ M = -\nabla^2 \log p(\hat\mu), \qquad M \leftarrow \tfrac{1}{2}(M + M^\top). $$

Symmetric eigendecomposition $M = V\Lambda V^\top$ with $\Lambda = \mathrm{diag}(\lambda_1,\dots,\lambda_d)$. Eigen-floor: $$ \tilde\lambda_i = \max(\lambda_i,; \lambda_{\max}\cdot\text{eigen_floor_rel}), \qquad \lambda_{\max} = \max_i \lambda_i. $$

Covariance, symmetric inverse square root, and log-determinant (all using floored eigenvalues): $$ M^{-1} = V,\mathrm{diag}(\tilde\lambda_i^{-1}),V^\top, \qquad M^{-1/2} = V,\mathrm{diag}(\tilde\lambda_i^{-1/2}),V^\top, \qquad \log\det M = \sum_i \log \tilde\lambda_i. $$

Newton step and Newton decrement (mode offset in SD units): $$ \Delta = M^{-1},\nabla\log p(\hat\mu), \qquad \sigma_{\text{offset}} = \sqrt{\max!\bigl((\nabla\log p)^\top M^{-1},\nabla\log p,;0\bigr)}. $$

Condition numbers: $$ \kappa_{\text{floored}} = \frac{\max_i \tilde\lambda_i}{\min_i \tilde\lambda_i}, \qquad \kappa_{\text{unfloored}} = \begin{cases} \dfrac{\max_i |\lambda_i|}{\min_i |\lambda_i|} & \min_i |\lambda_i| > 0,\ +\infty & \text{otherwise.} \end{cases} $$

Internal sample size $S_{\text{fit}} = \max(\text{draws},\,\text{fit\_quality\_draws})$; unconstrained draws $U_{\text{full}}$ are produced by .gdpar_geom_laplace_draws_unconstrained(mode, Lhalf, S_fit, seed), then fed to .gdpar_geom_laplace_fit_quality(geom_target, mode, M, logdetM, U_full). The fidelity label is .gdpar_geom_laplace_label(fq, all_pos, d, mode_offset_sd).

Returns

A list with class c("gdpar_geom_laplace", "list") and elements:

  • mode: numeric vector of length d, the located mode.
  • M: the symmetrized (but not eigen-floored) precision matrix $-\nabla^2\log p(\hat\mu)$; may carry a "method" attribute set by .gdpar_geom_observed_information().
  • cov: $M^{-1}$ computed with floored eigenvalues.
  • Lhalf: symmetric $M^{-1/2}$ computed with floored eigenvalues.
  • logdet: $\log\det M$ using floored eigenvalues.
  • eig: numeric vector of unfloored eigenvalues of $M$.
  • cond: floored condition number $\kappa_{\text{floored}}$.
  • cond_unfloored: $\kappa_{\text{unfloored}}$.
  • n_floored: integer count of eigenvalues strictly below floor_value.
  • floor_value: $\lambda_{\max}\cdot\text{eigen\_floor\_rel}$.
  • all_pos: logical, TRUE iff every unfloored eigenvalue is $&gt;0$.
  • mode_offset_sd: Newton decrement $\sigma_{\text{offset}}$.
  • grad_norm: gradient $L_2$ norm at the mode (from cl$grad_norm, or NA_real_ if missing).
  • logp: as.numeric(cl$logp %||% geom_target$log_prob(mode)).
  • converged: isTRUE(cl$converged).
  • method: attr(M, "method") %||% "unknown".
  • dim: integer dimension d.
  • draws: $S_{\text{returned}}\times d$ numeric matrix, where $S_{\text{returned}}=\text{draws}$ (a 0-row matrix when draws == 0L), carrying attribute "approximation" = "laplace".
  • fit_quality: the object returned by .gdpar_geom_laplace_fit_quality().
  • fit_quality_label: the label string returned by .gdpar_geom_laplace_label().

Notes

  • Input errors (class "gdpar_input_error", raised via gdpar_abort): draws is NA or negative; reference (when non-NULL) has length $\neq d$.
  • Geometry warnings (class "gdpar_geometry_warning", raised via gdpar_warn): (i) any unfloored eigenvalue $\le 0$ (reports the count and minimum eigenvalue, flags the approximation as "very_poor" via the label machinery); (ii) cond_unfloored > cond_warn when cond_warn is finite (reports both numbers, cautions that the eigen-floor may mask ill-conditioning).
  • fit_quality_draws is silently raised to at least 2L.
  • The stored M is symmetrized but not floored; only cov, Lhalf, logdet, and cond use the floored spectrum. eig reports the unfloored spectrum so callers can detect silent flooring via n_floored and floor_value.
  • The returned draws matrix is tagged with attr(., "approximation") <- "laplace" so downstream methods cannot silently consume it as exact MCMC output; the attribute survives subsetting because it is set on the final sliced matrix.
  • When climb = FALSE, converged is NA (propagated through isTRUE → FALSE), and grad_norm is computed explicitly as $\sqrt{\sum_i (\nabla\log p(\text{reference})_i)^2}$.
  • logp falls back to geom_target$log_prob(mode) when cl$logp is NULL.
  • seed = NULL is replaced by the fixed default 20260625L before validation.

print.gdpar_geom_laplace(x, ...)

Purpose S3 print method for objects of class "gdpar_geom_laplace". Renders a compact multi-line diagnostic summary covering dimension, fit-quality label, mode optimality, curvature conditioning, importance-sampling fidelity, and draw dimensions.

Arguments

  • x: a gdpar_geom_laplace object (as produced by gdpar_geom_laplace()).
  • ...: unused; absorbed for S3 generic conformance.

Mathematics

The "unreliable" flag on the fidelity line is raised when $$ \hat k_{\text{Pareto}} \ge 0.7, $$ using x$fit_quality$pareto_k. The expected log-drop reference shown is x$fit_quality$logdrop_expected, printed as E[d/2] (the theoretical Laplace expectation $\mathbb{E}[d/2]$ for a $d$-dimensional Gaussian).

Returns

Invisibly returns x (via invisible(x)).

Notes

  • Exported S3 method registered for class "gdpar_geom_laplace".
  • Output layout (in order):
    1. <gdpar_geom_laplace> dim {dim} | fit-quality: {fit_quality_label}
    2. mode: |grad| {grad_norm} | offset {mode_offset_sd} SD | PD {all_pos} | method {method}
    3. condition: floored {cond} | un-floored {cond_unfloored} | eigen-floored dirs {n_floored}/{dim}
    4. fidelity: IS-ESS {ess_is} ({100*ess_is_frac}%) [unreliable: Pareto-k >= 0.7] | Pareto-k {pareto_k} | log-drop {logdrop_mean} (E[d/2] = {logdrop_expected}) — the bracketed [unreliable: ...] segment is emitted only when isTRUE(fq$pareto_k >= 0.7).
    5. If fit_quality_label is not "good": NOTE: a '{label}' Laplace approximation -- a plug-in / INLA-style fallback, NOT exact MCMC.
    6. draws: {nrow(draws)} x {ncol(draws)}
  • Numeric formatting uses format(..., digits = 3) (or digits = 2 for the ESS percentage); integers (dim, n_floored, method) are printed raw via sep = "".
  • No side effects beyond console output; does not mutate x.

R/geometry_orchestrator.R

gdpar_geom_orchestrate_criteria()

Purpose Defines the multi-signal success criteria for the geometry-adaptive orchestrator. It establishes the conjunction of size-invariant sampler signals (acceptance rate, divergence rate, and E-BFMI) that a sampler level must meet to be considered successful, preventing false positives from a sampler that moves too little.

Arguments None.

Mathematics Implements a success gate requiring:

  • Acceptance rate in a healthy band: $0.5 \leq \hat{\alpha} \leq 0.999$
  • Divergence rate below a ceiling: $\hat{\delta} \leq 0.02$
  • Energy Bayesian Fraction of Missing Information (E-BFMI) above a floor: $\widehat{\text{E-BFMI}} \geq 0.3$

Returns A named list with four numeric elements:

  • accept_low (0.5): Lower bound of the healthy acceptance band.
  • accept_high (0.999): Upper bound of the healthy acceptance band.
  • divergent_rate_high (0.02): Maximum allowed divergence rate.
  • ebfmi_low (0.3): Minimum allowed E-BFMI.

Notes These are static, calibrated defaults. The decisive gate used in gdpar_geom_orchestrate applies a hysteresis margin (tightening these thresholds by 1 + hysteresis).

gdpar_geom_orchestrate_budget()

Purpose Defines the budget, stopping rules, and tuning knobs that bound the closed-loop operation of the geometry-adaptive orchestrator. It controls the maximum computational effort (rounds, fits, time) and the parameters for the successive-halving budget tiering and step-size tuning.

Arguments None.

Mathematics Sets the parameters for:

  • A cheap probe tier: $\text{probe\_warmup}=150, \text{probe\_iter}=150$
  • A full tier: $\text{full\_warmup}=500, \text{full\_iter}=500$
  • Base HMC integrator: $\epsilon = 0.25, L=25$
  • Tuning: if tune_epsilon=TRUE, a short search over $\epsilon$ is run with tune_iter=60 iterations per probe.

Returns A named list with 15 elements:

  • max_rounds (8L): Hard ceiling on diagnose–sample cycles.
  • max_levels (5L): Cap on distinct ladder levels tried.
  • probe_warmup (150L), probe_iter (150L): Budget for cheap probes.
  • full_warmup (500L), full_iter (500L): Budget for full runs.
  • epsilon (0.25), L (25L): Base leapfrog step size and trajectory length.
  • tune_epsilon (TRUE): Whether to run a coarse $\epsilon$ search.
  • tune_iter (60L): Iterations per tuning probe.
  • max_seconds (Inf), max_seconds_per_fit (Inf): Wall-time caps.
  • max_fits (40L): Global fit cap.
  • n_rediagnose (1L): Number of pilots per re-diagnosis.
  • stall_limit (2L): Consecutive no-progress rounds before stopping.
  • hysteresis (0.1): Margin for tightening the decisive success gate.

Notes The loop checks remaining budget before starting each fit and emits a certified limit if caps are exceeded. The n_rediagnose ≥ 2 enables a stability check on the re-diagnosis.

.gdpar_orch_ladder()

Purpose Defines the ordered hierarchy of sampler levels (the "ladder") that the orchestrator can select and escalate through. Each level is a geometry-aware HMC variant.

Arguments None.

Mathematics The hierarchy levels and their indices: 0. Euclidean diagonal

  1. Euclidean dense
  2. Riemannian
  3. Relativistic
  4. Sub-Riemannian

Note: Indices are not contiguous; levels 6 (tempering/reparametrisation) and -1 (reparametrise/eliminate) are deliberately excluded.

Returns A named list of five elements, each a sublist with:

  • key: Character string identifier.
  • index: Integer for ordering.
  • needs_fisher: Logical, TRUE only for sub_riemannian.

Notes Levels 6 and -1 are out-of-scope for the orchestrator; a diagnosis pointing to them short-circuits to a certificate.

.gdpar_orch_ladder_keys()

Purpose Returns the character vector of in-ladder level keys in ascending hierarchical order.

Arguments None.

Mathematics None.

Returns A character vector: c("euclidean_diagonal", "euclidean_dense", "riemannian", "relativistic", "sub_riemannian").

Notes Pure utility function; used for iteration over the ladder.

.gdpar_orch_pathology_to_key(pathology)

Purpose Maps a classified pathology (as returned by gdpar_geometry_diagnostic) to the corresponding ladder key, or NA for out-of-scope pathologies.

Arguments

  • pathology (character): The pathology classification string (e.g., "isotropic", "funnel").

Mathematics None; a direct lookup table.

Returns A character string (ladder key) or NA_character_ for:

  • "multimodal", "boundary", "flat_direction", "unknown".

Notes This is a one-to-one mapping except for the out-of-scope pathologies, which yield NA.

.gdpar_orch_next_above(index)

Purpose Finds the first ladder key strictly above a given hierarchy index. Used for monotone escalation.

Arguments

  • index (integer): The current ladder index.

Mathematics Given an index $i$, returns the key with the smallest index $j$ such that $j &gt; i$.

Returns A character string (ladder key) or NA_character_ if the input index is at or above the top of the ladder.

Notes Assumes the ladder defined by .gdpar_orch_ladder() is non-empty and indexed.

.gdpar_orch_pathology_scores(sig, grows, th)

Purpose Computes a continuous, size-invariant proximity score for each in-ladder pathology. This is the second, independent estimator in the combined level selection map, cross-checked against the discrete classifier.

Arguments

  • sig (data.frame row): A single row of size-invariant diagnostic signals (columns include condition_number, heavy_kurtosis, ebfmi_min).
  • grows (logical): Whether the difficulty curve indicates the problem grows with sample size (from diag$difficulty_curve$grows_with_n).
  • th (list): Thresholds from the diagnostic classifier (must include condition_high, heavy_kurtosis_high, funnel_ebfmi_low, heavy_cond_max).

Mathematics Computes five scores (one per pathology) from the signals, scaled by the classifier's thresholds:

Let:

  • $c = \text{sig\$condition\_number}$ (default 1 if non-finite)
  • $k = \text{sig\$heavy\_kurtosis}$ (default 0 if non-finite)
  • $e = \max(\text{sig\$ebfmi\_min}, 10^{-3})$

Ratios:

  • $r_c = c / \text{th\$condition\_high}$
  • $r_k = k / \text{th\$heavy\_kurtosis\_high}$

Scores:

  • Isotropic: $- \max(r_c, r_k)$
  • Anisotropic: $r_c \times \begin{cases} 0 & \text{if } \text{grows} \ 1 & \text{otherwise} \end{cases} \times \begin{cases} 1 & \text{if } r_k < 1 \ 0.25 & \text{otherwise} \end{cases}$
  • Funnel: $r_k \times (\text{th\$funnel\_ebfmi\_low} / e)$
  • Heavy tails: $r_k \times (\text{th\$heavy\_cond\_max} / \max(c, 1)) \times \min(1, e / \text{th\$funnel\_ebfmi\_low})$
  • Quasi-deterministic: $r_c \times \begin{cases} 2 & \text{if } \text{grows} \ 0.5 & \text{otherwise} \end{cases}$

Returns A named numeric vector of length 5 with the scores for each pathology: isotropic, anisotropic, funnel, heavy_tails, quasi_deterministic.

Notes The design ensures that a deep funnel (low E-BFMI) is not mistaken for heavy tails by damping the heavy-tail score when E-BFMI is depressed. All signals are made safe for non-finite values.

.gdpar_orch_select_entry(diag, level_map, entry_level, th)

Purpose Resolves the entry level for the orchestrator using the combined map: a transparent, calibrated rule-based classifier (primary) and a continuous proximity score (secondary). Implements the conflict-resolution strategy: on disagreement or low confidence, start at the lower of the two candidate ladder levels.

Arguments

  • diag (list): A gdpar_geometry_diagnostic object (or compatible list) with fields $pathology, $signals, $difficulty_curve$grows_with_n, $confidence.
  • level_map (list or NULL): User override mapping from pathology to ladder key.
  • entry_level (character or NULL): Direct user override for the starting level key.
  • th (list): Thresholds for the proximity scores (from the diagnostic's classifier).

Mathematics Let:

  • $p_d$: discrete pathology from diag$pathology.
  • $p_p$: proximity pathology (argmax of .gdpar_orch_pathology_scores).
  • $k_d$, $k_p$: corresponding ladder keys (or NA if out-of-scope).

Decision logic:

  1. If entry_level is provided, use it directly (after validation).
  2. Else if level_map maps $p_d$ to a key, use that.
  3. Else if $k_d$ is NA (out-of-scope), return NA (caller will short-circuit to a certificate).
  4. Else if ($p_d \neq p_p$ or confidence $&lt; 0.6$) and $k_p$ is not NA, start at $\min(\text{index}(k_d), \text{index}(k_p))$.
  5. Else use $k_d$.

Returns A named list with:

  • key: Character string (ladder key or NA).
  • discrete: The discrete pathology.
  • proximity: The proximity pathology.
  • conflict: Logical, whether discrete and proximity disagree.
  • confident: Logical, whether confidence $\geq 0.6$.
  • source: Character string indicating the decision source ("user_entry_level", "user_level_map", or "combined_map").

Notes The key may be NA for out-of-scope pathologies. The caller uses this to decide whether to emit a certified limit. The algorithm is designed to be robust at the funnel/heavy-tail border.

.gdpar_orch_dense_mass(geom_target, fisher, reference, warmup_draws)

Purpose (internal). Constructs a dense, symmetric positive-definite mass matrix $M \in \mathbb{R}^{d \times d}$ for the euclidean_dense geometry level. It implements a three-tier fallback chain: (1) evaluate the caller-supplied Fisher information function at the reference, (2) invert the sample covariance of warmup draws, (3) compute the SoftAbs Riemannian mass at the reference. The result is SPD-floored before return.

Arguments

Argument Type Meaning
geom_target gdpar_geom_target The geometry target; its $dim field supplies $d$, the parameter dimension.
fisher function or NULL Optional fisher(theta) returning the expected Fisher information matrix at $\theta$.
reference numeric vector of length $d$ Position at which position-dependent metrics are evaluated.
warmup_draws numeric matrix or NULL Matrix of warmup draws (rows = draws, columns $\geq d$); used to compute a sample covariance when fisher fails.

Mathematics

The fallback chain is:

  1. If fisher is a function, compute $M \leftarrow \texttt{fisher}(\texttt{reference})$, coerced via as.matrix.
  2. If step 1 returned NULL and warmup_draws has $n &gt; d$ rows, compute the sample covariance $\hat{\Sigma} = \texttt{cov}(\texttt{warmup\_draws})$ and then $M \leftarrow \hat{\Sigma}^{-1}$ via chol2inv of the Cholesky factor obtained from .gdpar_geom_chol_spd.
  3. If both steps 1 and 2 returned NULL, fall back to the SoftAbs Riemannian metric of geom_target and evaluate its $mass at the reference.

Finally, the result is floored to the nearest SPD matrix with eigenvalue lower bound $10^{-8}$:

$$M \leftarrow \texttt{.gdpar_geom_floor_spd}(M,; 10^{-8})$$

Returns A dense numeric matrix of class matrix, dimension $d \times d$, guaranteed symmetric positive-definite with eigenvalues $\geq 10^{-8}$.

Notes

  • Every intermediate computation is wrapped in tryCatch; errors silently yield NULL and the chain advances to the next tier.
  • The SoftAbs fallback (tier 3) is always available as long as geom_target supports Riemannian geometry, making the function total.
  • warmup_draws must have strictly more rows than $d$ for tier 2 to activate; a matrix with $n \leq d$ rows is silently ignored.
  • fisher must return a square matrix of dimension $d$; if it returns something non-coercible, the error is caught and the fallback proceeds.

.gdpar_orch_build_metric(key, geom_target, fisher, reference, speed, rest_mass, warmup_draws)

Purpose (internal). Constructs the metric object for a given ladder key (one of the five geometry levels). This is the dispatcher that translates a string key into the appropriate gdpar_geom_metric_* object. It is armoured: errors are caught and reported in-band rather than propagated as crashes, so the controller can record the failure and escalate.

Arguments

Argument Type Meaning
key character string Ladder key: one of "euclidean_diagonal", "euclidean_dense", "riemannian", "relativistic", "sub_riemannian".
geom_target gdpar_geom_target The geometry target to sample.
fisher function or NULL Optional expected Fisher information function.
reference numeric vector of length $d$ Reference position for position-dependent metrics.
speed numeric scalar Relativistic speed parameter (used only when key = "relativistic").
rest_mass numeric scalar Relativistic rest-mass parameter (used only when key = "relativistic").
warmup_draws numeric matrix or NULL Warmup draws for the dense-mass fallback (used only when key = "euclidean_dense").

Mathematics

No explicit formula; the function dispatches on key via switch:

key Metric constructed
euclidean_diagonal gdpar_geom_metric_euclidean(dim = d)
euclidean_dense gdpar_geom_metric_euclidean(M = .gdpar_orch_dense_mass(...))
riemannian gdpar_geom_metric_riemannian(geom_target, curvature = "fisher"/"softabs", fisher = fisher)
relativistic gdpar_geom_metric_relativistic(geom_target, curvature = "fisher"/"softabs", fisher = fisher, speed = speed, rest_mass = rest_mass)
sub_riemannian gdpar_geom_metric_subriemannian(geom_target, fisher = fisher, reference = reference)

The curvature argument for the riemannian and relativistic levels is "fisher" when fisher is a function, "softabs" otherwise.

Returns A list with three elements:

  • ok (logical): TRUE if the metric was successfully built, FALSE on error.
  • metric (metric object or NULL): the constructed metric, or NULL on failure.
  • reason (character): NA_character_ on success, the conditionMessage of the caught error on failure.

Notes

  • If key is not one of the five recognized strings, a stop(sprintf("unknown ladder key '%s'", key)) is raised inside the tryCatch and captured, yielding ok = FALSE.
  • For sub_riemannian, fisher must be a function; otherwise a descriptive error is raised and caught. This is because the sub-Riemannian level structurally requires the expected Fisher information to define the wall/floor split.
  • The tryCatch wrapper means the caller (the orchestrator) never receives a thrown error from this function; it always receives the list.

.gdpar_orch_success_gate(fit, total, criteria, margin = 0)

Purpose (internal). Evaluates the multi-signal success gate for a single HMC fit. The gate checks three independent criteria (acceptance rate, divergence rate, and E-BFMI) and returns whether the fit passes all of them. The margin parameter implements hysteresis: for the decisive full gate a positive margin tightens the thresholds (making it harder to pass), while the probe gate passes margin = 0.

Arguments

Argument Type Meaning
fit list or NULL The HMC fit object; must contain $accept_rate, $n_divergent, $ebfmi. A NULL fit immediately yields a failure.
total integer Total number of proposals (warmup + sampling) against which the divergence rate is computed.
criteria list Named list with fields $accept_low, $accept_high, $divergent_rate_high, $ebfmi_low (see gdpar_geom_orchestrate_criteria).
margin numeric scalar Hysteresis margin; default 0 (probe gate). A positive value tightens thresholds.

Mathematics

Let $a = \texttt{fit\$accept\_rate}$, $n_d = \texttt{fit\$n\_divergent}$, $N = \texttt{total}$, $e = \texttt{fit$ebfmi}$, and $m = 1 + \texttt{margin}$.

The three pass/fail conditions are:

$$ \texttt{pass_accept}: \quad a \geq \texttt{accept_low} \cdot m ;;\wedge;; a \leq \texttt{accept_high} $$

$$ \texttt{pass_div}: \quad \frac{n_d}{N} \leq \frac{\texttt{divergent_rate_high}}{m} $$

$$ \texttt{pass_ebfmi}: \quad e ;\text{is finite} ;\wedge; e \geq \texttt{ebfmi_low} \cdot m $$

The gate passes if and only if all three conditions hold:

$$ \texttt{pass} = \texttt{pass_accept} ;\wedge; \texttt{pass_div} ;\wedge; \texttt{pass_ebfmi} $$

Note the asymmetric application of the hysteresis factor $m$: it raises the lower bound on acceptance and E-BFMI, and lowers the upper bound on divergence rate, but does not widen the upper acceptance bound.

Returns A list with fields:

  • pass (logical): whether the fit passed the gate.
  • accept (numeric): the acceptance rate, or NA_real_ if fit is NULL.
  • divergent_rate (numeric): $n_d / \max(N, 1)$, or NA_real_ if fit is NULL.
  • ebfmi (numeric): the E-BFMI, or NA_real_ if fit is NULL.
  • reasons (character vector): human-readable reasons for failure, or "all criteria met" on success, or "no fit" if fit is NULL.

Notes

  • The divergence rate denominator is $\max(\texttt{total}, 1)$ to avoid division by zero.
  • Each criterion uses isTRUE(...) which returns FALSE for NA values, so NA acceptance or NA E-BFMI will fail the gate.
  • The hysteresis margin tightens both the lower acceptance bound and the upper divergence bound, and raises the E-BFMI floor, but leaves the upper acceptance bound (accept_high) unchanged. This means a fit can pass the tight probe gate and still fail the looser full gate only if it violates the upper acceptance bound.
  • When fit is NULL, all returned values are NA and pass is FALSE.

.gdpar_orch_score(fit, total)

Purpose (internal). Computes a scalar "health" score for a single HMC fit. Used by the orchestrator for best-so-far tracking (the best result) and by the no-progress detector (to determine whether the current round improved on the previous best). Higher scores are better.

Arguments

Argument Type Meaning
fit list or NULL The HMC fit object; must contain $accept_rate, $n_divergent, $ebfmi.
total integer Total number of proposals (warmup + sampling).

Mathematics

The score is a linear combination of three statistics:

$$ \text{score} = e + 0.2 , a - 5 , r_d $$

where $e = \texttt{fit$ebfmi}$ (clamped to 0 if not finite), $a = \texttt{fit\$accept\_rate}$, and $r_d = \texttt{fit\$n\_divergent} / \max(\texttt{total}, 1)$ is the divergence rate.

The E-BFMI term has unit weight, acceptance is weighted at $0.2$, and divergences are heavily penalised at weight $5$.

Returns A numeric scalar: the health score, or $-\infty$ if fit is NULL or fit$accept_rate is NULL.

Notes

  • Non-finite E-BFMI values (e.g., NaN, Inf) are silently replaced with 0 before scoring, so a pathological energy geometry degrades to a penalty rather than producing a non-finite score.
  • The $-\infty$ sentinel for a missing fit ensures that any valid fit always outranks a NULL fit in comparisons.
  • The penalty weight of 5 for divergences means that a single divergence in 100 proposals ($r_d = 0.01$) reduces the score by $0.05$, roughly offsetting an acceptance rate improvement of $0.25$.

.gdpar_orch_tune_epsilon(geom_target, metric, budget, base_seed, round)

Purpose (internal). Performs a short, coarse step-size ($\epsilon$) search for a given geometry level. The orchestrator selects the geometry, but a level should not be failed merely because the default step size is ill-matched. This function probes a handful of candidate $\epsilon$ values on a geometric grid and returns the one whose acceptance is closest to a healthy target, penalising divergences. It is a deterministic, bounded analogue of the dual-averaging step-size adaptation used inside Stan.

Arguments

Argument Type Meaning
geom_target gdpar_geom_target The geometry target to sample.
metric metric object The metric for the current level; may contain $suggested_epsilon.
budget list Budget object; uses $epsilon (default step), $tune_iter (iterations per probe), $L (number of leapfrog steps).
base_seed integer Base seed for the orchestrator.
round integer Current round index (used for deterministic seeding).

Mathematics

The candidate $\epsilon$ grid is constructed from the base step $\epsilon_0$:

$$ \epsilon_0 = \begin{cases} \texttt{metric$suggested_epsilon} & \text{if available} \ \texttt{budget$epsilon} & \text{otherwise} \end{cases} $$

The raw grid is ${2\epsilon_0,; \epsilon_0,; 0.5\epsilon_0,; 0.25\epsilon_0,; 0.1\epsilon_0}$, then values are capped at 5 via pmin(..., 5), duplicates removed, and sorted ascending:

$$ G = \text{sort}!\bigl(\text{unique}!\bigl(\min(\epsilon_0 \cdot {2, 1, 0.5, 0.25, 0.1},; 5)\bigr)\bigr) $$

For each $\epsilon \in G$, a short HMC fit of $2 \cdot \texttt{tune\_iter}$ total proposals (warmup = sampling = tune_iter) is run with a deterministic seed. The gap metric is:

$$ \text{gap}(\epsilon) = \bigl|a(\epsilon) - 0.8\bigr| + 5 \cdot r_d(\epsilon) $$

where $a(\epsilon)$ is the acceptance rate and $r_d(\epsilon)$ is the divergence rate over $2 \cdot \texttt{tune\_iter}$ proposals. The $\epsilon$ with the smallest finite gap is selected.

Returns A list:

  • epsilon (numeric): the chosen step size (defaults to base if all probes fail or yield non-finite gaps).
  • fits (integer): number of HMC fits actually run (each probe, including failed ones, counts as one).
  • seconds (numeric): cumulative wall-clock time in seconds spent on all probes.

Notes

  • If a probe's gdpar_geom_hmc call throws an error, it is caught and silently skipped (the grid iteration advances).
  • The target acceptance rate is hardcoded at $0.8$.
  • The penalty weight for divergences (5) matches the score function.
  • The seed for each probe is derived deterministically via .gdpar_orch_seed(base_seed, round, 100L + k), where k is the 1-based grid index. The offset of 100 ensures these seeds do not collide with other slots in the orchestrator.
  • The function is bounded: at most $|G| \leq 5$ probes are run, each of tune_iter iterations, so the cost is predictable and small relative to a full sampling phase.

.gdpar_orch_seed(base, round, slot)

Purpose (internal). Produces a deterministic, reproducible integer seed from a base seed, a round index, and a slot index. The entire adaptive trajectory of the orchestrator is a deterministic function of the base seed, because every child seed is derived through this function. This is the resumability guarantee: re-running with the same base seed retraces the identical sequence of random numbers.

Arguments

Argument Type Meaning
base integer The caller-provided base seed.
round integer The current round index (0-based or 1-based, depending on the caller).
slot integer A slot identifier distinguishing different uses within a round (e.g., pilot, tune, full sampling, or epsilon probes).

Mathematics

$$ \text{seed} = \bigl(\texttt{base} + \texttt{round} \times 7919 + \texttt{slot} \times 104729\bigr) \bmod 2147483646 + 1 $$

The constants $7919$ and $104729$ are primes (the 1000th and 10000th primes, respectively), chosen to decorrelate the linear combinations. The modulus $2147483646 = 2^{31} - 2$ ensures the result fits in a 32-bit signed integer range suitable for R's RNG. The + 1 shifts the result to $[1, 2147483647]$, avoiding the seed 0 (which some RNGs treat specially).

Returns An integer scalar in $[1, 2147483647]$.

Notes

  • The result is passed through as.integer(...), so it is of R type integer.
  • Because the formula is linear in round and slot, distinct (round, slot) pairs can collide if the base seed is adversarial, but the prime multipliers make this practically negligible for reasonable seed spaces.
  • The modulus is not the standard Mersenne prime $2^{31} - 1$ but rather $2^{31} - 2$, yielding a range of $[1, 2^{31} - 1]$ after the + 1 offset.

.gdpar_orch_prescription(final_diag, ledger, fisher_available, budget, best_score, tried_keys)

Purpose (internal). Derives a conjectured, falsifiable prescription from the accumulated evidence of an orchestrator run that ended in a certified limit. Each prescription item is tagged as a "conjectured" status (following the ORPHEUS framework, section 16.3) and carries a falsifiable re-run test. The function implements the "smallest change predicted to break the limit" heuristic.

Arguments

Argument Type Meaning
final_diag list The final diagnostic object; fields used include $pathology, $difficulty_curve$grows_with_n, $difficulty_curve$slope, $cost$tractability.
ledger list of lists The per-round decision and sampling record; each element may contain $level and $metric_ok.
fisher_available logical Whether a Fisher information function was supplied by the caller.
budget list The budget object used for the run.
best_score numeric The best health score achieved across all rounds (from .gdpar_orch_score).
tried_keys character vector The ladder keys that were actually tried during the run.

Mathematics

No explicit formula; the logic is a rule-based cascade of conditional prescriptions (see Returns). The key decision variables are:

  • patho $= \texttt{final\_diag\$pathology}$: one of "multimodal", "boundary", "flat_direction", "quasi_deterministic", or other curvature pathologies.
  • grows $= \texttt{final\_diag\$difficulty\_curve\$grows\_with\_n}$: whether the diagnostic difficulty grows with sample size.
  • sub_sampled: whether the sub-Riemannian level was actually built and sampled (not just attempted).
  • tract $= \texttt{final\_diag\$cost\$tractability}$: the cost verdict from the diagnostic.

Returns A list of prescription items. Each item is a list with fields:

  • action (character): the recommended action.
  • rationale (character): why this action is recommended.
  • falsifiable_test (character): a re-run test that would falsify the prescription.
  • status (character): always "conjectured".

The prescription cascade (first match wins for out-of-scope pathologies):

  1. "multimodal" pathology $\Rightarrow$ single item recommending tempering / a multimodal sampler; short-circuit return.
  2. "boundary" pathology $\Rightarrow$ single item recommending reparametrisation of the boundary-pinned parameter; short-circuit return.
  3. "flat_direction" pathology $\Rightarrow$ single item recommending reparametrisation or elimination of the flat direction (Option A); short-circuit return.
  4. No Fisher + sub-Riemannian not sampled + geometry contracts ($\texttt{grows} \lor \texttt{patho} = \texttt{"quasi\_deterministic"}$) $\Rightarrow$ item recommending supply of the expected Fisher and enabling the sub-Riemannian level.
  5. Sub-Riemannian was sampled $\Rightarrow$ item recommending setting the sub-Riemannian reference to the warmup mode rather than the origin.
  6. Best score is finite and positive + tractability is "expensive" $\Rightarrow$ item recommending increased per-fit budget and epsilon/L tuning.
  7. Fallback (no items above) $\Rightarrow$ item recommending increased budget (more rounds / iterations) and revisiting epsilon / L.

Notes

  • The out-of-scope pathologies (multimodal, boundary, flat_direction) short-circuit immediately; they name the proper remedy outside the geometry ladder, deliberately avoiding spurious geometry advice.
  • The sub_sampled flag is computed by scanning the ledger: it is TRUE only if some round has level = "sub_riemannian" and metric_ok = TRUE (the metric was actually built successfully).
  • The function uses a local add closure to append items, modifying the outer items list via <<-.
  • The budget and tried_keys arguments are accepted but not explicitly referenced in the current implementation (they may be used in future expansions or are present for API consistency).

.gdpar_orch_certificate(verdict, final_diag, ledger, culprit, prescription, reproducibility, budget_spent)

Purpose (internal). Assembles the certified-limit certificate object. This is the structured output when the orchestrator exhausts its budget without achieving a passing gate, or when the diagnostic identifies an out-of-scope pathology. The certificate has three rigour layers (algebraic, statistical, numerical) plus a conjectured prescription and reproducibility metadata.

Arguments

Argument Type Meaning
verdict character The overall verdict string (e.g., "certified_limit", "out_of_scope").
final_diag list The final diagnostic object; fields used include $pathology, $signals$condition_number, $difficulty_curve$slope, $difficulty_curve$grows_with_n.
ledger list of lists The per-round decision and sampling record (the statistical layer).
culprit character or list Identifies the primary cause or level at which the limit was reached.
prescription list of lists The prescription items from .gdpar_orch_prescription.
reproducibility list Reproducibility metadata (base seed, round structure, etc.).
budget_spent list Numerical accounting: total fits, total seconds, rounds completed, etc.

Mathematics

No formula; this is a pure data-assembly function.

Returns A list of class c("gdpar_geom_certificate", "list") with the following structure:

  • $verdict: the verdict string.
  • $algebraic: the algebraic rigour layer —
    • $pathology: the identified pathology class.
    • $condition_number: the final condition number from the signals.
    • $difficulty_slope: the difficulty-vs-$n$ slope.
    • $grows_with_n: logical, whether difficulty grows with sample size.
    • $culprit: the culprit identifier.
  • $statistical: the full ledger (per-level sampling diagnostics).
  • $numerical: the budget_spent accounting.
  • $prescription: the list of conjectured prescription items.
  • $reproducibility: the reproducibility metadata.

Notes

  • The class is set to c("gdpar_geom_certificate", "list"), enabling S3 dispatch (e.g., a print.gdpar_geom_certificate method if one exists).
  • The condition_number is extracted as the last row of final_diag$signals$condition_number, i.e., the most recent diagnostic signal.
  • The function is purely structural: it does no computation, validation, or transformation beyond field extraction and list assembly.

.gdpar_orch_attach_laplace(obj, geom_target, reference, draws, seed, inform)

Purpose (internal). Implements the opt-in Laplace fallback (RG.7 step 4, D99). When a run ends in a certified limit and the caller opted in via laplace_fallback = TRUE, this function attempts to attach a Laplace approximation to the certificate and relabels the status accordingly. The function is conservative: on the out-of-scope path it only attaches the Laplace when the curvature at the mode is genuinely positive-definite, to avoid presenting a Gaussian approximation when the Gaussian premise is violated.

Arguments

Argument Type Meaning
obj list The orchestrator result object; must have $status and will gain $laplace on success.
geom_target gdpar_geom_target The geometry target for the Laplace computation.
reference numeric vector of length $d$ The position (mode) at which the Laplace is computed.
draws integer Number of iid draws from the Laplace approximation to carry (0 = mode + precision only).
seed integer Seed for reproducibility of the Laplace draws.
inform function Callback for informational messages (e.g., message).

Mathematics

The function calls gdpar_geom_laplace(geom_target, reference, draws, seed) which computes the Laplace approximation around reference:

$$ q(\theta) \approx \mathcal{N}!\bigl(\mu = \texttt{reference},; \Sigma = [-\nabla^2 \log p(\texttt{reference})]^{-1}\bigr) $$

The returned object carries $all_pos (whether the Hessian is positive-definite at the reference) and $fit_quality_label (a human-readable fidelity diagnosis).

Returns The modified obj (the same list, mutated in place by adding $laplace and updating $status), or the original obj unchanged if the Laplace computation fails or the Gaussian premise is violated.

Status relabeling logic:

Original $status Condition New $status
"certified_limit" Laplace succeeds "certified_limit_laplace"
"out_of_scope" Laplace succeeds and $all_pos = TRUE "out_of_scope_laplace"
"out_of_scope" Laplace succeeds and $all_pos = FALSE unchanged (no attachment)
any Laplace throws an error unchanged (warning emitted)

Notes

  • If gdpar_geom_laplace throws an error, it is caught, a warning of class "gdpar_geometry_warning" is emitted via gdpar_warn, and the original obj is returned unmodified. This makes the function fail-safe.
  • On the "out_of_scope" path, if the curvature at the mode is not positive-definite ($all_pos = FALSE), the function deliberately does not attach the Laplace and prints an informational message explaining why. This prevents overreach: the Gaussian premise does not hold for multimodal, boundary, or flat-direction pathologies.
  • The inform callback is called with a formatted message indicating which fit_quality_label was attached and what the new status is.
  • When draws = 0 (the default for laplace_draws), the Laplace object contains only the mode and the precision matrix; no random draws are generated.
  • The function mutates obj by direct assignment (obj$laplace <- lap, obj$status <- ...) before returning it, so the caller receives the modified object.

gdpar_geom_orchestrate(target, geom_target = NULL, fisher = NULL, reference = NULL, level_map = NULL, entry_level = NULL, budget = NULL, criteria = NULL, speed = 10, rest_mass = 1, n_grid = NULL, checkpoint_dir = NULL, laplace_fallback = FALSE, laplace_draws = 0L, seed = 20260603L, verbose = TRUE, ...)

Purpose
Main exported orchestrator for geometry-adaptive parameter estimation. It implements a closed-loop control system that diagnoses the geometry of a sampling target, selects and evaluates parameterization levels on a ladder of increasing complexity, and applies a two-phase probe-and-full sampling strategy (successive-halving spirit) under a configurable computational budget. The function aims to find a parameterization level that passes a success gate, certifying intractability if none are found within the budget.

Arguments

  • target (object): The primary diagnostic target. Can be a gdpar_geometry_target object, a list with a model (e.g., a CmdStanModel) and dim, or another object coercible to a geometry target. Its structure determines how geom_target is derived if not supplied.
  • geom_target (list or NULL): The explicit engine sampling target, a gdpar_geom_target object. If NULL, derived from target via .gdpar_orch_derive_geom_target.
  • fisher (function or NULL): A function of parameter vector theta returning the Fisher information matrix (or its approximation). Used to construct a Riemannian metric at certain ladder levels. Must be NULL or a function.
  • reference (numeric vector or NULL): The reference point for parameterization anchoring. Its length must match the dimension d of geom_target. Defaults to a zero vector of length d.
  • level_map (list or NULL): A mapping from diagnosed pathology codes to specific ladder level keys. Used by .gdpar_orch_select_entry to override the default pathology-to-level mapping.
  • entry_level (character or NULL): A specific ladder level key to start the orchestration, bypassing the initial diagnosis-based selection.
  • budget (list or NULL): Computational budget controls. If NULL, defaults to gdpar_geom_orchestrate_budget(). Elements include max_rounds, max_fits, max_seconds, max_seconds_per_fit, max_levels, L, probe_warmup, probe_iter, full_warmup, full_iter, epsilon, tune_epsilon, hysteresis, stall_limit, n_rediagnose.
  • criteria (list or NULL): Success gate criteria. If NULL, defaults to gdpar_geom_orchestrate_criteria(). Used by .gdpar_orch_success_gate to evaluate sampler output.
  • speed (numeric scalar): Parameter controlling the speed of the dynamics (e.g., in the kinetic energy). Must be non-negative.
  • rest_mass (numeric scalar): Parameter controlling the rest mass (inverse scale) in the kinetic energy. Must be non-negative.
  • n_grid (integer or NULL): Number of grid points for the initial geometry diagnostic (passed to gdpar_geometry_diagnostic). NULL uses the diagnostic's default.
  • checkpoint_dir (character or NULL): Path to a directory. If non-NULL, the ledger is atomically saved here after each round via .gdpar_orch_checkpoint.
  • laplace_fallback (logical): If TRUE, attaches a Laplace approximation fallback to the output object when the orchestration terminates (either out_of_scope or certified_limit).
  • laplace_draws (integer): Number of Laplace draws to take if laplace_fallback is TRUE.
  • seed (integer): Master seed for reproducibility. Used to generate per-round, per-slot seeds via .gdpar_orch_seed(seed, round, slot).
  • verbose (logical): If TRUE, emits progress messages via gdpar_inform.
  • ... (additional arguments): Passed to gdpar_geometry_diagnostic.

Mathematics

The orchestration implements a monotonically increasing ladder index $i$ with a floor $f$ such that $f = \max{f, i}$ after each round. The ladder is defined by .gdpar_orch_ladder(), where each level has an index. The controller cycles through levels $k$ selected by:

  1. Initial diagnosis $D_0$ via gdpar_geometry_diagnostic.
  2. Selection $s = \text{select\_entry}(D_0)$.
  3. For each round $r = 1, \dots, R_{\max}$: a. Budget admission check: Stop if estimated remaining time exceeds max_seconds or if max_fits reached. b. Cycle guard: If state $(s, \text{"probe"})$ visited, escalate to next level above floor $f$. c. Build metric $M_s$ for level $s$. d. Step-size tuning (optional): $\epsilon \leftarrow \text{tune\_epsilon}(\text{target}, M_s)$. e. Probe: sample $n_p = \text{probe\_warmup} + \text{probe\_iter}$ draws with metric $M_s$, step size $\epsilon$, and leapfrog steps $L$. f. Gate probe: $G_p = \text{success\_gate}(\text{probe\_fit}, n_p, \text{criteria}, 0)$. g. If $G_p.\text{pass}$, full phase: sample $n_f = \text{full\_warmup} + \text{full\_iter}$ draws. Gate full: $G_f = \text{success\_gate}(\text{full\_fit}, n_f, \text{criteria}, \text{hysteresis})$. h. Update best score: $\text{score} = \text{score\_fit}(\cdot)$. i. If $G_f.\text{pass}$, return resolved. j. Else, update floor $f = \max(f, \text{index}(s))$, re-diagnose $D_r$ with fresh seed, choose next level $s'$:
    • If re-diagnosis unstable (checked with n_rediagnose >= 2), choose safe level above $f$.
    • Else, map pathology to candidate $c$. If $c$ is valid and $\text{index}(c) &gt; f$, choose $c$; else safe level.
    • If $s' == s$, increment stall counter; if stall $\ge$ stall_limit, stop. k. Loop until budget exhausted or top of ladder reached.
  4. If loop ends without resolution, output certified_limit with a prescription.

The metric construction at level $s$ uses fisher and reference to define a Riemannian metric (e.g., Euclidean, Fisher-Rao, or Hessian-based) depending on the level's definition.

Returns
A list of class c("gdpar_geom_orchestration", "list") with elements:

  • status (character): One of "out_of_scope", "resolved", "certified_limit".
  • level (character): Ladder level key at which resolution occurred (if "resolved").
  • level_index (integer): Index of the resolved level.
  • metric (list): The metric object used at the resolved level.
  • draws (matrix): Post-warmup draws from the final successful fit (if "resolved").
  • fit (list): The final successful HMC fit object.
  • gate (list): The success gate result for the final fit.
  • certificate (list): A certificate object (for "out_of_scope" or "certified_limit" statuses).
  • ledger (list): A list of per-round records with diagnostics, gates, acceptance rates, EBFMI, timing, etc.
  • diagnosis (list): The initial geometry diagnostic object (diag0).
  • best (list): Best fit found during orchestration, with elements score, draws, fit, key.
  • budget_spent (list): fits (integer count) and seconds (total elapsed time).
  • reproducibility (list): seed, gdpar_version, budget, criteria, initial_pathology, dim.
  • laplace_fit (list, optional): Attached if laplace_fallback = TRUE and status is "out_of_scope" or "certified_limit".

Notes

  • The function performs validation on seed, speed, rest_mass, and fisher (must be a function or NULL). Raises "gdpar_input_error" on failure.
  • Emits an initial message stating it is an "opt-in geometry-adaptive controller; it does not modify any fit."
  • Uses setTimeLimit(elapsed = cap, transient = TRUE) to enforce per-fit time limits; timeout is caught as an error and recorded.
  • Checkpointing writes the ledger atomically via .gdpar_orch_checkpoint (temp file then rename).
  • The floor_index acts as a monotone ratchet: the controller never considers levels with index $\le f$.
  • The stall counter resets when the selected level changes; triggers stop at budget$stall_limit consecutive stalls.
  • The warmup_draws are updated from the best probe or full fit to warm-start subsequent metric constructions (via reference update and passing draws to .gdpar_orch_build_metric).
  • The best object tracks the highest-scoring fit seen across all rounds, regardless of gate outcome.
  • The ... arguments are passed to gdpar_geometry_diagnostic, allowing custom diagnostic settings.
  • If geom_target is not a gdpar_geom_target, it is coerced via gdpar_geom_target(geom_target).
  • The reference is updated in-place to the last row of the best fit's draws (if finite) to potentially guide future metric constructions.

.gdpar_orch_derive_geom_target(target)

Purpose
Internal helper that constructs a gdpar_geom_target object from the primary target argument when geom_target is not supplied by the user. It handles two supported target types: a gdpar_geometry_target (from the diagnostic phase) or a list containing a Stan model (CmdStanModel) with dim and optionally data.

Arguments

  • target (object): The diagnostic target. Must be either:
    • A gdpar_geometry_target object (has a make method and default parameters).
    • A list with a model element (e.g., a CmdStanModel) and a dim element (integer, the unconstrained dimension). An optional data element is passed to gdpar_geom_target.

Mathematics
None.

Returns
A gdpar_geom_target object, created by calling gdpar_geom_target on the appropriate object(s).

Notes

  • Raises a "gdpar_input_error" if the target type is unrecognized or if a Stan-model target lacks dim.
  • For a gdpar_geometry_target, it calls target$make(target$default_n, target$default_difficulty) to instantiate the model, then passes the result to gdpar_geom_target.
  • For a Stan-model list target, it passes object = target$model, dim = target$dim, and data = target$data to gdpar_geom_target.
  • The error message for unrecognized targets suggests supplying geom_target explicitly.

.gdpar_orch_checkpoint(dir, ledger)

Purpose Internal helper that atomically persists the orchestrator's ledger (the running log of level attempts) to disk. It writes to a temporary file first, then renames, so that a crash mid-write cannot corrupt the canonical ledger file. Called by the orchestrator loop to create recoverable checkpoints.

Arguments

Argument Type Meaning
dir character(1) Directory in which to store the checkpoint. Created recursively if it does not exist.
ledger list (typically a list of per-level result records) The current state of the orchestration ledger to persist.

Returns invisible(NULL). Called purely for its side effect of writing orchestrate_ledger.rds into dir.

Notes

  • Atomic write: The function writes to a temporary .rds file (via tempfile with tmpdir = dir) and then atomically renames it to orchestrate_ledger.rds using file.rename. This ensures that readers of orchestrate_ledger.rds never observe a partially written file.
  • Directory creation: dir.create is called with recursive = TRUE and showWarnings = FALSE, so nested directories are created silently and an existing directory does not raise a warning.
  • No return value on purpose: invisible(NULL) signals that callers should not depend on a return value; the meaningful output is the file on disk.

print.gdpar_geom_orchestration(x, ...)

Purpose S3 print method for objects of class gdpar_geom_orchestration. Provides a concise, human-readable summary of the orchestrator's final status—whether it resolved successfully or issued a certificate of failure—including key diagnostics and budget consumption.

Arguments

Argument Type Meaning
x gdpar_geom_orchestration (list) The orchestration result object produced by gdpar_geom_orchestrate.
... (unused) Generic ellipsis for S3 method compatibility; not consumed.

Returns invisible(x). The object is printed as a side effect and returned invisibly so that it can be piped.

Notes

  • Branching on x$status:
    • If identical(x$status, "resolved"), the printout includes:
      • The resolution level and its index (x$level, x$level_index).
      • The dimensions of the draw matrix (nrow(x$draws) × ncol(x$draws)).
      • The acceptance rate (x$fit$accept_rate) and E-BFMI (x$fit$ebfmi), both formatted to 3 significant digits.
    • Otherwise (failure / certificate case), the printout includes:
      • The certificate verdict (x$certificate$verdict).
      • The chain of levels attempted, extracted by vapply over x$ledger pulling each record's $level field, deduplicated with unique, and joined with " -> ".
      • The number of prescription items in x$certificate$prescription.
  • Laplace fallback block: If x$laplace is not NULL, a secondary block is printed showing the Laplace approximation's fit-quality label, the (unfloored) condition number (x$laplace$cond_unfloored), and the Pareto-$k$ diagnostic (x$laplace$fit_quality$pareto_k), each formatted to 3 significant digits.
  • Budget summary: Always prints x$budget_spent$fits (number of model fits) and x$budget_spent$seconds (wall-clock time, 3 significant digits).
  • Side effect: Output is written to the console via cat.

print.gdpar_geom_certificate(x, ...)

Purpose S3 print method for objects of class gdpar_geom_certificate. Summarises the structured certificate emitted by gdpar_geom_orchestrate when the geometry-adaptive ladder cannot be resolved within the prescribed budget or when diagnosis points outside the ladder entirely. The certificate separates evidence into rigour layers and attaches a falsifiable prescription.

Arguments

Argument Type Meaning
x gdpar_geom_certificate (list) The certificate object. Expected to contain the named elements verdict, algebraic, statistical, numerical, prescription, and (optionally) reproducibility.
... (unused) Generic ellipsis for S3 method compatibility; not consumed.

Returns invisible(x). The certificate is printed as a side effect and returned invisibly.

Notes

  • Algebraic layer (x$algebraic): Printed fields are:
    • pathology — a label describing the detected geometric pathology (e.g., ill-conditioning).
    • condition_number — the condition number of the relevant matrix, formatted to 3 significant digits.
    • grows_with_n — logical flag indicating whether the pathology worsens with increasing sample size $n$.
  • Statistical layer (x$statistical): Only the count of recorded levels is printed (length(x$statistical)), not the per-level details themselves.
  • Numerical layer (x$numerical): Total number of fits and total seconds consumed, the latter formatted to 3 significant digits.
  • Prescription (x$prescription): Iterated with for; each element is expected to be a list with at least an $action field, which is printed as a bullet. The prescription represents conjectured, falsifiable fixes to the identified pathology.
  • Side effect: All output goes to the console via cat.

R/geometry_suite.R

gdpar_geometry_suite(which = NULL)

Purpose Exported entry point that assembles and returns the Block RG calibration suite — a named registry of eight synthetic target distributions (G0–G7), each engineered to exhibit one row of the posterior-geometry pathology taxonomy with a known ground truth. The suite is the falsifiable backbone against which gdpar_geometry_diagnostic is calibrated.

Arguments

Argument Type Meaning
which NULL or character vector Optional subset of target ids to select. When NULL (default), all eight targets are returned.

Returns A named list of gdpar_geometry_target objects. The names are exactly G0_isotropic, G1_anisotropic, G2_funnel, G3_heavy_tails, G4_quasi_deterministic, G5_multimodal, G6_boundary, G7_flat_direction. When which is non-NULL, the returned list is the sub-list registry[which], preserving the original names and order.

Notes

  • If which is supplied but is not a character vector, or if any element is not among names(registry), the function calls gdpar_abort with a message of the form "Argument 'which' must be a subset of: G0_isotropic, G1_anisotropic, …" and error class "gdpar_input_error".
  • The eight registry entries are built by the internal constructors .gdpar_geom_G0() through .gdpar_geom_G7() (only G0–G4 are defined in this section; G5–G7 are referenced but their definitions lie in the subsequent section).

.gdpar_geom_target(id, label, pathology, geometry_remedy, culprit, difficulty_scales_with_n, bounds, default_n, default_difficulty, n_grid, make)

Purpose Internal constructor helper that bundles the static metadata fields of a synthetic target together with its make instance-constructor and stamps the S3 class gdpar_geometry_target on the result. Every .gdpar_geom_G*() factory delegates to this helper to produce its registry entry.

Arguments

Argument Type Meaning
id character Short identifier of the target (e.g. "G0_isotropic").
label character Human-readable descriptive label.
pathology character Pathology category from the Block RG taxonomy.
geometry_remedy character The geometry level that remedies the pathology (the "remedy" column).
culprit character or NA_character_ Name(s) of the culprit parameter(s), or NA_character_ when no single parameter is implicated.
difficulty_scales_with_n logical Ground-truth flag: does the difficulty grow with the size knob n?
bounds named list or NULL Constrained-scale parameter bounds; NULL when the target is unconstrained.
default_n numeric/integer Default value of the size knob n.
default_difficulty numeric/integer Default difficulty value.
n_grid numeric/integer vector Suggested sweep of n values for measuring the difficulty-vs-n curve.
make function Instance constructor with signature make(n, difficulty) returning a target instance list.

Returns A list with elements id, label, pathology, geometry_remedy, culprit, difficulty_scales_with_n, bounds, default_n, default_difficulty, n_grid, make, assigned S3 class c("gdpar_geometry_target", "list").

Notes

  • The leading comment states the helper "validates the make() contract lightly so a malformed target surfaces early," but the code as written performs no validation: it only assembles the list and assigns the class. A malformed make would not be caught here.

.gdpar_geom_G0()

Purpose Builds the G0_isotropic registry entry: an isotropic standard normal serving as the easy negative control. The prescribed remedy is the Euclidean diagonal metric (the current default).

Arguments None.

Mathematics The make(n, difficulty) closure sets $d = \text{as.integer}(\text{difficulty})$ and defines a $d$-dimensional standard normal on the unconstrained scale. The log-density (up to an additive constant) and its gradient are:

$$ \log p(\theta) = -\frac{1}{2}\sum_{i=1}^{d}\theta_i^{2}, \qquad \nabla\log p(\theta) = -\theta . $$

The Stan program is:

data { int<lower=1> d; }
parameters { vector[d] theta; }
model { theta ~ std_normal(); }

Returns A gdpar_geometry_target with static fields:

Field Value
id "G0_isotropic"
label "Isotropic standard normal (easy control)"
pathology "isotropic"
geometry_remedy "euclidean_diagonal"
culprit NA_character_
difficulty_scales_with_n FALSE
bounds NULL
default_n 1
default_difficulty 5
n_grid c(1, 4, 16)

The make(n = 1, difficulty = 5) closure returns an instance list with components stan_code (character), stan_data = list(d = d), log_prob (function), grad_log_prob (function), dim = d, and param_names = paste0("theta[", seq_len(d), "]").

Notes

  • The argument n is accepted but ignored; only difficulty is used (it sets the dimension $d$).
  • stan_data contains only d.

.gdpar_geom_G1()

Purpose Builds the G1_anisotropic registry entry: a diagonal normal with a fixed, large condition number $\kappa$ that does not grow with n — a straight canyon. The prescribed remedy is a constant dense Euclidean metric.

Arguments None.

Mathematics The make(n, difficulty) closure fixes $d = 5$ and sets $\kappa = \text{as.numeric}(\text{difficulty})$. The coordinate variances are geometrically spaced so that the ratio of the largest to the smallest variance is exactly $\kappa$:

$$ \sigma_i^{2} = \exp!\Bigl(\tfrac{i-1}{d-1}\log\kappa\Bigr), \qquad i = 1,\dots,d . $$

Equivalently, vars <- exp(seq(0, log(kappa), length.out = d)). With $\sigma_i = \sqrt{\sigma_i^{2}}$ and inverse variances $\sigma_i^{-2}$, the log-density and gradient are:

$$ \log p(\theta) = -\frac{1}{2}\sum_{i=1}^{d}\frac{\theta_i^{2}}{\sigma_i^{2}}, \qquad \nabla_i\log p(\theta) = -\frac{\theta_i}{\sigma_i^{2}} . $$

The Stan program is:

data { int<lower=1> d; vector<lower=0>[d] sigma; }
parameters { vector[d] theta; }
model { theta ~ normal(0, sigma); }

Returns A gdpar_geometry_target with static fields:

Field Value
id "G1_anisotropic"
label "Anisotropic normal, fixed condition number"
pathology "anisotropic"
geometry_remedy "euclidean_dense"
culprit NA_character_
difficulty_scales_with_n FALSE
bounds NULL
default_n 1
default_difficulty 100
n_grid c(1, 4, 16)

The make closure returns an instance with stan_data = list(d = d, sigma = sigma), dim = 5L, and param_names = paste0("theta[", seq_len(d), "]").

Notes

  • n is ignored; the dimension is fixed at $d = 5$.
  • difficulty is interpreted as the condition number $\kappa$ (variance ratio), defaulting to 100.

.gdpar_geom_G2()

Purpose Builds the G2_funnel registry entry: Neal's funnel, in which a log-scale variable v governs the spread of the remaining coordinates, producing position-dependent curvature. The prescribed remedy is a position-dependent (Riemannian) metric.

Arguments None.

Mathematics The make(n, difficulty) closure fixes $d = 10$ and sets $\text{scale\_v} = \text{as.numeric}(\text{difficulty})$. The parameter vector is $\theta = (v,, x_1,\dots,x_{d-1})$ with the generative model

$$ v \sim \mathcal{N}(0,,\text{scale_v}^{2}), \qquad x_i \mid v \sim \mathcal{N}(0,, e^{v}) . $$

On the unconstrained scale, the log-density (up to an additive constant) is

$$ \log p(\theta) = -\frac{1}{2}\Bigl(\frac{v}{\text{scale_v}}\Bigr)^{2} - \frac{1}{2}(d-1),v - \frac{1}{2},e^{-v}\sum_{i=1}^{d-1}x_i^{2} . $$

The gradient components are

$$ \frac{\partial\log p}{\partial v} = -\frac{v}{\text{scale_v}^{2}} - \frac{1}{2}(d-1) + \frac{1}{2},e^{-v}\sum_{i=1}^{d-1}x_i^{2}, \qquad \frac{\partial\log p}{\partial x_i} = -x_i,e^{-v} . $$

The Stan program is:

data { int<lower=1> d; real<lower=0> scale_v; }
parameters { real v; vector[d - 1] x; }
model {
  v ~ normal(0, scale_v);
  x ~ normal(0, exp(v / 2));
}

Returns A gdpar_geometry_target with static fields:

Field Value
id "G2_funnel"
label "Neal's funnel (variable curvature)"
pathology "funnel"
geometry_remedy "riemannian"
culprit "v"
difficulty_scales_with_n FALSE
bounds NULL
default_n 1
default_difficulty 3
n_grid c(1, 4, 16)

The make closure returns an instance with stan_data = list(d = d, scale_v = scale_v), dim = 10L, and param_names = c("v", paste0("x[", seq_len(d - 1), "]")).

Notes

  • n is ignored; the dimension is fixed at $d = 10$ (so there are $d - 1 = 9$ x coordinates).
  • difficulty is interpreted as the prior standard deviation scale_v of the funnel variable v.
  • The gradient closure returns c(gv, gx) where gv is the scalar derivative w.r.t. v and gx is the length-$(d-1)$ vector of derivatives w.r.t. x.

.gdpar_geom_G3()

Purpose Builds the G3_heavy_tails registry entry: independent Student-$t$ distributions with a small number of degrees of freedom, exhibiting directional heaviness. The prescribed remedy is a non-inner-product (Finsler / relativistic) kinetic energy.

Arguments None.

Mathematics The make(n, difficulty) closure fixes $d = 5$ and sets $\nu = \text{as.numeric}(\text{difficulty})$. Each coordinate is independently $t_\nu(0, 1)$. Dropping the normalizing constant, the log-density is

$$ \log p(\theta) = -\frac{\nu+1}{2}\sum_{i=1}^{d}\log!\Bigl(1 + \frac{\theta_i^{2}}{\nu}\Bigr) , $$

and the gradient is

$$ \nabla_i\log p(\theta) = -\frac{(\nu+1),\theta_i}{\nu + \theta_i^{2}} . $$

The Stan program is:

data { int<lower=1> d; real<lower=0> nu; }
parameters { vector[d] theta; }
model { theta ~ student_t(nu, 0, 1); }

Returns A gdpar_geometry_target with static fields:

Field Value
id "G3_heavy_tails"
label "Heavy tails (independent Student-t)"
pathology "heavy_tails"
geometry_remedy "finsler_relativistic"
culprit NA_character_
difficulty_scales_with_n FALSE
bounds NULL
default_n 1
default_difficulty 2
n_grid c(1, 4, 16)

The make closure returns an instance with stan_data = list(d = d, nu = nu), dim = 5L, and param_names = paste0("theta[", seq_len(d), "]").

Notes

  • n is ignored; the dimension is fixed at $d = 5$.
  • difficulty is interpreted as the degrees-of-freedom parameter $\nu$; low $\nu$ (default 2) yields heavier tails.
  • The gradient uses log1p inside the log-density but the closed-form derivative avoids log1p entirely.

.gdpar_geom_G4()

Purpose Builds the G4_quasi_deterministic registry entry: a near-degenerate canyon in which $d - 1$ directions are pinned with variance $1/n$ while one direction stays free, so the condition number grows with n and the posterior collapses toward a lower-dimensional manifold. This is the eBird count case; the prescribed remedy is a sub-Riemannian (distribution-of-directions) treatment.

Arguments None.

Mathematics The make(n, difficulty) closure sets $d = \text{as.integer}(\text{difficulty})$ and $n_{\text{eff}} = \text{as.numeric}(n)$. The inverse variances are

$$ \frac{1}{\sigma_i^{2}} = \begin{cases} 1 & i = 1, \ n_{\text{eff}} & i = 2,\dots,d. \end{cases} $$

Thus the first direction has variance $1$ and the remaining $d - 1$ directions have variance $1/n_{\text{eff}}$, giving a condition number of $n_{\text{eff}}$. The log-density and gradient are

$$ \log p(\theta) = -\frac{1}{2}\sum_{i=1}^{d}\frac{\theta_i^{2}}{\sigma_i^{2}} = -\frac{1}{2}\theta_1^{2} - \frac{n_{\text{eff}}}{2}\sum_{i=2}^{d}\theta_i^{2}, \qquad \nabla_i\log p(\theta) = -\frac{\theta_i}{\sigma_i^{2}} . $$

The Stan program is:

data { int<lower=1> d; real<lower=0> n_eff; }
parameters { vector[d] theta; }
model {
  theta[1] ~ normal(0, 1);
  theta[2:d] ~ normal(0, inv_sqrt(n_eff));
}

Returns A gdpar_geometry_target with static fields:

Field Value
id "G4_quasi_deterministic"
label "Quasi-deterministic canyon (condition number grows with n)"
pathology "quasi_deterministic"
geometry_remedy "sub_riemannian"
culprit NA_character_
difficulty_scales_with_n TRUE
bounds NULL
default_n 100
default_difficulty 5
n_grid c(10, 100, 1000)

The make closure returns an instance with stan_data = list(d = d, n_eff = n_eff), dim = d, and param_names = paste0("theta[", seq_len(d), "]").

Notes

  • This is the only target in this section for which difficulty_scales_with_n is TRUE: n is the genuine size knob driving the condition number.
  • difficulty controls the dimension $d$ (default 5), while n controls the inverse variance of the pinned directions (default 100).
  • The n_grid c(10, 100, 1000) spans three orders of magnitude, reflecting the structural (not merely quasi-deterministic) nature of the pathology.
  • The inv_var vector is constructed as c(1, rep(n_eff, d - 1)); when $d = 1$ this reduces to c(1) (no pinned directions), though the default difficulty = 5 gives $d = 5$.

.gdpar_geom_G5()

Purpose

Factory function that constructs and registers the "G5_multimodal" geometry target — a $d=3$ dimensional target whose density is an equal mixture of two multivariate normal modes centred at $\pm\text{sep}$ along every coordinate. This target exercises multimodality; the prescribed remedy is tempering.

Arguments

This function takes no arguments.

Mathematics

The target density (up to additive constants) is the equal-weight mixture

$$ p(\theta) = 0.5,\mathcal{N}(\theta \mid -\mathrm{sep}\cdot\mathbf{1}_d,, I_d) + 0.5,\mathcal{N}(\theta \mid +\mathrm{sep}\cdot\mathbf{1}_d,, I_d), $$

with $d = 3$. The R-side log_prob omits the normalising constant $-\frac{d}{2}\log(2\pi)$ from each component and computes

$$ \ell_1 = -\frac{1}{2}\sum_{i=1}^{d}(\theta_i + \mathrm{sep})^2, \qquad \ell_2 = -\frac{1}{2}\sum_{i=1}^{d}(\theta_i - \mathrm{sep})^2, $$

then applies the log-sum-exp stabilisation with $m = \max(\ell_1, \ell_2)$:

$$ \log p(\theta) = \log(0.5) + m + \log!\bigl(e^{\ell_1 - m} + e^{\ell_2 - m}\bigr). $$

The gradient uses the posterior responsibilities (mixture weights)

$$ w_1 = \frac{e^{\ell_1 - m}}{e^{\ell_1 - m} + e^{\ell_2 - m}}, \qquad w_2 = 1 - w_1, $$

giving

$$ \nabla_\theta \log p(\theta) = w_1\bigl[-(\theta + \mathrm{sep})\bigr] + w_2\bigl[-(\theta - \mathrm{sep})\bigr]. $$

The Stan code declares d and sep as data, theta as a vector[d] parameter, and adds log_mix(0.5, normal_lpdf(theta | -sep, 1), normal_lpdf(theta | sep, 1)) to the target.

Returns

The return value of .gdpar_geom_target(...), i.e. a gdpar_geometry_target object. The embedded make closure (with signature make(n = 1, difficulty = 4)) returns a list containing:

Element Type / Structure
stan_code Character scalar (Stan model string)
stan_data list(d = 3L, sep = <difficulty>)
log_prob Function: function(theta) → numeric scalar
grad_log_prob Function: function(theta) → numeric vector of length d
dim Integer scalar 3L
param_names Character vector c("theta[1]", "theta[2]", "theta[3]")

The registered metadata is: id = "G5_multimodal", label = "Two well-separated normal modes", pathology = "multimodal", geometry_remedy = "tempering", culprit = NA_character_, difficulty_scales_with_n = FALSE, bounds = NULL, default_n = 1, default_difficulty = 4, n_grid = c(1, 4, 16).

Notes

  • difficulty is coerced to numeric via as.numeric(difficulty) and used directly as sep; larger values place the modes farther apart, increasing multimodality severity.
  • The n argument of make is accepted but unused — the geometry is structural and does not scale with n (difficulty_scales_with_n = FALSE).
  • The R-side log_prob and grad_log_prob omit normalising constants that the Stan normal_lpdf includes; they are self-consistent (gradient matches the log-density) but not identical to the Stan log-density up to a constant.
  • culprit is NA_character_ (no single parameter is responsible).

.gdpar_geom_G6()

Purpose

Factory function that constructs and registers the "G6_boundary" geometry target — a $d=4$ dimensional target where one coordinate $b$ is constrained to $(0,1)$ with a $\mathrm{Beta}(\alpha, 1)$ prior ($\alpha &lt; 1$), pushing mass against the lower boundary, while the remaining $d-1 = 3$ coordinates are standard normal nuisances. The R closure operates on the unconstrained scale $\phi_b = \mathrm{logit}(b)$. This exercises boundary pathology; the prescribed remedy is boundary reparametrisation.

Arguments

This function takes no arguments.

Mathematics

On the constrained scale, the density is

$$ p(b, z) = \alpha, b^{\alpha-1} \cdot \prod_{j=1}^{d-1}\phi(z_j), \qquad b \in (0,1),; z_j \in \mathbb{R}, $$

where $\phi$ is the standard normal density and $d = 4$. Transforming to the unconstrained scale $\phi_b = \mathrm{logit}(b)$, i.e. $b = \mathrm{logit}^{-1}(\phi_b) = 1/(1+e^{-\phi_b})$, the Jacobian is $db/d\phi_b = b(1-b)$. The unconstrained log-density (dropping the constant $\log\alpha$) is

$$ \log p(\phi_b, z) = \alpha \log b + \log(1 - b) - \frac{1}{2}\sum_{j=1}^{d-1} z_j^2, $$

where $\log(1-b)$ combines the Jacobian term $\log b + \log(1-b)$ with the $(\alpha-1)\log b$ from the Beta density, yielding $\alpha\log b + \log(1-b)$.

The gradient with respect to $\phi_b$ is

$$ \frac{\partial \log p}{\partial \phi_b} = \alpha - (\alpha+1),b, $$

and with respect to each $z_j$ is $-z_j$, so

$$ \nabla \log p = \bigl(\alpha - (\alpha+1)b,; -z_1,; \ldots,; -z_{d-1}\bigr). $$

The Stan code declares d and alpha as data, b as a real<lower=0, upper=1> parameter, z as a vector[d-1] parameter, with b ~ beta(alpha, 1) and z ~ std_normal().

Returns

The return value of .gdpar_geom_target(...), i.e. a gdpar_geometry_target object. The embedded make closure (with signature make(n = 1, difficulty = 0.3)) returns a list containing:

Element Type / Structure
stan_code Character scalar (Stan model string)
stan_data list(d = 4L, alpha = <difficulty>)
log_prob Function: function(theta) → numeric scalar
grad_log_prob Function: function(theta) → numeric vector of length d
dim Integer scalar 4L
param_names Character vector c("phi_b", "z[1]", "z[2]", "z[3]")

The registered metadata is: id = "G6_boundary", label = "Mass pinned against a (0, 1) bound", pathology = "boundary", geometry_remedy = "boundary_reparam", culprit = "b", difficulty_scales_with_n = FALSE, bounds = list(b = c(0, 1)), default_n = 1, default_difficulty = 0.3, n_grid = c(1, 4, 16).

Notes

  • difficulty is coerced to numeric and used as alpha; values below 1 (default 0.3) concentrate Beta mass near $b=0$.
  • The inv_logit helper is defined inside make and closes over nothing external.
  • log1p(-b) is used for numerical stability in computing $\log(1-b)$.
  • The R-side log_prob drops the constant $\log\alpha$; the gradient is exact for the unnormalised log-density.
  • n is accepted but unused (structural geometry, difficulty_scales_with_n = FALSE).
  • bounds is set to list(b = c(0, 1)), recording the constrained range of the culprit parameter.
  • culprit is "b" (the constrained parameter), not "phi_b" (the unconstrained working parameter).

.gdpar_geom_G7()

Purpose

Factory function that constructs and registers the "G7_flat_direction" geometry target — a 2-dimensional target where the sum $a+b$ is identified by a unit-variance normal likelihood centred at 0, but the contrast $a-b$ is only weakly identified through wide priors $a, b \sim \mathcal{N}(0, \tau^2)$ with large $\tau$. This produces a near-flat direction (small Hessian eigenvalue); the prescribed remedy is reparametrisation/elimination of the redundant parameter.

Arguments

This function takes no arguments.

Mathematics

The log-density (up to additive constants) is

$$ \log p(a, b) = -\frac{1}{2}(a+b)^2 - \frac{1}{2\tau^2}(a^2 + b^2), $$

where the first term is the log-likelihood $\mathcal{N}(0 \mid a+b, 1)$ and the second term is the sum of the two prior log-densities $\mathcal{N}(a \mid 0, \tau^2) + \mathcal{N}(b \mid 0, \tau^2)$.

The gradient is

$$ \nabla \log p = \begin{pmatrix} -(a+b) - \tau^{-2},a \[4pt] -(a+b) - \tau^{-2},b \end{pmatrix}. $$

The Hessian is

$$ H = \begin{pmatrix} -1 - \tau^{-2} & -1 \ -1 & -1 - \tau^{-2} \end{pmatrix}, $$

whose eigenvalues are $-\tau^{-2}$ (eigenvector $(1,-1)/\sqrt{2}$, the flat contrast direction) and $-2 - \tau^{-2}$ (eigenvector $(1,1)/\sqrt{2}$, the identified sum direction). As $\tau \to \infty$, the smaller eigenvalue $\to 0$, making the contrast direction increasingly flat.

The Stan code declares tau as data, a and b as real parameters, with target += normal_lpdf(0 | a + b, 1), a ~ normal(0, tau), b ~ normal(0, tau).

Returns

The return value of .gdpar_geom_target(...), i.e. a gdpar_geometry_target object. The embedded make closure (with signature make(n = 1, difficulty = 100)) returns a list containing:

Element Type / Structure
stan_code Character scalar (Stan model string)
stan_data list(tau = <difficulty>)
log_prob Function: function(theta) → numeric scalar
grad_log_prob Function: function(theta) → numeric vector of length 2
dim Integer scalar 2L
param_names Character vector c("a", "b")

The registered metadata is: id = "G7_flat_direction", label = "Reparametrisation redundancy (flat a - b)", pathology = "flat_direction", geometry_remedy = "reparam_eliminate", culprit = c("a", "b"), difficulty_scales_with_n = FALSE, bounds = NULL, default_n = 1, default_difficulty = 100, n_grid = c(1, 4, 16).

Notes

  • difficulty is coerced to numeric and used as tau; larger values widen the priors and flatten the contrast direction.
  • inv_tau2 <- 1 / tau^2 is precomputed once in make and captured by the log_prob and grad_log_prob closures for efficiency.
  • The R-side log_prob omits normalising constants from both the likelihood and the priors; the gradient is exact for the unnormalised log-density.
  • n is accepted but unused (structural geometry, difficulty_scales_with_n = FALSE).
  • culprit is a length-2 character vector c("a", "b") because both parameters jointly contribute to the redundancy.
  • bounds is NULL (both parameters are unconstrained).

print.gdpar_geometry_target(x, ...)

Purpose

S3 print method for objects of class gdpar_geometry_target. Produces a human-readable multi-line summary of the geometry target's identity, pathology, remedy, culprit, and difficulty-scaling flag.

Arguments

Argument Type Meaning
x gdpar_geometry_target The geometry target object to print. Expected to contain elements id, label, pathology, geometry_remedy, culprit, difficulty_scales_with_n.
... Any Unused; present for S3 generic compatibility.

Returns

Invisibly returns x (via invisible(x)).

Notes

  • Dispatches on class gdpar_geometry_target via the S3 print generic.
  • The culprit field is checked with all(is.na(x$culprit)); if all entries are NA, the string "none" is printed. Otherwise, the culprit values are collapsed with paste(x$culprit, collapse = ", ").
  • Output format (four lines):
    1. <gdpar_geometry_target> <id>
    2. <label>
    3. pathology: <pathology> | remedy: <geometry_remedy>
    4. culprit: <culprits or "none"> | difficulty scales with n: <difficulty_scales_with_n>
  • No validation is performed on x; missing elements would trigger R's standard subsetting error.
  • Exported (has @export roxygen tag).

← Part IV — Exhaustive Function Reference (4/7) · gdpar Wiki Home · Part IV — Exhaustive Function Reference (6/7) →

gdpar Wiki

Home

Clone this wiki locally