User-defined GMRF latent blocks with tgmrf()

What `tgmrf()` is

tgmrf() lets a user script define a Gaussian Markov Random Field (GMRF) latent block and plug it into tulpa’s inference layers as a first-class latent term. It is the latent-side dual of LikelihoodSpec: that interface lets a downstream package own the observation model; tgmrf() lets a single R script own one latent block.

The contract is two closures and an init vector:

my_block <- tgmrf(
  Q     = function(theta) { ... returns dgCMatrix ... },
  prior = function(theta) { ... returns scalar log-density ... },
  init  = c(...)
)

The closures are called numerically: no autodiff, no DSL, no codegen. For a Gaussian latent block

\[ \log p(z \mid \theta) \;=\; \tfrac{1}{2}\,\log\det Q(\theta) - \tfrac{1}{2}(z-\mu)^\top Q(\theta)(z-\mu) + \text{const} \]

the score and Hessian in $z$ are closed-form:

\[ \partial_z \log p(z \mid \theta) = -Q(\theta)(z-\mu), \qquad \partial_z^{\,2} \log p(z \mid \theta) = -Q(\theta). \]

So tulpa’s Laplace inner step only needs Q(theta) and mu(theta) at numeric theta. Hyperparameter gradients (d/d theta) come from finite differences, typically 2-5 extra Q evaluations per outer step, which is cheap when Q is sparse. The same block reaches every inference tier: Laplace, IMH-Laplace, NUTS, VI.

How a block becomes a fitted model: the algorithm

A tgmrf object carries everything the engine needs and nothing it does not. The construct is a templated GMRF latent block: templated because the same machinery integrates any precision family the closures describe, latent because the block contributes a vector $z$ that lives under the linear predictor rather than in the likelihood. When you call tgmrf(), the constructor evaluates Q(init) exactly once. That single call does three jobs at registration time. It catches user errors early (a Q that throws, returns a dense matrix, or returns a non-square or asymmetric matrix all stop here with a specific message). It infers n_latent from nrow(Q(init)). And it captures the sparsity pattern, the set of stored nonzero positions, which the inner solver reuses at every later theta so the symbolic Cholesky factorisation is computed once, not per grid point.

The block records theta_dim and theta_names from init (named elements become the canonical hyperparameter names; unnamed ones get theta_1, theta_2, …), validates the optional mu, graph, bounds, and obs_idx arguments against the inferred shape, and stores the two closures by reference. Nothing is differentiated, compiled, or parsed. The returned object has class c("tgmrf", "tulpa_latent_block") and a type = "tgmrf" registry key that the nested-Laplace dispatcher recognises.

At fit time the integration follows a fixed sequence. The steps below describe what tulpa_nested_laplace() does with a tgmrf block; the IMH, VI, and NUTS adapters share the inner machinery and differ only in how they explore theta.

Build the outer grid over theta. If you supplied bounds, the driver lays a per-axis grid between lower and upper. The default resolution is 5 points per axis; the IMH, VI, and NUTS adapters expose this as pilot_axis_points. With $d$ hyperparameters the default grid is $5^d$ cells, so a two-parameter block is a 25-cell grid. You can also hand the driver an explicit grid matrix.
For each grid point theta_k, run the inner Laplace solve. The driver assembles the joint mode of $(\beta, z)$ given theta_k by Newton iteration. Each Newton step needs the score and Hessian in $(\beta, z)$. The $z$-part comes from Q(theta_k) and mu(theta_k) in closed form; the likelihood part comes from the chosen family (or a model-supplied likelihood). Because the sparsity pattern was captured at registration, the Cholesky of the system reuses the same symbolic factor at every grid point.
Score the grid point. At the inner mode the driver evaluates the Laplace approximation to the log marginal likelihood $\log p(y, \hat{z}_k \mid \theta_k)$ and adds the user prior(theta_k). This single scalar per grid point is what $log_marginal holds.
Integrate over the grid. The per-cell log-marginals become normalised integration weights ($weights), and the driver reports posterior moments of every hyperparameter as weighted means and SDs ($theta_mean, $theta_sd) plus weighted quantiles.
Diagnose the outer approximation. With the grid integrated, the driver fits a Gaussian proposal to the hyperparameter posterior and importance-samples it to compute the outer Pareto-$\hat{k}$ ($pareto_k); see the convergence section.

