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

Directional Evidence Guided Search-Space Reduction for Exact DAG Learning
Oct 6, 2026Upala Junaida Islam, Abdelmonem Elrefaey, Rong PanScore-Based Causal Discovery Bayesian Network Structure Learning
School of Computing and Augmented Intelligence, Arizona State University, Tempe, AZ 85281, USA.
Decentralized Proximal Stochastic Gradient Langevin Dynamics
May 1, 2026Mohammad Rafiqul Islam, Lingjiong ZhuMarkov Chain Monte Carlo Stochastic Gradient Langevin Dynamics
Department of Mathematics, Florida State University, Tallahassee, Florida, United States of America
SVI-DAG: A Structured Variational Inference Approach to Bayesian Causal Discovery
Aug 5, 2026Shrenik ZinageCausal Discovery Variational Inference
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.

Oct 6, 2026J/K move · Enter open · S save

Directional Evidence Guided Search-Space Reduction for Exact DAG Learning

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

Open paper AskPDF arXiv

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.

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