GGM

Implementation details of the Gaussian graphical model in bgms. Source: src/models/ggm/.

The GGM is parameterized by the $p \times p$ precision matrix $\boldsymbol{\Theta} = \boldsymbol{\Sigma}^{-1}$. The off-diagonal entries $\theta_{ij}$ encode partial associations: $\theta_{ij} = 0$ if and only if variables $i$ and $j$ are conditionally independent given all other variables. The model maintains an upper-triangular Cholesky factor $\boldsymbol{\Phi}$ such that $\boldsymbol{\Phi}^\top \boldsymbol{\Phi} = \boldsymbol{\Theta}$ throughout sampling. For complete data, the likelihood depends on the data only through the sufficient statistic $\mathbf{S} = \mathbf{X}^\top \mathbf{X}$.

Priors are specified on the partial-association scale: the slab (interaction_prior) acts on $-\theta_{ij}/2$ and the scale prior (precision_scale_prior) on $\theta_{ii}/2$, matching the convention of sample_ggm_prior(). Both are polymorphic prior objects (see Parameter Priors); the defaults are a Normal slab and an exponential diagonal.

Sampling overview

The GGM supports three update methods:

NUTS (default) — Gradient-based sampling of all continuous parameters simultaneously. Uses the constrained Cholesky parameterization, whose null-space coordinates keep excluded edges exactly zero under edge selection.
Adaptive Metropolis — Element-wise updates of individual $\theta_{ij}$ entries via rank-1 Cholesky updates.
Gibbs — An exact conjugate row-block sweep; with edge selection, the between-model step switches to the conjugate edge birth/death proposal.

The sampler is selected by the update_method argument in bgm(). Edge indicators $\gamma_{ij}$ are updated by Metropolis–Hastings moves regardless of the continuous-parameter sampler. Within each iteration, edge indicator updates run first, then the continuous sampler runs with the graph held fixed. See Sampler Hierarchy for the full dispatch logic.

Free-element Cholesky parameterization

NUTS requires an unconstrained parameter vector. The GGM achieves this through a column-wise parameterization of $\boldsymbol{\Phi}$ that separates free parameters from constrained (zero-precision) entries.

Column layout

The parameter vector $\theta$ is packed column by column:

\[ \theta = (\psi_1,\; f_2, \psi_2,\; f_3, \psi_3,\; \ldots,\; f_p, \psi_p), \]

Notation, fixed before anything else, because three symbols carry this page:

$\theta$ with no subscripts is the packed sampler position, the one vector the sampler moves through. The NUTS and Constrained Cholesky Parameterization pages call that same vector $q$; this page cannot, because $q$ is needed below as the Cholesky column index, mirroring the code.
$\theta_{ij}$, doubly subscripted, is always an entry of the precision matrix $\boldsymbol{\Theta}$, never a position coordinate.
$q$ is a column index of the Cholesky factor $\Phi$, nothing else.

For each column $q$:

$\psi_q = \log(\Phi_{qq})$ — the log-transformed diagonal entry, ensuring $\Phi_{qq} > 0$.
$f_q \in \mathbb{R}^{d_q}$ — free parameters for column $q$’s off-diagonal entries, with $d_q$ equal to the number of included edges in column $q$.

Null-space basis

When column $q$ has $m_q$ excluded edges ($\gamma_{iq} = 0$), the constraint matrix $A_q$ encodes the zero-precision requirements $\theta_{iq} = \sum_l \Phi_{li}\,\Phi_{lq} = 0$ (see Constrained Cholesky Parameterization). The off-diagonal entries of $\boldsymbol{\Phi}$ for column $q$ must lie in the null space of $A_q$. A Givens QR factorization of $A_q^\top$ produces an orthonormal null-space basis $N_q$, and the off-diagonal entries are recovered as $x_q = N_q f_q$.

This maps $d_q$ unconstrained parameters to the $(q - 1)$-dimensional off-diagonal slot while satisfying all $m_q$ constraints exactly. When $m_q = 0$ (no excluded edges in column $q$), $N_q$ is the identity and $f_q = x_q$ directly.

Active and full layouts

NUTS always integrates in the active coordinates $(f_q, \psi_q)$, of dimension $p + \lvert E\rvert$ — the null-space rotation $N_q f_q$ maps free parameters to Cholesky entries, so the sampler never sees the excluded dimensions, with or without edge selection. The full (zero-padded) layout of dimension $p + p(p-1)/2$ exists only as a fixed-slot frame: the mass matrix is estimated there so entries keep their positions when edge moves change the active set (see NUTS — Mass matrix), and sample storage uses the $p(p+1)/2$ upper triangle of $\boldsymbol{\Theta}$ itself.