The adapters change only step 1’s exploration of theta. IMH localises the grid argmax, builds a finite-difference Hessian there, and proposes $\theta' \sim N(\text{mode}, \text{scale}^2 H^{-1})$ with an exact Metropolis accept/reject against log_marginal(theta) + prior(theta). VI runs L-BFGS on log_marginal(theta) and fits a Gaussian at the optimum. NUTS runs leapfrog integration on the same target with a finite-difference gradient. Every one of them calls the same inner Laplace solve to evaluate a candidate theta. The inner solve is the shared kernel, the outer loop is the only thing that varies. That is why a block written once reaches all four tiers without modification: the closures describe the latent prior, and the engine owns everything from the Newton system up.

A worked example: periodic AR1

The smallest non-trivial GMRF tulpa does not ship is a periodic AR1 (useful for diurnal, seasonal, or phase-wrapped data). The precision is tridiagonal with wrap-around:

\[ Q_{ii} \;=\; \tfrac{1+\rho^2}{\sigma^2}, \qquad Q_{ij} \;=\; -\,\tfrac{\rho}{\sigma^2}\;\text{ for } |i-j| = 1 \;(\mathrm{mod}\,n). \]

Reparameterise to unconstrained theta = (log_sigma, atanh_rho):

periodic_ar1 <- function(n) {
  tgmrf(
    Q = function(theta) {
      sigma <- exp(theta[1]); rho <- tanh(theta[2])
      d <- rep((1 + rho^2) / sigma^2, n)
      o <- rep(-rho / sigma^2, n)
      M <- Matrix::bandSparse(n, k = c(-1, 0, 1),
                              diagonals = list(o, d, o))
      M[1, n] <- M[n, 1] <- -rho / sigma^2     # wrap
      methods::as(methods::as(M, "generalMatrix"), "CsparseMatrix")
    },
    prior = function(theta) {
      dnorm(theta[1], 0, 1, log = TRUE) +     # weak on log_sigma
      dnorm(theta[2], 0, 1, log = TRUE)       # weak on atanh_rho
    },
    init   = c(log_sigma = 0, atanh_rho = atanh(0.3)),
    bounds = list(lower = c(log(0.3), atanh(0.0)),
                  upper = c(log(3.0), atanh(0.9))),
    name   = "periodic_ar1"
  )
}

The optional bounds argument sets the outer grid for nested-Laplace and the hard-wall region for IMH and NUTS. The prior closure is a log-density on theta. If you don’t have prior beliefs, function(theta) 0 is fine; it’s a flat prior on the unconstrained scale.

The closure shape is the contract. Q takes one numeric vector theta and returns a square symmetric sparse matrix; prior takes the same theta and returns one finite scalar; init is the numeric vector both closures expect. The reparameterisation to log_sigma and atanh_rho is deliberate: it maps the natural constraints $\sigma > 0$ and $|\rho| < 1$ to the whole real line, so the outer grid, the L-BFGS search, and the leapfrog steps all operate on unconstrained coordinates and never propose an invalid precision. The bounds you give are in those unconstrained units too. Use names on init and they propagate through every fit as the reported parameter names.

Simulate from the block

n <- 40L
theta_true <- c(log(0.8), atanh(0.6))
blk        <- periodic_ar1(n)

# z | theta_true ~ N(0, Q(theta_true)^{-1})
Q_true <- blk$Q(theta_true)
L      <- Matrix::Cholesky(Q_true)
z_true <- as.numeric(Matrix::solve(L, rnorm(n), system = "Lt"))

# y_i ~ Poisson(exp(beta_0 + z_i))
beta0 <- 0.3
y     <- rpois(n, exp(beta0 + z_true))
X     <- matrix(1, n, 1L)

Tier 2 (structured): nested Laplace

fit_lap <- tulpa_nested_laplace(
  y = y, n_trials = rep(1L, n), X = X,
  prior  = blk,
  family = "poisson"
)
fit_lap$theta_mean    # posterior mean over the outer grid

Tier 2 (structured): Pathfinder VI

