Generalizing Metropolis-Hastings for DAG Sampling

Daniele Nikzad; Mikko Koivisto

cs.LGOct 21, 2026

Generalizing Metropolis-Hastings for DAG Sampling

Authors: Daniele Nikzad, Mikko Koivisto

Organizations: Department of Computer Science, University of Helsinki, Helsinki, FI-00014, Finland

Abstract

A fundamental challenge in Bayesian network structure learning lies in efficiently approximating the posterior distribution of directed acyclic graphs (DAGs). For large, high-dimensional datasets, the task remains beyond the capabilities of current methods. This work advances Markov chain Monte Carlo algorithms by introducing novel moves that simultaneously resample the parents for multiple nodes of the DAG. At the core of these moves is duet sampling, a new generalization of the Metropolis-Hastings algorithm that selects two proposal distributions from a pool at each iteration. One is used for generating a candidate state of the Markov chain and another for evaluating the return to the current state; both are ingredients of an appropriately crafted acceptance criterion. Pairing different proposal distributions provides additional freedom and enables constructing Markov chains that mix more efficiently than their basic Metropolis-Hastings counterparts, at no additional computational cost. Experiments show that our new moves, both individually and in combination, outperform the previous state of the art, which relies on a more restrictive pool of moves. This advantage persists even when letting the previous methods leverage duet sampling to enhance their performance.

Figures & tables

Panel of Figure 1 — Figure 1: Illustration of an EC move transforming the current DAG $G$ into the proposal $G^{\prime}$ . First, a node pair is sampled. In this example the pair is $(2,4)$ . Since there is a directed edge from $2$ to $4$ , the corresponding edge constraint is $c=\mathrm{rem}$ . The nodes $2$ and $4$ are then orphaned, yielding the intermediate DAG $G_{0}$ . Next, a new parent set $P^{\prime}_{2}$ is sampled from those that ensure acyclicity and satisfy the edge constraint $c$ (in this case, node 4 must not be a parent of 2), producing $G^{\prime}_{2}$ . Finally, a new parent set $P^{\prime}_{\text{4}}$ is sampled from those that ensure acyclicity and satisfy the edge constraint $c$ (node 2 must not be a parent of 4), generating the proposal $G^{\prime}$ .

Panel of Figure 2 — Figure 2: Illustration of an AC move transforming the current DAG $G$ into the proposal $G^{\prime}$ . First, a node pair is sampled. In this example, the pair is $(\text{3},\text{2})$ . Since 2 is an ancestor of 3 in $G_{0}$ , the corresponding ancestor constraint is $c=\text{rev}$ , thus the move forces 3 to be an ancestor of 2 in $G^{\prime}$ . The nodes 3 and 2 are then orphaned, yielding the intermediate DAG $G_{0}$ . Next, a new parent set $P^{\prime}_{3}$ is sampled from those that maintain acyclicity and satisfy the ancestor constraint $c$ (in this case, node 2 and its descendants in $G_{0}$ must not be parents of $3$ ). This produces the DAG $G^{\prime}_{3}$ . Finally, a new parent set $P^{\prime}_{2}$ is sampled from those that guarantee acyclicity and satisfy the ancestor constraint $c$ (node 3 must be an ancestor of 2), generating the DAG $G^{\prime}$ .

Panel of Figure 3 — Figure 3: Illustration of a USR move transforming the current DAG $G$ into the proposal $G^{\prime}$ . First, a permutation $\sigma$ is sampled. Here $\sigma=(1,3,2)$ . The nodes in $\sigma$ are orphaned, yielding the intermediate DAG $G_{0}$ . Next, the new parent sets of the nodes in $\sigma$ are sequentially sampled from those that preserve acyclicity, according to the order given by $\sigma$ , and incorporated into $G_{0}$ . In detail, we obtain the DAG $G^{\prime}_{\sigma,1}$ after sampling $P^{\prime}_{1}$ , the DAG $G^{\prime}_{\sigma,2}$ after sampling $P^{\prime}_{3}$ , and finally $G^{\prime}$ after sampling $P^{\prime}_{2}$ . In this particular case, the only candidate parent set of node 2 is the empty parent set as it is the only one that preserves acyclicity.

Network	#nodes	#edges	MID
H ailfinder	56	66	4
A ndes	223	338	6
P igs	441	592	2

Table 1: Summary of the benchmark networks used in the experiments.

	#steps