`GraphConstraintStructure`

The GraphConstraintStructure in graph_constraint_structure.h precomputes the index maps needed by the parameterization. For each column $q$ it stores:

excluded_indices / included_indices — which rows are constrained vs. free
m_q / d_q — number of constraints / free parameters
theta_offsets[q] — start of column $q$’s block in the active $\theta$ vector
full_theta_offsets[q] — start of column $q$’s block in the full position vector

The structure is rebuilt whenever the edge indicators change (the constraint_dirty_ flag triggers a lazy rebuild at the next gradient evaluation).

Gradient computation

The gradient engine (GGMGradientEngine in ggm_gradient.h/.cpp) computes $\log p(\theta)$ and $\nabla_\theta \log p(\theta)$ jointly.

Forward map

The forward map $\theta \to \boldsymbol{\Phi} \to \boldsymbol{\Theta}$ processes columns left to right:

For each column $q$, build $A_q$ from earlier (finalized) columns.
QR-factorize $A_q^\top$ via Givens rotations to get $N_q$.
Compute $x_q = N_q f_q$ and set $\Phi_{qq} = e^{\psi_q}$.
Accumulate the log-Jacobian determinant from the QR factors.

The log-Jacobian has three terms:

\[ \log|J| = p \log 2 + 2\sum_q \psi_q + \sum_{q=0}^{p-2}(p - 1 - q)\,\psi_q - \sum_q \sum_j \log|R_{q,jj}|, \]

combining the $\boldsymbol{\Phi} \to \boldsymbol{\Theta}$ Jacobian, the $\psi_q \to \Phi_{qq}$ exponential transform, and the null-space QR rotation.

Backward pass

The log-posterior being differentiated is

\[ \log p = \tfrac{n}{2}\log|\boldsymbol{\Theta}| - \tfrac{1}{2}\operatorname{tr}(\boldsymbol{\Theta}\mathbf{S}) + \delta \log|\boldsymbol{\Theta}| + \log p_{\text{prior}}(\boldsymbol{\Theta}) + \log|J|, \]

where $\delta \log|\boldsymbol{\Theta}|$ is the determinant tilt and the prior term sums the polymorphic slab and scale priors over the included entries. The backward pass runs in two phases:

Phase 1 (data + priors → $\bar{\boldsymbol{\Phi}}$) — Computes the adjoint of the Cholesky factor from the data likelihood, the tilt, the interaction prior on off-diagonals, and the scale prior on diagonals; see Parameter Priors. Throughout, a bar denotes the reverse-mode adjoint: $\bar{A} = \partial \log p / \partial A$, the sensitivity of the log-posterior to that quantity.
Phase 2 (right-to-left column sweep) — Propagates $\bar{\boldsymbol{\Phi}}$ through the stored Givens rotations to extract $\partial\log p / \partial f_q$ and $\partial\log p / \partial \psi_q$. Reverse-mode differentiation through the QR factorization captures cross-column dependencies automatically.

Determinant tilt

The prior on $\boldsymbol{\Theta}$ is multiplied by $|\boldsymbol{\Theta}|^{\delta}$ with $\delta \ge 0$ (the delta argument of bgm(); NULL auto-resolves to $0.5 \log p$). The tilt softly repels the chain from the boundary of the positive-definite cone, where the null-space geometry degenerates. It enters every update path: the NUTS gradient adds $\delta$ times the log-determinant adjoint, the element-wise Metropolis and edge moves add $\delta \cdot \Delta \log|\boldsymbol{\Theta}|$ to their acceptance ratios (an $O(p)$ rank-2 determinant-lemma computation), and the row-block Gibbs absorbs it as a shift in the $\xi$ Gamma shape (see below). The tilt also reshapes the graph marginal of the joint spike-and-slab prior, which is what the edge-prior normalizing-constant corrections and the hierarchical specification account for.

Element-wise Metropolis

When the adaptive-Metropolis sampler is selected, each off-diagonal entry $\theta_{ij}$ is updated individually via Metropolis–Hastings. The algorithm exploits a structural property of the Cholesky parameterization: changing one off-diagonal entry $\theta_{ij}$ while holding the Cholesky factor consistent requires adjusting the diagonal $\theta_{jj}$ deterministically. This reduces a constrained positive-definite update to a one-dimensional proposal, and the entire update — proposal, likelihood ratio, and Cholesky maintenance — runs without permutation or refactorization of the precision matrix.