tulpa_tgmrf_vi() runs L-BFGS on log_marginal(theta) and fits a Gaussian at the optimum. Same tier as the nested-Laplace posterior but gives you a full covariance matrix rather than a discrete grid.

fit_vi <- tulpa_tgmrf(
  y = y, n_trials = rep(1L, n), X = X, block = blk,
  family  = "poisson", mode = "vi",
  n_draws = 800L
)
fit_vi$mode_theta
fit_vi$elbo

Tier 1 (exact): IMH-Laplace

The Tier-1 upgrade composes the Laplace body with an independence-MH bias correction. The VI fit makes a near-optimal proposal: high acceptance, low MC error.

fit_imh <- tulpa_tgmrf(
  y = y, n_trials = rep(1L, n), X = X, block = blk,
  family = "poisson", mode = "imh",
  n_iter = 1200L, warmup = 400L
)
fit_imh$means
fit_imh$mean_accept

Tier 1 (exact): NUTS over `theta`

For higher-dim theta (5+), the gradient-based sibling tulpa_tgmrf_nuts() mixes better than independence-MH. For theta_dim = 2 IMH is generally faster; NUTS is the right call when the posterior is multi-modal or strongly non-Gaussian.

fit_nuts <- tulpa_tgmrf(
  y = y, n_trials = rep(1L, n), X = X, block = blk,
  family    = "poisson", mode = "nuts",
  n_iter    = 200L, warmup = 100L,
  max_depth = 4L
)
fit_nuts$means

Comparing posteriors

The three estimators target the same posterior p(theta | y):

data.frame(
  parameter = blk$theta_names,
  true      = theta_true,
  laplace   = fit_lap$theta_mean,
  vi_mode   = fit_vi$mode_theta,
  imh_mean  = fit_imh$means
)

Laplace and VI agree to within grid resolution / Gaussian-fit curvature; IMH lifts the structured-tier bias at the cost of n_iter inner Laplace solves.

The mathematics

Fix notation first. The model has three layers. The data is $y = (y_1, \dots, y_N)$ with a per-observation likelihood from family. A fixed-effects design matrix $X$ of dimension $N \times p$ carries coefficients $\beta \in \mathbb{R}^p$. The latent block contributes a vector $z \in \mathbb{R}^{n_{\mathrm{lat}}}$ where $n_{\mathrm{lat}}$ is n_latent. An index vector obs_idx maps each observation to a latent slot; when it is NULL the driver assumes $N = n_{\mathrm{lat}}$ and uses row order. The linear predictor for observation $i$ is \[ \eta_i \;=\; (X\beta)_i \;+\; z_{\,\mathrm{obs\_idx}[i]}, \] and the mean is $g^{-1}(\eta_i)$ for the family’s link $g$. The hyperparameter vector is $\theta \in \mathbb{R}^d$ with $d$ equal to theta_dim, living on the unconstrained scale your closures define.

The GMRF prior

A Gaussian Markov Random Field is a multivariate normal whose precision matrix $Q$ is sparse, and the sparsity is the conditional-independence structure: $z_i$ and $z_j$ are conditionally independent given the rest exactly when $Q_{ij} = 0$. The block density is \[ p(z \mid \theta) \;=\; (2\pi)^{-n_{\mathrm{lat}}/2}\, |Q(\theta)|^{1/2}\, \exp\!\Big(-\tfrac{1}{2}(z-\mu(\theta))^\top Q(\theta)(z-\mu(\theta))\Big), \] where $\mu(\theta)$ is the mean vector (your mu closure, or zero when mu is NULL) and $|Q(\theta)|$ is the determinant. Taking logs, \[ \log p(z \mid \theta) \;=\; \tfrac{1}{2}\log|Q(\theta)| - \tfrac{1}{2}(z-\mu)^\top Q(\theta)(z-\mu) + \text{const}. \] The const absorbs the $(2\pi)$ factor, which does not depend on $z$ or $\theta$ and drops out of every gradient. Differentiating in $z$ gives the closed-form score and Hessian quoted above: $\partial_z \log p = -Q(z-\mu)$ and $\partial_z^2 \log p = -Q$. This is why you never write gradient code for the latent block. The Hessian in $z$ is just $-Q(\theta)$, a constant in $z$, so the latent part of the Newton system is exact, not an approximation.

