A Scalable MCMC Algorithm for Bayesian Inference on Binary Model Spaces

1. Problem Statement

The paper addresses the challenge of Bayesian model inference on binary model spaces. In many statistical applications (e.g., variable selection, graphical models, Ising models, decision trees), the model space consists of binary vectors $m \in \{0, 1\}^k$, where each bit represents the inclusion or exclusion of a specific component (e.g., a variable or an edge in a graph).

• The Core Difficulty: The size of the model space grows exponentially ($2^k$). Traditional Markov Chain Monte Carlo (MCMC) methods, such as Reversible Jump MCMC (RJ-MCMC) or standard Birth-Death (BD) processes, suffer from slow exploration (poor mixing) because they typically propose local moves (changing only one bit per iteration).
• Computational Bottleneck: For large-scale problems (e.g., $k \approx 10^6$ potential edges in a 1000-node graph), BD algorithms require millions of iterations to traverse the space, leading to prohibitive computation times and memory issues.
• Goal: Develop an MCMC algorithm that is rejection-free, capable of global updates (flipping multiple bits simultaneously), and maintains high accuracy in targeting the posterior distribution while being orders of magnitude faster than state-of-the-art methods.

2. Methodology: The Multiple Jump MCMC (MJ-MCMC)

The authors propose the Multiple Jump MCMC (MJ-MCMC) algorithm. The core innovation is transforming a continuous-time Birth-Death (BD) process into a discrete-time Markov chain that allows for simultaneous updates of multiple model components.

2.1. Theoretical Foundation

• Starting Point: The algorithm begins with a standard BD process (Lemma 2.1) defined by transition rates $Q(m, m')$. In a BD process, the chain moves from state $m$ to a neighbor $m_i$ (differing by one bit) with rate $q_i(m) = \min(1, \frac{p(m_i|y)}{p(m|y)})$. This process is rejection-free but slow.
• The Transformation: The authors construct a Discrete-Time Markov Chain (DTMC) where, in a single iteration $s$, every bit $i \in \{1, \dots, k\}$ is flipped independently with probability $q_i(m)\epsilon_s$.
  • Here, $\epsilon_s \in (0, 1)$ is a user-defined sequence.
  • The transition probability from $m$ to $m'$ is the product of probabilities for bits that flip and bits that stay the same:
    $$P_{\epsilon_s}(m, m') = \prod_{i \in H_{m,m'}} q_i(m)\epsilon_s \prod_{i \notin H_{m,m'}} (1 - q_i(m)\epsilon_s)$$
    where $H_{m,m'}$ is the set of indices where $m$ and $m'$ differ.
• Key Feature: Unlike BD processes which move to a single neighbor, MJ-MCMC can jump to any state in the model space in a single iteration (provided the posterior probability is non-zero).

2.2. Variants and Convergence

The paper analyzes two variants based on the sequence $\{\epsilon_s\}$:

1. Homogeneous Case ($\epsilon_s = \epsilon$):
   • The chain is a time-homogeneous DTMC.
   • Theorem 3.1: As $\epsilon \to 0$, the stationary distribution of the MJ-MCMC converges to the invariant distribution of the original BD process. The error is $O(\epsilon)$.
   • Trade-off: Small $\epsilon$ ensures accuracy but slows mixing; large $\epsilon$ speeds up mixing but may introduce bias (though experiments show this bias is negligible for sparse graphs).
2. Inhomogeneous Case (Decaying $\epsilon_s$):
   • The sequence $\epsilon_s$ decays to 0 but its sum diverges ($\sum \epsilon_s = \infty$).
   • Theorem 3.2: Under these conditions, the limiting distribution of the MJ-MCMC is exactly the same as the BD process's invariant distribution (the true posterior), regardless of the initial speed of mixing.

2.3. Extension to Joint Posterior

While the primary algorithm targets the marginal posterior $p(m|y)$, Lemma 3.4 provides a strategy to sample from the joint posterior $p(m, \theta|y)$. After generating a sample $m^{(s)}$ via MJ-MCMC, one simply samples the associated parameters $\theta^{(s)}$ from the conditional posterior $p(\theta|m^{(s)}, y)$.

3. Key Contributions

1. Algorithmic Innovation: Introduction of MJ-MCMC, a rejection-free algorithm that allows global updates (flipping multiple bits simultaneously) while maintaining the theoretical guarantees of a BD process.
2. Theoretical Guarantees: Rigorous proofs (Theorems 3.1 and 3.2) demonstrating that the algorithm targets the correct posterior distribution under specific conditions on the step-size sequence $\epsilon_s$.
3. Scalability: The algorithm is designed to traverse the entire model space in a single iteration, making it exceptionally scalable for high-dimensional problems ($k \approx 10^6$).
4. General Applicability: The framework is not limited to one model type; it is applicable to any Bayesian inference problem on binary model spaces, including Gaussian Graphical Models (GGMs), Ising models, and variable selection.

4. Results

The authors evaluated MJ-MCMC primarily on Undirected Gaussian Graphical Models (GGMs), comparing it against the state-of-the-art BD-MPL algorithm (Mohammadi et al., 2024).

4.1. Simulation Study

• Setup: $p=1000$ variables, $n=400$ observations, sparse graphs (0.2% edge density).
• Performance:
  • Speed: MJ-MCMC solved 1000-variable problems in under 30 seconds. This is 100 to 200 times faster than the BD-MPL algorithm.
  • Accuracy: Despite the massive speedup, the edge inclusion probabilities estimated by MJ-MCMC were nearly identical to BD-MPL. The Area Under the Precision-Recall Curve (AUC-PR) scores were comparable across various $\epsilon$ values (including high values like 0.95).
  • Robustness: Even with large $\epsilon$ (which theoretically introduces bias), the algorithm performed well on sparse graphs, likely because the posterior behaves approximately factorizable in these settings (Lemma 3.3).

4.2. Real-World Data Application

• Dataset: GSE15907 microarray dataset (623 genes, 653 mouse immune cells).
• Results:
  • MJ-MCMC achieved the same quality of inference as BD-MPL (correlation of edge inclusion probabilities > 0.988).
  • Efficiency: MJ-MCMC was up to 10 times faster in wall-clock time and required up to 600 times fewer MCMC iterations to reach convergence.

5. Significance and Impact

• Solving the "Curse of Dimensionality" in Model Selection: The paper demonstrates that high-dimensional Bayesian model selection (e.g., learning graphs with thousands of nodes) is computationally feasible where it was previously intractable due to slow mixing of traditional MCMC.
• Paradigm Shift: It challenges the necessity of local moves in MCMC. By proving that global, rejection-free jumps can target the correct posterior, it opens new avenues for designing efficient samplers.
• Practical Utility: The ability to solve 1000-node graphical models in seconds rather than hours makes Bayesian structure learning accessible for real-time or large-scale genomic and network analysis.
• Future Directions: The authors suggest extending the framework to mixed graphical models, decision trees, and continuous model spaces, as well as optimizing computational costs further by updating subsets of rates.

In summary, the paper presents a breakthrough in Bayesian computation, offering a simple yet powerful algorithm that drastically reduces the computational cost of exploring complex binary model spaces without sacrificing statistical accuracy.