Algorithm	Hailfinder	Andes	Pigs
USR	$2.6\text{\times}{10}^{8}$	$1.2\text{\times}{10}^{7}$	$3.1\text{\times}{10}^{7}$
EC & AC	$7.7\text{\times}{10}^{8}$	$7.3\text{\times}{10}^{7}$	$1.1\text{\times}{10}^{8}$
FBM	$1.7\text{\times}{10}^{12}$	$1.7\text{\times}{10}^{12}$	$1.5\text{\times}{10}^{12}$
REV	$1.9\text{\times}{10}^{9}$	$6.1\text{\times}{10}^{7}$	$6.2\text{\times}{10}^{8}$
MBR	$1.6\text{\times}{10}^{9}$	$7.5\text{\times}{10}^{7}$	$2.1\text{\times}{10}^{8}$

Table 2: Number of Markov chain steps taken by the single-move algorithms in 15 000 s.

Panel of Figure 4 — Figure 4: Traceplots of MCMC runs for data sets sampled from large benchmark networks. Columns correspond to data sets. Rows correspond to (groups of) algorithms. For each algorithm, ten independent runs are displayed. For each run, the log-unnormalized posterior probability (log-score) is shown as a function of running time for 100 evenly spaced DAG samples.

	Hailfinder			Andes			Pigs
	$\ell_{*}=-497\,506.9$			$\ell_{*}=-934\,406.0$			$\ell_{*}=-3\,321\,511.0$
Algorithm	$\Delta^{\text{min}}$	$\Delta^{\text{max}}$	$\mathrm{WEPD}$	$\Delta^{\text{min}}$	$\Delta^{\text{max}}$	$\mathrm{WEPD}$	$\Delta^{\text{min}}$	$\Delta^{\text{max}}$	$\mathrm{WEPD}$
USR-rev	0	0	0.0074	0	0	$0.027$	0	0	$0.025$
USR-rnd	0	0	$0.095$	0	0	$0.044$	0	0	0.018
EC	0	$5\text{\times}{10}^{2}$	$1.00$	0	0	$0.29$	0	0	$0.034$
AC	0	$1\text{\times}{10}^{3}$	$1.00$	0	0	$0.080$	0	0	$0.032$

Table 3: Comparison of the highest scores found and the differences of edge probability estimates across 10 independent MCMC runs. Smaller values are better.

Appendix figures & tables4 assets

Supplementary material from the paper’s appendix.

Appendix

Algorithm 1 iMCMC (single step) 1: input target density

\pi(x)

2: input proposal distribution

q(v|x)

3: input involution

f(x,v)

4: input current state

x

5: sample

v\sim q(v|x)

6: propose

(x^{\prime},v^{\prime}):=f(x,v)

7: compute the acceptance probability

\displaystyle\alpha:=\min\left\{1,\,\frac{\pi(x^{\prime})\,q(v^{\prime}|x^{\prime})}{\pi(x)\,q(v|x)}\left|\frac{\partial f(x,v)}{\partial[x,v]}\right|\right\}

8: output

\displaystyle\begin{cases}x^{\prime}&\text{with probability }\alpha\,,\\ x&\text{with probability }1-\alpha\end{cases}

Algorithm 2 Duet sampling (single step) 1: input target density

\pi(x)

2: input proposal distributions

\{q_{k}\}_{k\in L}

3: input pairing

\{\beta_{kl}\}

\{q_{k}\}

4: input current state

x

5: sample

(x^{\prime},k,l)\sim\beta_{kl}(x)q_{k}(x^{\prime}|x)

\triangleright

New state

x^{\prime}

got already proposed 7: compute the acceptance probability

\displaystyle\alpha:=\min\left\{1,\,\frac{\pi(x^{\prime})\,\beta_{lk}(x^{\prime})\,q_{l}(x|x^{\prime})}{\pi(x)\,\beta_{kl}(x)\,q_{k}(x^{\prime}|x)}\right\}

8: output

\displaystyle\begin{cases}x^{\prime}&\text{with probability }\alpha\,,\\ x&\text{with probability }1-\alpha\end{cases}

Appendix

Table 8

	Ancestor relationship between $i$ and $j$ in $G^{\prime}$
Descendants	$i\rightsquigarrow j$	$j\rightsquigarrow i$	$i\not\rightsquigarrow j$ , $j\not\rightsquigarrow i$
$D_{i}(G^{\prime}_{j})$	$D_{i}(G_{0})\cup D_{j}(G_{0})$	$D_{i}(G_{0})$	$D_{i}(G_{0})$
$D_{j}(G^{\prime}_{i})$	$D_{j}(G_{0})$	$D_{i}(G_{0})\cup D_{j}(G_{0})$	$D_{j}(G_{0})$

Appendix

Table 4: Descendants of nodes

i

and

j

in the intermediate DAGs

G^{\prime}_{i}

and

G^{\prime}_{j}

along the forward path from

G

G^{\prime}

Network	#nodes	#edges	MID
A sia	8	8	2
S achs	11	17	3

Appendix

Table 5: Summary of the benchmark networks used in the experiments.