The precision parameterisation

The closure Q(theta) is the precision family. Whatever structure you encode (AR1, periodic AR1, a random walk, a custom spatial GMRF) lives entirely in how Q maps $\theta$ to nonzero entries. The reparameterise step matters here. The natural parameters of a precision are often positive scales or bounded correlations, and the integrator works best on the whole real line, so the convention is to define theta on an unconstrained scale and transform inside Q. For the periodic AR1 above, $\theta = (\log\sigma, \operatorname{atanh}\rho)$ so that $\sigma > 0$ and $|\rho| < 1$ hold automatically. The Jacobian of that transform is absorbed into your prior(theta): because the prior is a log-density on the unconstrained $\theta$, you supply whatever density you want there directly and the engine does not add a change-of-variables term. A flat function(theta) 0 is a flat prior on the unconstrained scale, which is weakly informative on the natural scale.

The inner Laplace solve

For a fixed $\theta$ the engine needs the marginal likelihood $p(y \mid \theta) = \int p(y \mid \beta, z)\,p(\beta)\,p(z \mid \theta)\, \mathrm{d}\beta\,\mathrm{d}z$. The Laplace approximation replaces the integrand by a Gaussian centred at its mode. Let $w = (\beta, z)$ and $\ell(w) = \log p(y \mid w) + \log p(z \mid \theta) + \log p(\beta)$ be the log joint. Newton iteration finds the mode $\hat{w}$ by repeatedly solving \[ H(w)\,\Delta w \;=\; -\,\nabla \ell(w), \qquad H(w) \;=\; -\,\nabla^2 \ell(w), \] where $H$ is the negative Hessian of the log joint. The $z$-block of $H$ is $Q(\theta)$ plus the likelihood’s Fisher information, the $\beta$-block is $X^\top W X$ for the family’s working weights $W$, and the cross term couples them through obs_idx. The Laplace approximation to the marginal is \[ \log p(y \mid \theta) \;\approx\; \ell(\hat{w}) \;+\; \tfrac{1}{2}\,d_w\log(2\pi) \;-\; \tfrac{1}{2}\log|H(\hat{w})|, \] with $d_w = p + n_{\mathrm{lat}}$ the joint dimension. The engine adds your prior(theta) to this to form log_marginal(theta) + prior(theta), which is the unnormalised log posterior of $\theta$ that every outer layer targets. The sparsity pattern captured at registration is what makes the $\log|H|$ and the linear solve cheap: a sparse Cholesky on a pattern that does not change across grid points reuses one symbolic factor.

Outer integration over `theta`

The nested-Laplace layer treats $\theta$ as the remaining integration problem. It evaluates $L(\theta_k) = \log p(y \mid \theta_k) + \log p(\theta_k)$ at each grid point $\theta_k$, normalises to weights $w_k \propto \exp(L(\theta_k))$ with $\sum_k w_k = 1$, and reports posterior summaries as weighted statistics over the grid. The posterior mean of hyperparameter $j$ is $\sum_k w_k\,\theta_{k,j}$ and the SD comes from the weighted second moment. Derived quantities such as the natural-scale $\sigma = e^{\theta_1}$ or $\rho = \tanh(\theta_2)$ are summarised by evaluating the transform at each grid cell and taking weighted quantiles, not by transforming the mean. The IMH, VI, and NUTS layers replace the fixed grid by an adaptive exploration of the same $L(\theta)$ surface but integrate against the identical inner Laplace solve. Every symbol here is something your two closures and the family definition fully determine.

Convergence

There are two convergence questions and they belong to different layers. The first is whether the inner Newton solve found the joint mode of $(\beta, z)$ at each $\theta$. The second is whether the outer exploration of $\theta$ (the grid, the L-BFGS optimum, or the MCMC chain) has actually characterised the hyperparameter posterior. These are diagnosed separately.

For the inner solve, tulpa_nested_laplace() returns $n_iter, the number of Newton iterations spent per grid point, and $log_marginal, the per-cell scores. The inner budget is control$max_iter (default 50L) with tolerance control$tol (default 1e-6). If a cell hits max_iter without converging, or if log_marginal is non-finite at any cell, the fit is unreliable at that $\theta$: tighten bounds so the grid stays where $Q(\theta)$ is well-conditioned, or raise max_iter. The IMH and NUTS adapters add an eager check: they evaluate the inner solve at the pilot grid argmax before sampling and stop with a specific message if it is non-finite, because a failure there is a bug in the closures rather than numerical infeasibility at the edge.