Roverato constraint structure

The name is the code’s own, after Roverato (2002), whose analysis of Cholesky parameterizations of graph-constrained precision matrices this update exploits. The structure is this:

Since $\boldsymbol{\Theta} = \boldsymbol{\Phi}^\top\boldsymbol{\Phi}$, each off-diagonal precision entry ($i < j$) is

\[ \theta_{ij} = \sum_{l=1}^{i} \Phi_{li}\,\Phi_{lj}, \]

which is linear in the off-diagonal entries of column $j$ of $\boldsymbol{\Phi}$, with coefficients from the already-fixed column $i$. The diagonal entry $\theta_{jj} = \sum_l \Phi_{lj}^2$ is quadratic in those same entries. Taken together, proposing a new value for a single off-diagonal Cholesky entry $\Phi_{ij}$ changes $\theta_{ij}$ and $\theta_{jj}$ while all other precision entries remain fixed, and the diagonal adjustment is deterministic. The resulting $\boldsymbol{\Theta}$ is automatically positive definite because the Cholesky factor stays valid — only one off-diagonal entry changes and diagonal entries remain positive.

The standard approach exploits this by permuting the precision matrix so that the target pair $(i,j)$ occupies the bottom-right corner, then refactorizing — at $O(p^3)$ per element. The bgms implementation avoids this entirely: the same relationship between $\theta_{ij}$, $\theta_{jj}$, and $\Phi_{ij}$ can be read from the cached covariance matrix $\boldsymbol{\Sigma} = \boldsymbol{\Theta}^{-1}$ using cofactor identities, replacing an $O(p^3)$ factorization with an $O(1)$ lookup.

The update algorithm

For a single off-diagonal element $\theta_{ij}$:

Read the constraint relationship from $\boldsymbol{\Sigma}$ and $\boldsymbol{\Theta}$ in $O(1)$.
Propose $\theta_{ij}^* \sim \mathcal{N}(\theta_{ij},\, \sigma_{ij}^2)$ and set $\theta_{jj}^*$ deterministically.
Likelihood ratio via the matrix-determinant lemma — the log-determinant ratio and trace update both use $O(1)$ covariance entries, not a fresh factorization; the tilt adds $\delta$ times the same log-determinant ratio.
Accept or reject via Metropolis–Hastings.
On accept — update $\boldsymbol{\Theta}$, advance the running $\log|\boldsymbol{\Theta}|$ by the determinant-lemma ratio, and apply a Sherman–Morrison–Woodbury rank-2 update to $\boldsymbol{\Sigma}$ ($O(p^2)$), mirroring the result to exact symmetry.

The Cholesky factor itself is not maintained per accept: the sweep reads only $\boldsymbol{\Sigma}$ and the incrementally maintained log-determinant, and one refresh_cholesky() call at the end of the sweep restores the factor, an exact $\boldsymbol{\Sigma}$, and an exact log-determinant for the between-model step and sample recording that follow. A cheap consistency probe on the touched rows of $\boldsymbol{\Sigma}\boldsymbol{\Theta} = \mathbf{I}$ guards the SMW update through near-singular passages: when accuracy is lost, the cache is repaired immediately from $\boldsymbol{\Theta}$.

The proposal standard deviations $\sigma_{ij}$ are tuned during warmup by Robbins–Monro updates targeting an acceptance rate of 0.44 (stage 3b under NUTS uses the same moves through tune_proposal_sd(), consuming the schedule’s Robbins–Monro weight).

Row-block Gibbs

Under update_method = "gibbs", do_one_gibbs_step() sweeps the rows $i = 1, \ldots, p$, drawing each row’s active entries from their exact joint full conditional. Eligibility requires a Normal or Cauchy slab and a Gamma-family scale prior; neither the determinant tilt nor a non-unit Gamma shape gates it (both are handled inside the draw).

Write $\beta = \Theta_{N_i, i}$ for the active entries of row $i$ (neighbours $N_i$, $q = |N_i|$), $A = \Theta_{-i,-i}$, and $C = (A^{-1})_{N_i, N_i}$, computed from the cached covariance by a Schur complement. With slab standard deviation $\sigma$ (on the $-\theta_{ij}/2$ scale), diagonal Gamma rate $\beta_0$, per-edge Cauchy mixture weights $\omega_k$ (fixed at 1 for a Normal slab), and $s_{ii} = S_{ii}$:

$M = (\beta_0 + s_{ii})\, C + \operatorname{diag}\!\big(1 / (4 \sigma^2 \omega_k)\big)$,
$\beta \mid \cdot \sim \mathcal{N}\!\big(-M^{-1} S_{N_i, i},\; M^{-1}\big)$, drawn by triangular solves against the Cholesky factor of $M$,
$\xi \sim \text{Gamma}\!\big(\tfrac{n}{2} + \delta + 1,\; \tfrac{\beta_0 + s_{ii}}{2}\big)$ and $\theta_{ii} = \xi + \beta^\top C\, \beta$.

The determinant tilt enters purely as the $+\delta$ shift in the $\xi$ shape (since $|\boldsymbol{\Theta}| = |A|\,\xi$). For a Gamma shape $\alpha \neq 1$, the extra $(\theta_{ii}/2)^{\alpha-1}$ prior factor does not factorize over $(\beta, \xi)$; the $\alpha = 1$ draw is then used as an independence proposal and corrected by a Metropolis–Hastings accept. For a Cauchy slab, the per-edge scale-mixture weights $\omega_{ij}$ are refreshed from their inverse-Gamma full conditionals each sweep (excluded edges keep prior draws so an edge add can condition on a valid weight). Accepted rows update $\boldsymbol{\Theta}$ and $\boldsymbol{\Sigma}$ through the same rank-2 SMW machinery as the Metropolis path.

Edge selection

When edge selection is active, each edge $(i,j)$ has a binary indicator $\gamma_{ij}$, updated by add/delete moves (Gottardo & Raftery, 2008; van den Bergh et al., 2026) that reuse the Roverato constants and in-place machinery described above. In the proposals below, $c_3$ denotes the scaling constant, precomputed and stored by the model, that relates a proposal on the Cholesky scale to the resulting precision entry.

Delete move ($\gamma_{ij} = 1 \to 0$)

Propose to set $\theta_{ij} = 0$ and $\theta_{jj}$ to its constrained value. The log-acceptance ratio sums:

The $O(1)$ likelihood ratio via the matrix-determinant lemma, plus $\delta$ times the same log-determinant ratio (the tilt)
The prior odds $\log\{(1 - \pi_{ij}) / \pi_{ij}\}$
$-\log J$ under the default hierarchical specification (absent under the joint specification)
The reverse proposal density $+\log\{\phi(\theta_{ij}^{\text{curr}} / c_3;\, 0,\, \sigma_{ij}) / c_3\}$
Minus the slab density of the removed entry, evaluated on the partial-association scale with its Jacobian: $-\log p_{\text{slab}}(-\theta_{ij}^{\text{curr}}/2) + \log 2$
The scale-prior ratio on the slaved diagonal, $\log p_{\text{diag}}(\theta_{jj}^*/2) - \log p_{\text{diag}}(\theta_{jj}/2)$

Add move ($\gamma_{ij} = 0 \to 1$)

Sample $\epsilon \sim \mathcal{N}(0,\, \sigma_{ij})$ and propose $\theta_{ij}^* = c_3 \epsilon$, with $\theta_{jj}^*$ set from the constraint curve. The log-acceptance ratio is the mirror image: likelihood ratio plus tilt, prior odds $\log\{\pi_{ij} / (1 - \pi_{ij})\}$, $+\log J$ under the hierarchical specification, plus the slab density of the new entry (with the $-\log 2$ Jacobian), minus the proposal density, plus the same diagonal-prior ratio.

Rao-Blackwellized inclusion draws

Each edge indicator update — the Metropolis edge move and the conjugate birth/death alike — already forms the log acceptance probability of its proposal. Before the accept decision, the sampler also records the one-step draw

\[ J = \gamma + (1 - 2\gamma)\,\alpha, \]

with $\gamma$ the pre-move indicator state and $\alpha = \exp(\min(0, \log \text{accept}))$: the conditional probability that the edge ends the sweep included, given the rest of the state ($J = \alpha$ for a proposed birth, $1 - \alpha$ for a proposed death). $J$ must be formed here because it uses the pre-move state, which differs from the stored post-sweep indicator on accepted moves. Averaging $J$ over iterations estimates the posterior inclusion probability with lower variance than averaging the binary indicators.

