Mixed MRF
Implementation details of the mixed Markov random field in bgms. Source: src/models/mixed/mixed_mrf_model.h and four .cpp files (mixed_mrf_model.cpp, mixed_mrf_gradient.cpp, mixed_mrf_likelihoods.cpp, mixed_mrf_metropolis.cpp).
Block structure
The mixed MRF handles \(p\) discrete and \(q\) continuous variables. The parameter space is organized into five blocks — two main-effect blocks and three interaction blocks:
| Block | Parameters | Size |
|---|---|---|
| Thresholds | Discrete category thresholds | \(\text{thresholds}\) |
| Means | Continuous variable means | \(q\) |
| Discrete-discrete | Pairwise discrete interactions (symmetric) | \(p(p-1)/2\) |
| Cross-type | Discrete-continuous interactions | \(p \times q\) |
| Continuous-continuous | Cholesky factor of the precision matrix | \(q(q+1)/2\) |
All five blocks are updated jointly by NUTS. The continuous-block precision matrix \(\boldsymbol{\Theta}\) is parameterized through its upper Cholesky factor \(\boldsymbol{\Phi}\) with \(\boldsymbol{\Theta} = \boldsymbol{\Phi}^\top \boldsymbol{\Phi}\) (diagonal entries on the log scale, the same convention as GGM Internals), so that positive definiteness is guaranteed by construction. The gradient includes the Jacobian of the Cholesky-to-precision mapping.
Parameter vector layout
The model distinguishes three vector types:
Active vector (parameter_dimension())
Contains the parameters that NUTS updates, excluding edges set to zero by edge selection. The Cholesky block is always fully included:
[ thresholds | discrete interactions (active) | means | cross-type (active) | Cholesky (all) ]The mass matrix and gradient are dimensioned to this vector.
Full vector (full_parameter_dimension())
Contains all parameters including inactive edges:
[ thresholds | discrete interactions (all) | means | cross-type (all) | Cholesky (all) ]Used for storage and initialization.
Storage vector (storage_dimension())
Identical layout to the full vector. This is what gets saved to ChainResult::samples at each iteration.
Sampling strategy
Each MCMC iteration updates the full active parameter vector (all five blocks) in a single NUTS step. The logp_and_gradient() method computes the log-pseudoposterior and gradient for the entire vector, including the Cholesky block with its log-determinant Jacobian.
When edge selection is active (or the initial graph is sparse), the model reports has_constraints() = true and the NUTS sampler switches to constrained (RATTLE) integration, using logp_and_gradient_full() and projection callbacks to enforce zero entries for excluded edges (see Constrained Leapfrog (RATTLE)).
Edge indicator updates (when active) operate on all three interaction blocks.
The model also provides a component-wise adaptive Metropolis fallback (do_one_metropolis_step()), used when the sampler is set to "adaptive-metropolis". In Metropolis mode, each parameter is updated individually with Robbins–Monro-tuned proposal SDs targeting an acceptance rate of 0.44. Off-diagonal precision updates use a rank-2 Cholesky update; diagonal updates use rank-1.
Gradient computation
The gradient (mixed_mrf_gradient.cpp) covers all five parameter blocks in the NUTS vector. It involves:
- Discrete conditionals — Category probabilities via Variable Helpers, contributing gradients for thresholds and discrete pairwise interactions
- Cross-type terms — Contributions from discrete-continuous interactions to both discrete and continuous conditional distributions
- Continuous means — Gradient of the continuous-variable contribution to the pseudolikelihood
- Cholesky block — The precision gradient \(\bar{\boldsymbol{\Theta}}\) is computed from the GGM conditional contribution and the conditional-mean coupling, then mapped to the Cholesky gradient \(\bar{\boldsymbol{\Phi}} = \boldsymbol{\Phi}(\bar{\boldsymbol{\Theta}} + \bar{\boldsymbol{\Theta}}^\top)\). Diagonal entries use the log-scale Jacobian: \(\partial\ell / \partial\psi_j = \bar{\Phi}_{jj} \Phi_{jj} + (q - j + 1)\), where \(\psi_j = \log \Phi_{jj}\) and \(j\) runs over the \(q\) continuous variables (\(j = 1, \ldots, q\)). Throughout, a bar denotes the reverse-mode adjoint \(\bar{A} = \partial \ell / \partial A\) (as in GGM Internals).
The log-prior gradient adds the interaction-prior penalties (default Cauchy) on off-diagonal precision entries and the scale-prior penalties (default Gamma(1, 1)) on diagonal entries; see Parameter Priors.
Discrete pseudolikelihood
Write \(x_s\) for discrete variable \(s\), \(x_{-s}\) for the remaining discrete variables, and \(\mathbf{D}\) for the \(n \times q\) column-centered continuous data matrix. The discrete full conditionals integrate out the continuous block: the model uses \(\log p(x_s \mid x_{-s})\) rather than conditioning on the observed continuous values. The continuous block uses the exact GGM likelihood \(\log p(y \mid x) \propto \tfrac{n}{2}\log|\boldsymbol{\Theta}| - \tfrac{1}{2}\operatorname{tr}(\boldsymbol{\Theta} \mathbf{D}^\top \mathbf{D})\).
For a cross-model overview and the derivation of the marginal form, see Pseudolikelihood.
After integrating out the continuous block, the discrete conditionals use the marginal interaction matrix
\[ \mathbf{M} = \mathbf{A}_{xx} + 2\,\mathbf{A}_{xy}\,\boldsymbol{\Sigma}_{yy}\,\mathbf{A}_{xy}^\top, \]
where \(\mathbf{A}_{xx}\) is the discrete-discrete interaction block, \(\mathbf{A}_{xy}\) the cross-type (discrete-continuous) interaction block, and \(\boldsymbol{\Sigma}_{yy} = \boldsymbol{\Theta}^{-1}\) the conditional covariance of the continuous variables. (In the notation of the mixed-model guide page, these are the \(\boldsymbol{\Omega}^{xx}\) and \(\boldsymbol{\Omega}^{xy}\) blocks of the interaction matrix.) \(\mathbf{M}\) must be recomputed whenever the precision or cross-type parameters change. The \(\partial \mathbf{M}/\partial\boldsymbol{\Sigma}_{yy}\) chain rule contributes to the precision gradient via \(-2\boldsymbol{\Sigma}_{yy} \mathbf{A}_{xy}^\top\bar{\mathbf{M}} \mathbf{A}_{xy}\boldsymbol{\Sigma}_{yy}\).
Missing data
The mixed MRF tracks missing indices separately for discrete and continuous variables (missing_index_discrete, missing_index_continuous). During imputation:
- Missing discrete values are sampled from their full conditional categorical distributions
- Missing continuous values are sampled from their conditional normal distributions given the current precision matrix
Output assembly
The R output builder (build_output_mixed_mrf() in build_output.R) maps the flat storage vector back to the original variable ordering using fill_mixed_symmetric(). This function places discrete-discrete, continuous-continuous, and cross-type interactions into the correct cells of the \((p+q) \times (p+q)\) output matrix, respecting the original column order (which may interleave discrete and continuous variables).