For the outer approximation on a structured-tier fit, the relevant number is the Pareto-$\hat{k}$, the iid-fit counterpart of Rhat. The nested-Laplace fit reports $pareto_k and $pareto_k_is_ess: the integrator fits a Gaussian proposal to the hyperparameter posterior, importance-samples it (control$k_samples, default 200L, each one extra inner solve), and fits a generalised Pareto tail to the importance ratios. A $\hat{k} < 0.7$ means the Gaussian grid characterises the hyperparameter posterior well and the nested integration is reliable. A $\hat{k} \geq 0.7$ means the posterior is too skewed or heavy-tailed for the Gaussian grid, and you should escalate to the exact tier (tulpa_tgmrf_imh() or tulpa_tgmrf_nuts()). The diagnostic is computed for a single-block, single positive-scale-axis grid and left NA (with the grid’s quadrature effective sample size as the fallback) for multi-axis or bounded grids. diagnostic_summary(fit) surfaces this for any non-chain fit, and tulpa_psis() is the underlying PSIS core if you want to compute it on your own importance ratios.

For the exact-tier MCMC fits the diagnostics are the usual ones. The IMH fit reports $mean_accept, the post-warmup acceptance rate; a rate in a healthy band (the VI fit makes a near-optimal proposal, so a good model sees high acceptance) means the proposal matches the target. A collapsed acceptance rate is the signal to switch to NUTS. The NUTS fit reports $mean_accept plus $tree_depth per draw and the adapted $epsilon; draws repeatedly hitting max_depth indicate a step size that is too small for the posterior geometry. Both produce a $draws matrix you can feed to mcmc_diagnostics(fit) for Rhat and effective sample size, or to check_diagnostics(fit) for a pass/fail summary. The cleanest single check is to run the structured-tier fit, read its Pareto-$\hat{k}$, and only reach for the MCMC tiers when the diagnostic says the Gaussian approximation is biased.

Tuning

The split between statistical arguments and tuning knobs follows the engine convention. tulpa_nested_laplace() carries statistical arguments (y, n_trials, X, prior, family, phi) at the top level and puts every numerical or performance knob in a single control = list(). The tgmrf adapters expose their tuning as named arguments because each adapter has a small, specific set.

The first decision is the tier. Start with tulpa_nested_laplace() for a fast structured fit and read its $pareto_k. If the diagnostic is below 0.7 the structured posterior is trustworthy and you are done; if it is above 0.7, escalate. tulpa_tgmrf_imh() is the cheapest exact tier and the right default escalation for low theta_dim (2-4) with a roughly Gaussian posterior, since its proposal is the VI Gaussian, so acceptance stays high. tulpa_tgmrf_nuts() is for higher theta_dim (5+) or a strongly non-Gaussian or multi-modal posterior, where independence-MH acceptance collapses and gradient information pays for itself.

Inside the structured fit, the knobs that matter are the inner Newton budget and the diagnostic cost. control$max_iter (default 50L) and control$tol (default 1e-6) govern the inner solve; raise max_iter only if cells fail to converge. control$n_threads (default 1L) sets OpenMP threads inside each inner solve, useful when n_latent is large. control$diagnose_k (default TRUE) and control$k_samples (default 200L) control the Pareto-$\hat{k}$ computation; set diagnose_k = FALSE to skip the extra inner solves when you do not need the accuracy gate. The bounds you set on the block itself define the grid extent and the hard wall for IMH and NUTS. Set them generously enough to contain the posterior mass but tight enough to exclude $\theta$ where $Q$ goes singular.