The draws fill a rb_inclusion array in ChainResult through a get_vectorized_rb_inclusion() hook on BaseModel, aligned to the indicator vector order; the OMRF, mixed MRF, and bgmCompare samplers implement the same hook for their indicator updates. Alongside the per-iteration draws, ChainResult keeps per-edge accumulators (rb_counts) that sum the acceptance probabilities of birth and death proposals separately on the acceptance-probability scale — never forming $1 - \alpha$ per draw — so posterior inclusion odds stay finite even when the averaged draws round to 0 or 1. On the R side the draws surface as fit$raw_samples$rb_inclusion and back the default "rb" estimator of extract_posterior_inclusion_probabilities(); the accumulators back extract_inclusion_bf(). Accept/reject behaviour and the RNG stream are unchanged by the recording.

Iteration order

The order of edge updates is randomized at each iteration via prepare_iteration(), which shuffles the flat edge index with arma_randperm(). After all edge indicator updates complete, one refresh_cholesky() restores exact state, the constraint_dirty_ flag is set, and the constraint structure is rebuilt lazily at the next gradient evaluation. This ensures NUTS operates on an up-to-date graph.

Conjugate edge birth/death (Gibbs path)

The Gibbs sampler replaces the random-walk edge proposal with a full-conditional birth/death move (update_edge_indicator_conjugate()), enabled by set_conjugate_edge_proposal(true). The move works in the cofactor coordinate of the Roverato parameterization, where toggling the edge preserves $|\boldsymbol{\Theta}|$ — the determinant tilt and the likelihood determinant term cancel exactly — so the full conditional of the edge entry is Gaussian and can be proposed from directly. The acceptance ratio then reduces to the inclusion odds times $p_{\text{slab}}(0) / q(0)$, evaluated at the spike: it does not depend on the proposed value, and no tuning is needed. A Cauchy slab enters through its current mixture weight $\omega_{ij}$; a Gamma shape $\alpha \neq 1$ adds the independence-proposal correction $(\theta_{jj}/\theta_{jj}^{(0)})^{\alpha - 1}$; and under the hierarchical specification the same $\pm \log J$ Z-ratio term applies as in the random-walk moves.

Warm starts and initialization

With data present, the precision matrix is initialized from a regularized MLE: $\boldsymbol{\Theta}_0 = n\,(\mathbf{S} + r_0 \mathbf{I})^{-1}$ with ridge $r_0 = \operatorname{tr}(\mathbf{S}) / (p\, n)$, shrinking toward the identity at a scale-appropriate rate. (The ridge is written $r_0$ here to keep it apart from the determinant tilt $\delta$, a different quantity.) If the initial graph excludes edges, the excluded entries are zeroed and the diagonal is inflated to diagonal dominance so the start stays positive definite. If the inversion fails, the identity is kept.

Prior-only mode

With $n = 0$ (used by sample_ggm_prior()), the likelihood terms vanish and the same machinery samples the prior alone. One safeguard becomes necessary: with data, the $\tfrac{n}{2}\log|\boldsymbol{\Theta}|$ term vetoes proposals that leave the positive-definite cone, but without data that anchor is gone, and floating-point drift in the incrementally maintained covariance can hand a non-positive-definite proposal a finite acceptance probability. The Metropolis moves therefore run an explicit positive-definiteness check on the proposal when $n = 0$ and reject failures outright.

Log-likelihood

The GGM log-likelihood for $n$ observations is:

\[ \ell(\boldsymbol{\Theta}) = \frac{n}{2} \log |\boldsymbol{\Theta}| - \frac{1}{2} \operatorname{tr}(\mathbf{S}\,\boldsymbol{\Theta}) \]

The log-determinant is computed from the Cholesky factor as $\log|\boldsymbol{\Theta}| = 2\sum_q \log \Phi_{qq}$ (and advanced incrementally between refreshes). The trace term uses the precomputed sufficient statistic $\mathbf{S}$.

Missing data imputation

When na_action = "impute", the GGM imputes missing values using the conditional normal distribution implied by the current precision matrix. After imputation, the sufficient statistic $\mathbf{S}$ is recomputed for the next likelihood evaluation.

References

Gottardo, R., & Raftery, A. E. (2008). Markov chain Monte Carlo with mixtures of mutually singular distributions. Journal of Computational and Graphical Statistics, 17(4), 949–975. https://doi.org/10.1198/106186008X386102

Roverato, A. (2002). Hyper inverse Wishart distribution for non-decomposable graphs and its application to Bayesian inference for Gaussian graphical models. Scandinavian Journal of Statistics, 29(3), 391–411. https://doi.org/10.1111/1467-9469.00297

van den Bergh, D., Clyde, M. A., Raftery, A. E., & Marsman, M. (2026). Reversible jump MCMC with no regrets: Bayesian variable selection using mixtures of mutually singular distributions. Manuscript in Preparation.

GGM