tgmrf: User-defined GMRF latent block

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

User-defined GMRF latent block

Description

tgmrf() lets a user define a Gaussian Markov random field as a latent block in a tulpa model by supplying two R closures: a precision-matrix factory Q(theta) and a log-prior on theta prior(theta). The resulting object plugs into a formula via latent(tgmrf(...)) and is consumed by every tulpa inference tier (Laplace, EM+Laplace, VI, nested_laplace+CCD, NUTS, IMH-Laplace).

This is the latent-side dual of LikelihoodSpec. LikelihoodSpec lets a model package own the observation model; tgmrf() lets a script own a single latent block.

Usage

tgmrf(
  Q,
  prior,
  init,
  mu = NULL,
  graph = NULL,
  bounds = NULL,
  obs_idx = NULL,
  name = NULL
)

## S3 method for class 'tgmrf'
print(x, ...)

## S3 method for class 'tgmrf'
summary(object, ...)

Arguments

Q	A function ⁠function(theta)⁠ returning a sparse precision matrix of class Matrix::sparseMatrix (typically dgCMatrix). The matrix must be square and symmetric. Called once at registration with theta = init to infer n_latent and capture the sparsity pattern.
prior	A function ⁠function(theta)⁠ returning a finite numeric scalar – the log-prior density at theta.
init	Numeric vector of starting values for theta. Names, if present, become the canonical theta names.
mu	Optional function ⁠function(theta)⁠ returning a numeric vector of length n_latent. Default NULL is equivalent to a zero mean.
graph	Optional sparse matrix specifying the upper bound on Q's sparsity pattern. If supplied, the registration check verifies that the nonzero pattern of Q(init) is a subset.
bounds	Optional list with components lower and upper, each a numeric vector of length length(init). Used by nested_laplace + CCD to build the outer grid; unused by NUTS / VI.
obs_idx	Optional integer vector mapping each observation to a latent slot in ⁠[1, n_latent]⁠. If NULL (default), the fit-time driver assumes N == n_latent and uses seq_len(N) – i.e. one observation per latent slot, in row order.
name	Optional character; cosmetic label used by print() / summary().
x	A tgmrf block.
...	Passed from summary() to print().
object	A tgmrf block.

Details

For a Gaussian latent block z ~ N(mu(theta), Q(theta)^{-1}):

\log p(z\mid\theta) = \tfrac{1}{2}\log\det Q(\theta) - \tfrac{1}{2}(z - \mu(\theta))^\top Q(\theta)(z - \mu(\theta)) + \mathrm{const}

the gradient -Q(z - mu) and Hessian -Q wrt z are closed form, so the user never writes gradient code. The Laplace inner step needs only Q(theta) and mu(theta) at numeric theta. NUTS and VI additionally use a forward finite-difference gradient on theta, costing dim(theta) extra Q calls per outer step – cheap for the typical dim(theta) <= 5.

Value

An object of class c("tgmrf", "tulpa_latent_block") with components:

• Q, prior, mu – the user closures (or NULL for mu).

• init, theta_names, theta_dim – hyperparameter metadata.

• n_latent – inferred from Q(init).

• pattern – captured sparsity pattern as a dgCMatrix of 1s.

• graph – user-supplied pattern, validated against pattern if given.

• bounds, name – passed through.

summary() prints the block, its log prior at init and its bounds, and returns object invisibly.

See Also

latent() for the formula slot that consumes a tgmrf block.

Examples

# Periodic AR(1) block, wrap-around tridiagonal precision.
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
      methods::as(M, "CsparseMatrix")
    },
    prior = function(theta) {
      stats::dnorm(theta[1], 0, 1, log = TRUE) +
        stats::dnorm(theta[2], 0, 1, log = TRUE)
    },
    init = c(log_sigma = 0, atanh_rho = 0),
    name = "periodic_ar1"
  )
}

blk <- periodic_ar1(20)
print(blk)

tulpa documentation built on Oct. 19, 2026, 5:09 p.m.