Panel of Figure 5 — Figure 5: The MAD of the best and worst runs of each algorithm on low-dimensional datasets sampled from benchmark networks. Columns correspond to datasets, while rows correspond to different algorithms.

Bayesian Network Structure Learning

See all 13 papers

Oct 6, 2026

Directional Evidence Guided Search-Space Reduction for Exact DAG Learning

Learning a directed acyclic graph (DAG) from observational data is a challenging combinatorial problem due to the exponential growth in the number of candidate parent-set configurations. Existing exact score-based methods often require computationally intensive combinatorial search, whereas constraint-based methods can become unreliable or computationally demanding as graph size and conditioning-set complexity increase. We develop a non-parametric hybrid framework, referred to as DECO (Directional Evidence-guided Configuration Optimization), that extracts dependency and directional evidence from observation data to construct admissible parent sets prior to exact optimization. It reduces the optimization search space by eliminating empirically unsupported parent configurations while preserving flexibility for all plausible edge orientations. Theoretical analysis establishes an exponential reduction in the admissible parent-set configuration space and quantifies how bounded edge-level omission affects the probability of retaining the true parent structure. Experiments on benchmark Bayesian networks and synthetic discrete and continuous DAGs demonstrate substantial search-space reduction while achieving competitive structure-recovery performance, with favorable structural Hamming distance across many evaluated settings. These results show that directional evidence can provide an effective preprocessing mechanism for reducing the computational burden of exact DAG learning without requiring a fixed parametric structural~model.

Upala Junaida Islam, Abdelmonem Elrefaey, Rong Pan

School of Computing and Augmented Intelligence, Arizona State University, Tempe, AZ 85281, USA.

Score-Based Causal Discovery Bayesian Network Structure Learning Causal Discovery Causal Structure Learning Bayesian Networks

Ask PDF

May 1, 2026

Decentralized Proximal Stochastic Gradient Langevin Dynamics

We propose Decentralized Proximal Stochastic Gradient Langevin Dynamics (DE-PSGLD), a decentralized Markov chain Monte Carlo (MCMC) algorithm for sampling from a log-concave probability distribution constrained to a convex domain. Constraints are enforced through a shared proximal regularization based on the Moreau-Yosida envelope, enabling unconstrained updates while preserving consistency with the target constrained posterior. We establish non-asymptotic convergence guarantees in the 2-Wasserstein distance for both individual agent iterates and their network averages. Our analysis shows that DE-PSGLD converges to a regularized Gibbs distribution and quantifies the bias introduced by the proximal approximation. We evaluate DE-PSGLD for different sampling problems on synthetic and real datasets. As the first decentralized approach for constrained domains, our algorithm exhibits fast posterior concentration and high predictive accuracy.

Mohammad Rafiqul Islam, Lingjiong Zhu

Department of Mathematics, Florida State University, Tallahassee, Florida, United States of America

Markov Chain Monte Carlo Stochastic Gradient Langevin Dynamics Langevin Dynamics Log-Concave Sampling Posterior Sampling Constrained Sampling

Ask PDF

Aug 5, 2026

SVI-DAG: A Structured Variational Inference Approach to Bayesian Causal Discovery

Bayesian causal discovery seeks to determine the posterior distribution of causal theories, which are interpreted as directed acyclic graphs (DAGs) that explain the observed data. The resulting posterior allows systematic reasoning regarding epistemic uncertainty within these theories. Nonetheless, finding such graphs is difficult due to identifiability problems and limited observational data. Furthermore, precisely approximating posterior over graphs is challenging given vast range of potential DAGs. Recent Bayesian approaches have addressed some of these challenges, yet they remain limited as they fail to encode dependencies between edges, and lack principled ways to incorporate domain knowledge as inductive biases during the search process. To overcome these limitations, we propose SVI-DAG, a structured variational inference approach to Bayesian causal discovery using observational data and prior beliefs that uses normalizing flows to model dependencies between edges, supporting expressive and multimodal posterior learning over DAGs. To mitigate mode seeking behaviour in evidence lower bound optimization and promote mode coverage, we use stein variational gradient descent to update the node potentials using a kernel in acyclicity space. We evaluate SVI-DAG against 5 state-of-the-art Bayesian DAG learning methods and demonstrate superior performance in uncertainty quantification while remaining competitive in terms of structural accuracy.

Shrenik Zinage

Massachusetts Institute of Technology · This work was conducted while the author was a PhD student at Purdue University, independently of the author’s dissertation research.

Causal Discovery Variational Inference

Ask PDF

Abstract

Bayesian Network Structure Learning

Explore similar work

Directional Evidence Guided Search-Space Reduction for Exact DAG Learning

Decentralized Proximal Stochastic Gradient Langevin Dynamics

SVI-DAG: A Structured Variational Inference Approach to Bayesian Causal Discovery