The adapters share pilot_axis_points (default 5L), the per-axis resolution of the pilot grid that initialises every adapter and feeds the mode and Hessian estimates. Raise it when the default grid is too coarse to locate the mode. tulpa_tgmrf_imh() adds n_iter (default 2000L), warmup (default n_iter %/% 2L), thin, scale (the proposal scale multiplier, default 1.0, lower it if acceptance is too high, raise it if too low), and fd_step (default 0.05, the finite-difference step for the mode Hessian). tulpa_tgmrf_vi() adds n_draws (default 1000L, the variational draws returned), max_lbfgs (default 100L), and lbfgs_tol (default 1e-6). tulpa_tgmrf_nuts() adds epsilon (default 0.2 / sqrt(theta_dim), the initial leapfrog step), max_depth (default 6L, the tree-depth cap so each draw is at most $2^6$ leapfrog steps), target_accept (default 0.65 for the warmup step-size adaptation), and fd_gradient_step (default 0.02, the central-difference step for the gradient). The defaults are tuned for the common theta_dim <= 5 case; reach for the knobs when the diagnostics in the previous section flag a specific failure, not pre-emptively.

Which tier when

Tier	Function	Cost	When to reach for it
2 (structured)	`tulpa_nested_laplace()`	`n_grid` inner Laplace solves	Quick fit, diagnostic, debias-target
2 (structured)	`tulpa_tgmrf_vi()`	L-BFGS + `n_draws` extra evals	Want covariance instead of grid; IMH proposal
1 (exact)	`tulpa_tgmrf_imh()`	`n_iter` Laplace solves	Low theta dim, near-Gaussian posterior
1 (exact)	`tulpa_tgmrf_nuts()`	`n_iter * 2^depth * 2 * theta_dim` Laplace solves	High theta dim or non-Gaussian

The general rule: start with tulpa_nested_laplace() for a sanity check, run tulpa_tgmrf_imh() for exact-tier posterior moments, and reach for tulpa_tgmrf_nuts() only when IMH acceptance collapses.

How this compares to other latent-block interfaces

Several frameworks let a user inject a custom latent structure, and they trade off along the same axis: how much of the structure you write versus which inference tiers you reach.

INLA’s rgeneric interface lets you define a GMRF in R by supplying functions for the precision and a few related quantities. The callback runs in R and there is no autodiff, so the block is restricted to INLA’s Laplace machinery: the structured tier only, with no path to exact MCMC and no debias step on a non-Gaussian residual. tgmrf() matches the ergonomics (two R closures, no compilation) but the same block reaches every tulpa tier, including the exact IMH and NUTS samplers, because the outer layers all drive the same inner Laplace solve.

INLA’s cgeneric interface moves the precision factory to a C function, which is faster than rgeneric but still has no autodiff and still supports the Laplace tier only: there is no exact-MCMC support for a cgeneric block. The tgmrf_cpp() constructor (noted below) is the closer analog on the speed axis, and it keeps the full-tier reach that cgeneric lacks.

Stan takes the opposite approach: you write the entire model in its DSL, which is then parsed and compiled to autodiffable C++. That gives exact HMC over everything, but it asks you to specify the whole model: there is no notion of plugging one latent block into an engine that owns the rest, and you pay full HMC cost on every block including the Gaussian-latent ones a Laplace approximation would handle cheaply.

TMB is the closest analog: you write a templated C++ snippet and get autodiff through CppAD, with a Laplace approximation over the random effects. The tgmrf() R path is lighter (no C++, closed-form latent score from $Q$ rather than autodiff) and the tgmrf_cpp() path is the direct counterpart for hot inner loops. The distinguishing move across all four comparisons is tulpa’s nested-approximation-plus-debias design: the structured tier gives a cheap Laplace answer, and the same block escalates to an exact-tier correction when the diagnostic says the approximation is biased, without rewriting anything.

What `tgmrf()` does not require

No gradient code. The score in z is closed-form from Q.
No Hessian code. The Hessian in z is just -Q.
No autodiff. Hyperparameter gradients are finite-difference.
No C++. The R closure path reaches every tier.

The tgmrf_cpp() constructor offers a templated-C++ backend with the same S3 contract for cases where the R-to-C boundary on Q(theta) becomes the bottleneck (large n_latent, hot CCD grids). It takes a .cpp file that defines the Q, mu, and log-prior kernels as templated C++ and registers them through the TULPA_REGISTER_TGMRF macro; inst/examples/tgmrf_periodic_ar1.cpp is a worked example. The user-visible API does not change; the same downstream consumers (formula parser, inference layers, methods) treat the two paths identically.