Parallel Chains

bgms runs multiple MCMC chains in parallel using Intel TBB (Threading Building Blocks), accessed through the RcppParallel package. The implementation is in src/mcmc/execution/chain_runner.h and chain_runner.cpp. Within each chain, the iteration loop drives the samplers described in Sampler Hierarchy under the adaptation schedule described in Warmup Schedule.

Design

The entry point is run_mcmc_sampler(), which takes a prototype model, a prototype edge prior, and a sampler configuration. It follows a clone-and-dispatch pattern:

1. Clone — Create independent copies of the model and edge prior for each chain via model.clone() and edge_prior.clone(). Each clone has its own internal state and RNG.
2. Seed — Each clone is seeded with config.seed + chain_id to ensure reproducibility while avoiding correlated chains.
3. Dispatch — All chains are dispatched to threads. Each thread runs run_mcmc_chain() with its own model, prior, and result container.
4. Collect — Results are gathered and returned as a vector of ChainResult objects.

After a chain’s main loop finishes, three post-loop steps run before the chain returns. If the Z-ratio trust gauge is enabled (options(bgms.zratio_gauge_sweeps = ...) under the hierarchical graph prior), the chain runs the configured number of assessment sweeps on the frozen kernel: it reruns its own edge-selection pass (the same code path that ran during sampling) and, for each non-trivial edge move (one with a non-empty mediating block), also evaluates the exact block-local reference for the same ratio, so the fast decision and the exact one can be compared. No stored samples are touched. Next the sampler’s adaptation-averaged step size and diagonal metric are retained on the ChainResult (final_step_size, final_inv_mass; NaN and empty for non-gradient samplers) so a refit can warm-start them (see Sensitivity and Refits). Then model.collect_chain_diagnostics(chain_result) copies run-level diagnostic state (Z-ratio counters, frozen constants, gauge results) into the ChainResult.

Model and prior cloning

Both BaseModel::clone() and BaseEdgePrior::clone() perform deep copies. For the model, this includes the full parameter state, Cholesky factors (GGM), residual matrices (OMRF), and all internal buffers. For the SBM edge prior, the cluster allocations and block probability matrix are cloned.

The clone methods return unique_ptr, so each chain owns its objects exclusively with no shared mutable state between threads.

TBB thread control

Before dispatching, bgms sets the maximum parallelism via TBB’s global control:

tbb::global_control control(
    tbb::global_control::max_allowed_parallelism,
    no_threads
);

This limits the TBB thread pool to no_threads (the cores argument in R), respecting the user’s request even when the system has more cores available.

RcppParallel dispatch

The MCMCChainRunner struct inherits from RcppParallel::Worker and implements operator()(begin, end). Each unit of work is one chain:

void MCMCChainRunner::operator()(std::size_t begin, std::size_t end) {
    for (std::size_t i = begin; i < end; ++i) {
        ChainResult& chain_result = results_[i];
        BaseModel& model = *models_[i];
        BaseEdgePrior& edge_prior = *edge_priors_[i];
        model.set_seed(config_.seed + static_cast<int>(i));

        try {
            run_mcmc_chain(chain_result, model, edge_prior, config_,
                           static_cast<int>(i), pm_, warm_eps, warm_metric);
        } catch (std::exception& e) {
            chain_result.error = true;
            chain_result.error_msg = e.what();
        } catch (...) {
            chain_result.error = true;
            chain_result.error_msg = "Unknown error";
        }
    }
}

The try/catch is what turns a thrown chain into the error and error_msg fields the R output builder reads, rather than an exception crossing a worker thread. warm_eps and warm_metric are the per-chain warm-start step size and diagonal metric, NaN and empty unless a refit supplied them.

Dispatch uses RcppParallel::parallelFor(0, no_chains, runner).

Sequential fallback

When no_threads == 1, bgms skips TBB and runs chains sequentially in a simple loop. This avoids TBB overhead for single-chain runs and simplifies debugging. The branch condition checks threads, not chain count, so multiple chains with no_threads = 1 still run sequentially.

Per-chain RNG seeding

Each chain’s RNG is a SafeRNG wrapping dqrng::xoshiro256plusplus. The seed for chain \(k\) is config.seed + k, ensuring:

• Reproducibility: the same seed and chain count always produce the same results
• Independence: different chains explore different random number sequences

Progress and interrupts

The ProgressManager class handles progress display and user interrupt detection across threads.

• Progress bar — Every chain counts its iterations into atomic per-chain counters. The manager aggregates them and displays a single bar (or per-chain bars, depending on progress_type).
• User interrupts — Each chain checks pm.shouldExit() between iterations. If the user presses Ctrl-C, the flag is set and all chains exit their sampling loops. The partial results are still returned (with a warning) rather than discarded.
• Thread contract — The manager must be constructed on the R main thread, and all R API interaction — the interrupt check, console output, and the R-level progress callback — happens only on that thread. The main thread always participates in the chain pool, so its update() calls drive the display; worker threads only bump their atomic counters and never touch the R interpreter.

ChainResult storage

Each chain writes its output into a preallocated ChainResult:

Field	Shape	Content
samples	param_dim × n_iter	Parameter samples (post-warmup)
indicator_samples	n_edges × n_iter	Edge indicators (if edge selection)
rb_inclusion_samples	n_edges × n_iter	Rao-Blackwellized inclusion draws (if edge selection)
rb_counts	n_edges × 4	Per-edge odds accumulators [n01, n10, n0_visits, n1_visits] on the alpha scale (if edge selection)
allocation_samples	n_vars × n_iter	SBM cluster allocations (if SBM prior)
inclusion_parameter_samples	n_iter	Sampled inclusion probability (if Beta-Bernoulli edge prior on a continuous or mixed model)
treedepth_samples	n_iter	NUTS tree depth (if NUTS)
divergent_samples	n_iter	Divergence flags (if NUTS)
energy_samples	n_iter	Hamiltonian energy (if NUTS)
accept_prob_samples	n_iter	NUTS mean per-trajectory acceptance probability (if NUTS)
am_accept_prob_samples	n_iter	Mean per-iteration acceptance across updated components (if adaptive Metropolis)
zratio_addc, zratio_counters, zratio_gauge_*	—	Z-ratio engine state and trust-gauge results (if hierarchical graph prior)

Storage is preallocated with “not sampled” sentinels — -1 for the integer buffers, NaN for the floating ones — so an interrupted run returns recognizable placeholders rather than uninitialized values for iterations that never ran.

The ChainResult also carries error status and user interrupt flags per chain. These are checked by the R output builder to decide whether to warn or error.