Expressiveness of Parametrized Distributions over DAGs for Causal Discovery
Authors: Simon Rittel
Sebastian Tschiatschek
Venue: Accepted by TMLR ()

Abstract: Bayesian approaches for causal discovery can in principle quantify uncertainty in the prediction of the underlying causal structure, typically modeled by a directed acyclic graph (DAG). Various semi-implicit models for parametrized distributions over DAGs have been proposed, but their limitations have not been studied thoroughly. In this work, we focus on the expressiveness of parametrized distributions over DAGs in the context of causal structure learning and show several limitations of candidate models in a theoretical analysis and validate them empirically in supervised settings. To overcome them, we propose using mixture models of the considered distributions over DAGs.

Published in Transactions on Machine Learning Research (10/2025)
Expressiveness of Parametrized Distributions over DAGs for Causal Discovery
Simon Rittel
imonnrittel@lmu.cd
Department of Statistics, LMU Munich, Germany Munich Center for Machine Learning, Germany UniVie Doctoral School Computer Science, Austria
Scbastian Tschiatschek
scbusskan.tschiaasccceeee.
Paculty of Computer Science, University of Vienna, Austria Reseorch Network Data Science, University of Vienma, Austria
Reviewed on OpenReview: /ttp::///penrevve// net/ forum? JJON6JRR
Abstract
Bayesian approaches to causal discovery can-in caaneiple-quantine uncertainty in the prediction of the underlying causal structure, typically modeled by a directed acyelic graph (DAG). Various semi-implicit models for parametrized distributions over DAGs have been proposed, but their limitations have not been studied thoroughly. In this work, focus on the expressiveness of parametrized distributions over DAGs in the context of causal structure learning and show several limitalions of candidate models in theoretical analysis and validate them empirically in relevant supervised settings. To overcome these limitations, we propose using mixture models of the considered distributions over DAGs, demonstrating improved flexibility and performance.
1 Introduction
Causal discovery, also known causol structure learning (CSL), is the task of uncovering cause and effect relations among modeled variables based on observed data Glyyoor el. al. 2019 Squires & Uhler 2022) The inferred structures govern the translation of a causal estimand into statistical one that can be measured from data and, hence, provide a basis for causal inference (Pearl 2009 Lundberg et al. 2021). Errors in causal discovery that occur, e.g., due to limited data availability as well as modeling choices, can thus imply different statistical estimands and bias the analysis of causal queries. This dependence highlights the importance of quantifying the uncertainty of inferred causal structures for causal inference.
To this end, Bayesian CSL goes beyond identifying only single, potentially incorrect, causal structure and rather models the uncertainty over the true causal structure. Tt incorporates prior knowledge in the form of a prior distribution and in principle enables the computation of a posterior distribution over possible causal structures when new evidence is presented. Typical assumptions for causal discovery are causal sufficiency, the absence of selection bias, and acyclicity of the causal graph, allowing to represent. the causal structure by a directed acyclic graph (DAG). Given asymptotic identifiability of the causal graph, two sources of crror for its prediction can be distinguished: (i) The apprecimation error results from the finite amount of evidence for the estimation of the possibly nonlinear relations among the observed variables and decreases with increasing data size; (ii) the model error results from the choices for the functional relationships among variables and the distribution over the causal graph and cannot be overcome with more data.
In our work, we focus on the model error and investigate the importance of the model choice for the distri- bution over causal structures. This is particularly important as novel Bayesian CSL algorithms introduced in recent works often differ in multiple model choices and only rarely identify for all of them which key detail is responsible for the reported increase in performance in comparison to competing models. While branch

---

Published in Transactions on Machine Learning Research (10/2025)
of research focuses exclusively on the functional relationships and parametric forms that allow identifica-
tion of cause and effect pairs (Shimizu et al. 2006; Hoyer et. al. 2008 Zhang & Hyvärinen 2009 Loh &
Bühhmann 2014 Immer et al. 2023). distributions over causal graphs, e.g., DAGs, have so far received only little attention.
Contributions In our work; we investigate and compare the expressiveness of distributions over DAGs
used in recent Bayesian CSL algorithms (Cundy et al. 2021; Lorch et al. [021] Charpentier et al.) 2022;
Deleu et al. 2022 Rittel & Tschiatschek 2023 Annadani et al, 2023; Toth et al. 2024). We highlight
their limitations in assigning equal probabilities to graphs of the same Markov equivalence class (Example 1 as well as capturing dependencies between edges in the graphs (Example 2) and compare their ability to match synthelic graph distributions. In addition, we provide experimental and theoretical evidence that the probabilistic models proposed in [Riteel & Tschiatschek/ (2023); Deleu et al. (2022); Toth et al. (2024) and mixtures of them are more expressive than simpler particle distributions. We believe that our work will be helpful to researchers and practitioners alike by demonstrating shortcomings of recently proposed distributions over DAGs that can limit the applicability of Bayesian CST and proposing mixture models effective countermeasure. In addition, we provide a recipe for efficiently approximating the probability mass function of given graph, directed path, or subgraph induced by generative model for the causal DAG by marginalization over an auxiliary random structure using importance samples.
Structure of the paper The remainder of this paper is structures as follows. In (section 2 we introduce our notation and background on causal discovery. We proceed in |section 3 with the presentation of the considered distributions over DAGs and discuss their theoretical limitations. In section 4 we outline efficient method to evaluate generative models based on importance sampling and derive analytically the minimal statistical divergences for particle distributions with K graphs. Finally, we apply both findings for supervised learning of different target distributions in section 5 and conclude our paper in section 6
2 Preliminaries
Probability and random variables We denote scalar random variables by lower case symbols, e.g., Y, random vectors and random matrices by upper case bold symbols, e.g., Y. Single elements of a random vector or matrix are written as Y. The probability of a discrete random variable y taking the value y is expressed by the probability mass function Py(v) == Ply = y)- When clear from the context, we omit the random variable in the subscript of a probability mass function. For random variables with parameterized probability mass functions with parameters 0, we write (1)).. With a slight abuse of terminology, we refer to a distribution Py by its induced probability mass function Py (b) Py(t)). To distinguish appraximations from the true target distribution, we denote model distributions by gy(v)= = (ry})) or simply go(v)
Functional causal models A functional causal model (FCM) is defined as a triple MM (X, (E, Pe), f} consisting of a set of endogenous random variables X, a set of exogenous noise variables with joint prob- ability distribution Pe and set of deterministic functions f, all three indexed by [D] (1....,)) D). Each endogenous variable X2 of the model is generated as a function of a subset of the endogenous variables X and its exogenous noise Ed, i.e., X == fa(X, Ea). The distribution of X. is implicitly defined as the push- forward mcasure of Pe through the causal mechanism fa. The structure induced by the dircct functional dependencies is often restricted to be acyclic such that it can be represented by a directed acyclic graph DAG or eqiivivvllently its adjacency matrix G E c A = {0, DxD with a one-to-one correspondence between random variables and nodes. Let rd denote the index set [D] \ {d), then the d-th column of G encodes the parents (Xa) C X-d of node/random variable Xa, i.e., the subset of X that has direct influence on Xa via fa. Ancestors of a random variable, (Xa( have directed path to X4 in the causal graph G. If they are not themselves parents of Xs, then their cansal effect on X is mediated by at least one parent of X.. Children CXa( and descendants Dea(X.) are causally affected by changes of the value of the corresponding node Xa and are defined as inverse relations to parents and ancestors.
Causal structure learning The objective of causal structure learning (CSL) is to infer the underlying causal graph G from observed random variables X that encodes the causal effects implied by the FCM
2

---

Published in Transactions on Machine Learning Research (10/2025)
(G
G)
G
(G) G(G)
D
G(
X1
X2
X3
G
X1
X2
X3
(a)
(b)
(c)
G
X
Figure 1: Modeling errors in Bayesian causal discovery due to approximations g of (a) the marginal likelihood function and the model for (b) the prior and (c) the posterior over the causal DAG G. For different graphs they either imply undorestimation or overestimation.
X1 X2 X3
Figure 2: Graphs of the MEC from Example and its CPDAG (below).
Mx- Throughout this work, we assume causal sufficiency, i.e., all endogenous variables X are observable, and that the cxogenous noise variables e are mutually independent. This implics that all dependencies and independenicies between the observed values of the random variables X result from their causal effects over the functions f and not from some latent common causes, i.., an unobserved shared ancestors of them. In addition, assume that all samples in the data set D (X J==I are generated i.i.d. from the FCM Mx without any selection bias for the generated samples, c.g., there is no conditioning on unobserved confounders.
Assuming faithfulness, the Markov equivalence class (MEC) of causal graph G' can be consistently identified from D using conditional independence (CT) tests (Spirtes et al.) 2000). This equivalence class contains all DAGs that entail the same observed independence relations and can be compactly represented by a completely partially directed acyelic graph (CPDAG), an acyelic mixed graph that has the same adjacencies, X, Xx, and unshielded colliders, X, - Xx < -j,, as the underlying true graph G, but with some of its edges remaining undirected. While CI tests can be easily parallelized and their required overall number can be sequentially restricted by the individual test results as in the PC algorithm (Spirtes et al. 2000), the combination of the uncertainty attached to each test is non-trivial. Moreover, CI testing for continuous random variables lacks statistical power against alternatives and suffers from the curse of dimensionality (Shah & Peters 2020
Bayesian causal structure learning An alternative to constraint- CSL with CI tests are score- based algorithms that optimize a scalar quantity that is typically derived from the likelihood of the observed data pe (D))). Bayesian causal structure learning is also based on the likelihood but models the full data- generating process for the joint distribution:
N xx mn o, G)((O)
II
(1)
that includes the uncertainty over model parameters θ well as prior distribution over plausible causal graphs GG). The posterior distribution over the causal p(G[D) that quantifies the uncertainty over the true causal graph is then obtained by Bayes' formula:
DG = D(D)), do P(G)
(2)
p(D)
where the so-called evidence p(D) results from marginalization of Eq. over O and G. Averaging the full likelihood p(D(@, G) over the (conditional) prior distribution of model parameters yields the marginal likelihood denoted by pDD(G). In contrast to a plain maximum likelihood estimate pewn. (DIG), it avoids overfitting to the noise of the data (Koller & Friedman, 2009) and motivates the stochastic modeling of O. In general, computing the evidence involves marginalization over all DAGs, G E G, and parameters of the likelihood function, + E RR, and is intractable. Therefore, approximate methods for Bayesian inference are
3

---

Published in Transactions on Machine Learning Research (10/2025)
needed. In the following, we focus variational inference where the joint distribution D) is approxi- mated by generative model (G))), and variational family of distributions Q = {1 (G & E }P)) is specified for the intractable posterior p(G(D). Both components involve parametric distributions over DAGs.
Bayesian model error Given a fixed data set D, deviations of a modeled posterior (G(D) to the true posterior p(G can be backtracked to three different model errors as depicted in Figure (a) The modeled marginal likelihood (D[G) may induce a bias, with the particular case of the maximum likelihood estimate leading an overconfident prediction. Throughout this work, we assume a flexible model for the marginal likelihood, and focus on the role of the distribution over DAGs instead.
(b) The model for the prior distribution (G) directly constrains the belief that domain experts can express over the causal structure, but the subjective bias vanishes in the asymptotic limit as long as it assigns cach DAG some positive probability mass. (c) Lastly, the true posterior may not be in the variational family 2, i.e., it cannot be expressed by (G
3 Distributions over DAGs
In this section we discuss different possibilities for representing distributions over DAGs. While the number of DAGs for a given number of variables D is finite, it is super-exponential in D OEIS Foundation Inc. 2023: Stanley 1973). Consequently, using categorical distributions that allow to distribute the probability mass arbitrarily among all graphs quickly becomes infeasible, even for small D. However, for many applications, it is not necessary to specify arbitrary probabilities for all DAGs, i.e., a smaller degree of freedom may suffice. For causal discovery, the posterior distribution expected concentrate on graphs that are similar to the ground truth graph according to which the observed data was generated. This motivates the development of probabilistic models for graphs that are ficxible cnough to model any possible DAG with D nodes and to additionally distribute the probability mass among some candidate graphs but require substantially fewer parameters than a full categorical distribution over all G E G.
Definition 1 (Expressiveness of parametric families of distributions). Let Q be family of distributions parametrized by & and let p be some target distribution over the same sample space Q. Then the family Q is (D, A) exprresiie for a distribution p, if there exists q E Q that has a statistical divergence D equal or less than 2 0, i.c., m (gg llp) A.
In the following, we present different models for distributions over DAGs proposed in recent works and investigate their expressiveness. For clarity and comprehensibiiiity, we introduce these distributions by their generative model and provide the corresponding probability mass functions in Appendix A) As summary, we visualize their generative models in Figure 3 and provide a overview of the figures, equations and the number of learnable parameters for each model in (Table 1 To avoid limitations by specifie functional relationships or subjective prior knowledge, consider independence relations general desideratum. In particular, begin the investigation of the expressiveness of the distributions over DAGs with the following MEC which will serve as rumning example.
Example 1. For parametric linear model with exogenous Gaussian noise, the true causal graph is iden- tifiable only up to its MEC (Peters et al. 2017). Consider the simple case with three variables and the chain graph X1 X2 Xx as the MEC. Its edges cannot be oriented even in the limit of infinite data. To represent the corresponding uncertainty over the true DAG, all graphs of this MEC depicted in Figure 2 the common cause, ( and the two causal chains, (2) and G13), should be assigned the probability 5.
3.1 Independent edges
The arguably simplest probabilistic model for distribution over directed graph A E A not necessarily of a product over independent Bernoulli probabilities, each of them modeling a possible directed edge X, + Xy:
A ~ (A) with Aig gon, (Ais)
(3)
4

---

Published in Transactions on Machine Learning Research (10/2025)
Table 1: Overview of the different considered distributions over DAGs alongside their number of learnable parameters which depends on the number of variables D, particles K, hidden neurons Hv as well embed- ding and key size HE and Hk (their used default values are reported in section D.2
Graph model
Figure
Equations
Number of learnable parameters
Generative model
Distribution
ProDAG
3a
22
D - 1)
Binary adjaceney matrix
36
51
23
D - 1)
Graph particles
3c
6
[24]
K (D(D - +)+
DPM-DAG
[8]
26
D D(D
ARCO-DAG
[S1
266
Hy + 1) (Hy + +1D++ DDDDDDDD
Mixture models
3e
[8] (
25
K[M +1)
GFlowNet-DAG
[20]
2DH + 7(H, 16[HeHH Hi)) + 4HH
where e are the parameters of the Bernoulli distributions over the edges. While self-loops can be directly ruled out by setting Vi e E [D] A. = O, the random graphs A drawn from this distribution can still have cycles. While rejection sampling can be used to constrain the outcome space to DAGs G, it is highly sample- inefficient. Modeling a random weighted adjacency matrix W E W = R instead of a binary one allows to define a generative model over unique projection of continuous samples to g (Thompson el al. 2024).
W ~ g W)
with Wij ~ qo., (Wi)),
G = P (W),
(4)
that is depicted in Figure& 3a It is based on the idea of a continuous relaxation of the discrete optimization problem over 9 and enforees acyclicity of a single weighted adjacency matrix using a nonnegative, differ- zero for any acyclic graph and to some posit.ive
entiable constraint function h W + RJ that evaluates value quantifying the deviation, e.g., number of closed cycles (Zheng et al.) 2018 Yu et al. 2019 Bello et al. 2022 Nazarct et al. 2024) otherwise. To avoid sampling only fully connected DAGs, thresholding by parameter ò can be applied element-wise to the projected continuous adjacency matrix subsequently. Since the resulting projection Proja,s has to be computed for each sample W, the probability of some target graph under this model cannot be evaluated directly, but estimated using Monte Carlo sampling.
A probabilistic model that allows for analytic odds to compare the probability of two graphs and avoids rejection sampling, can be derived by introducing h within an exponential prefactor to the unconstrained
probability q (A) of Eq. (3) i.e.,
2 (G) x exp(->n(G)) gG)
For a sufficiently high prefactor A, the resulting distribution (e.x (G) assigns only negligible probability mass to any cyclic graph. We provide the corresponding graphical model in Figure 3a For an independent factorization over the edges of the graph, this comes at the cost of expressiveness, since the resulting model can only learn single ordering of the nodes (Rittel & Tschiatschek 2024). In Example 1 each chain implies total order and the common cause partial order, all of them being incompatible with each other. Hence, such distribution with independent edge probabilities can only concentrate its probability mass on one of the three graphs resulting in a skewed uncertainty measure.
3.2 Particle distribution
Lorch el. al. (2021) circunvent the limitation of being constrained to a single order by modeling the graph posterior p( G by a particle distribution consisting of set K particles, each representing a single DAG ((4), The term "particle distribution" originates from the idea of approximating a continuous density function by discrete probability masses. In the context of approximating discrete distribution, it rather refers to simplified model that constrains the support to fewer possible outcomes, e.g., |9] K DAGs. The corresponding generative model is visualized in Figure 3d and can be expressed
k qw (k),
G G (k)
(6)
where qu is categorical distribution with K possible outcomes and weights w defining the probability of each particle/DAG.
5

---

Published in Transactions on Machine Learning Research (10/2025)
+
A
G
G
G
(a)
(1)
Figure 3: Generative models of the considered candidate distributions over DAGs in section 3
(a) Projec-
tion of sampled continuous adjacency matrix W, (b) constrained binary random matrix A, (c) K graph distribution extending the former model by a categorical distribution over k matrices, (d) order-based model that masks A by some random pernutation matrix II, (e) K order-based model, (f) sequence- based model that adds edges autoregressively until the finalization signal f is sampled or the DAG is complete.
3.3 Mixture models
More expressive distributions over DAGs can be defined by replacing the K particles from section 3.2, each corresponding to specific graph, by K distributions over DAGs, i.e.,
~ Iw (k),
G ~ (a() (G).
(7)
Consoquently, we no longer refer to the resulting distribution particle distribution but as a probabilistic mixture model.
In principle, Eq. (5) could be used as the base model allowing to model K different total orders the nodes instead of single [ In the running Example 1 the stated MEC consists only of graphs, but each of them is admissible under difffernt total order. This illustrates that here at least K = 3 required. In case the distribution concentrates its probability mass graphs that are admissible for low number of orders, e.g., sparse graphs with multiple components, such model may suffice. Although smaller than the nunber of possible DAGs, the number of total orders and permutations D! is still super-esponential Hence, for a higher number variables only a very small fraction of all permutations can be represented by the resulting mixture model. Evidently, having a single probabilistic model for each total order constitutes a poor trade- off between expressiveness and efficiency due the increased number of trainable model parameters. To address this limitation, we consider probabilistic models that explicitly allow for different orderings in the following without enumcrating all of them.
3.4 Probabilistic models over orders
Several works apply the idea of an order-based search Friedman & Koller 2003 Teyssier & Koller 2005
in probabilistic generative model (Cundy al.] 2021 Charpentier et al. 2022 Wang et al. 2022 Rittel
& Tschiatschek] 2023] Annadani et al.) 2023).. Their shared underlying idea is learn total order of the variables that imposes acyclicity constraint the (sampled) adjacency matrix A. Following the description of DPM-DAG (Rittel & Tschiatschek 2023). we outline its generative model and provide the corresponding graphical model in Figure 3d The generative process starts by drawing total order that defines a permutation II of the variables. A permuted upper-triangular matrix of ones M then defines a random acyclicity matrix (TT) that is used to mask a sample of an unconstrained adjacency matrix as modeled in Eqq. (5) i.e.,
~ go (I),
A A (A),
G MII A
(8)
M
This masking by (T) can be interpreted probabilistic projection of A to G. In contrast to Eq. (4) gp.o(G can be efficiently evaluated as II and A are independent random matrices with analytic probability
This correspond the model Lorch al. 2021 without collapsing the K distributions K single graphs.

---

Published in Transactions on Machine Learning Research (10/2025)
X1
X2
X
X2
X1
X2
X
X2
X
X
X3
X
(a) GO)
(b) G
(c) ((a)
(d) (
Figure 4: Graphs of Example 2 that all share the highlighted edge between X₃ and X2.
mass functions. For details we refer to soction A.6 The distribution over the permutation can be modeled by the Plackett-Luce (PL) distribution (Plackett 1975), but comes with a limitation for causal discovery (Rittel & Tschiatschek 2023; Toth et al.) 2024) that can be demonstrated with Example 1 Both chains, GO2 and G imply total order of the three variables, but reversed ones. In the PT, model, their permutations can only receive the same probability mass in the case of uniform weights for all three variables, yet the probability for the two causal chains is then upper-bounded by & and cannot take the value of J. Some remaining probability mass then concentrated on other graphs that do not belong to the MEC. This motivated the ARCO-DAG model (Toth ct al. 2024), an autoregressive model over causal orders that computes the weights for the categorical sampling without replacement at each sampling stage conditionally the previous drawn sequence:
qu (T)
Pepieo (I(d) d=1
II
with 2b(() = fo ([T(1))
(9)
The function fw RDD + [D] takes permutation matrix where some rows are still zeros input and predicts the weights for sampling the next variable in the order. For their experiments, the authors of ARCO applied a multilayer perceptron with single hidden layer with H > D neurons. To avoid marginalizing over the set of possible parents under a sampled permutation as in Toth et al. (2024), the same distribution over an unconstrained graph may be used yielding Eq. (8) but with an autoregressive distribution over the permutation. In the case of Example where the probability mass should be split equally on all graphs of the MEC X, X2 X3, this is consistent with setting PAi3(1) PA.: (1) 0 and PA., (1) = 1 otherwise. To further increase the expressiveness of DPM- and ARCO-DAG, mixtures of both be considered. The resulting graphical model is depicted in Figure 3c
3.5 Fully autoregressive model
In all covered candidate models for distributions over DAGs, the unmasked edges are still sampled from independent Bernoulli distributions as stated in Eq. (5) or relaxations of them, e.g., Gunbel-Softmax dis- tribution Maddison et al. (2017); Jang et al. (2017) Besides acyclicity, the edges of the causal graph drawn from these distributions are not coupled limiting the models' expressiveness.
Example 2 (Coupled odges). Consider the four DAGs depicted in Figure with G(I) being the truc causal graph. Assune the functional relationship between X, and X₂ is identifiable and, hence, the cumulative probability of all graphs that contain this edge should be very high in the posterior. If the structural equation for variable X3, (X,, c3), depends only both variables X, and X2, but not a single one alone, c.g., fs X,X2 Cs. then the corresponding edges are coupled. Consequently, the likelihoods of 1+X7+X3 G(2), (3)) and ((4) are almost the same, but all lower than the one of the true graph G(1), Favoring sparser graphs, e.g., by the prior, can then lead to posterior that concentrates its probability mass on (() and (2) that misses X₃ or X3 as a single cause of X2. The graphs (3) and (()) with only X1 or X, as causes of X3 should be assigned small probability in this setting.
To enable arbitrary dependencies between edges, Deleu et al. (002)) propose to construct DAGs sequentially by adding edges to DAGs with fewer edges, starting with the empty graph Go. Their generative model, to which we refer to as GFlowNeetDAG. is motivated by the following rationale: Denote E := |G| E No as the number of edges of DAG. Then a sampled (possibly empty) sequence of distinct edges S (Si...., Sg) E S uniquely defines DAG G over mapping 9
G, although E! different sequences yield the same

---

Published in Transactions on Machine Learning Research (10/2025)
graph G. The underlying sampling of the GFlowNet-DAG model is based transformer architecture -> - of the generative process as well as the probabilities for the potential new edges, == (S Each construction step i begins with sampling of a binary signal f, that indicates whether the graph prototype G,-1 is final, [ 1, and the generative process stops. Otherwise, == 0, a new edge S, is sampled and added to the current graph Ge
ls G
RDD that autoregressively computes the parameters defining the probability for discontinuation
f ~ (f S, ~ (o( G, = g((:)
(10)
During the construction of the DAG, acyclicity is enforced by masking any edges that would create cycle, i.e, the transitive closure of the adjacency matrix of the graph that decodes ancestral relations. After the addition of new edge, the ancestral mask is updated and the iterative sampling continues until the stop signal f; = 1 is sampled or the DAG is fully connected. The corresponding graphical model is shown in Figure 31, The discontinuation signal fi is a necessary requirement, since it guarantees that each DAG can be assigned an arbitrary probability mass provided that the function to has the capacity to model the transition probabilities between different states exactly.
4 Evaluation
In contrast to Bayesian causal discovery algorithms where an observed data set the basis for unsupervised learning of the causal graph, in this work we focus on the comparison of different models for distributions over DAGs w.r.t. their expressiveness for known target distributions as outlined in Def. I, This supervised setting allows us to shield off several sources of error, including any bias due the choice of functional relationships, prior distributions over the data, the size of the observed training data.
To demonstrate the limitations of the candidate models for distribution over DAGs from section 3 we fit each candidate model (G to a specified target distribution Pg- The target distribution is either derived from the MEC in Example 1, by the coupling of edges from Example 2 or a synthetically generated distribution that arises from concentrating the probability mass around a larget graph based on the structural Hamming distance (SHD).
For training of the parameters & with gradient descent, take the forward Kullback-Leibler (KL) divergence between the larget distribution pG and the candidate distribution qa the loss function and approximate it using samples from the target distribution. Due to the limited support of particle distributions, we evaluate the fitted candidate distribution with the reverse KL divergence (KG (lP)) and the total variation distance Drv(ugllas (UGlle))
In contrast to its backward formulation, the forward KL divergence estimates the log- ligelikelihood ratio using samples from the target distribution. During supervised training with gradient descent, this ensures finding the region of the support where the probability mass is concentrated. To illustrate this important detail consider the GFlowNet model where the generation of a DAG starts with the empty graph. If the probability mass is concentrated around graphs with a moderate number of edges, almost no signal for learning is provided for graphs with only one or two of its edges.
4.1 Particle distributions
For particle distributions, we can compute their optimal statistical divergences general discrete distribu- tions over finite of outcomes.
Lemma (Minimal statistical distances for particle distributions). The reverse Kullback- divergence, the total variation, the Hellinger distance, and the Bhattacharyya distance between discrete target distri- bution PG and a particle representation qG, as defined in Eq. (6) are minimized, if the particles of qa represent the graphs with the highest probability mass in the target distribution pG and have their normalized
8

---

Published in Transactions on Machine Learning Research (10/2025)
probability mass:
K
(G (G) = Σ k=1
G = ((k)
(c(k))
(11)
p G((())
with Ak E (1....,K} K}:
G(h)
arg max GEAX
P(G),
(G
(12)
Σ==ı P(G())
The minimal reverse KL divergence and total variation distance are given by:
K
min PRI.(allPcc) = log
Σ
p
((())
(13)
qEQ
k=1
K
min Drv Prolllll =1 Σ
p(G
(14)
gEC
We present the full derivation of all minima and the proof of Lemma in Appendix B
4.2 Semi-implicit generative models
A latent generative model for G defines sampling procedure that starts with sampling some auxiliary random objects Y ~ 9Φ. A function then maps the drawn random samples Y to DAG, i.c., G G= oYY), and induces a joint distribution for gGG, Y). For DAGs, Y can for instance be continuous adjacency matrix as in Eq. (4) a combination of a permutation and an unconstrained binary adjacency matrix as in EEq (8) or a sequence of edges S as Eq. (10) For such an implicit distribution, the probability of some given graph G*, c.g., sampled from the target distribution, cannot be casily evaluated since g is in general not injective. The joint distribution has be marginalized either by analytic iengration/explicee summation over Y
approximated by samples of Y, i..,
4s (G')) ((-
dY
(15)
Direct sampling of G over has relatively high variance for graphs with low probability. Moreover, op- timization of the parameters of a candidate distribution by gradient descent with samples from the target distribution p is costly, since the update requires computing the probability for all drawn samples from the candidate distribution 9. As remedy, we propose use importance sampling to efficiently estimate the marginal probability of G" using discrete samples Y from the proposal distribution o that only generate the target graph G*
Go Y*) Cg Y.)
2 (G))
Σ
-
(16)
Y":g("+)=G-
Lemma 2 (Sample-efficinecy of the IS estimator). Assume that q( Y) is defined by generative model which produces random DAG G a function g for some sampled discrete variables Y and partitions the space of DAGs. The imporlance-sarpll eslimalor @is(G*) derived from Eg. (16) allows direclly drawing samples of G and has smaller variance than standard Monte Carlo sampling, i.e., its effective sample size is greater than N, the number of drawn samples.
Note that continuous adjaceney matrices do not directly partition the discrete space of DAGs allowing to analytically evaluate the probability of a target graph. Hence, we focus investigation order- and sequenced-basee based models where Y is discrete and the optimal proposal distribution & be derived directly from the distribution i6 by constraining Y to Y*, c.g., evaluating only admissible permutations and unmasked parts of an adjacency matrix. Due to the parameter sharing of Φ, the proposal distribution is not constant over training and is derived from the candidate (target) distribution at every optimization step. We provide the proof of Lemma 2 and details on the implementation in Appendix CJ While we present the proposed importance sampling evaluate the probability of a full graphs, its application directed paths subgraphs is ftrwig.
9

---

Published in Transactions on Machine Learning Research (10/2025)
5 Experiments
In our experiments, compare the following models from section 3 and Table 1
in supervised sctting:
K graph porticle distributions. Distributions over a fixed number of DAGs.
Reinforced probabilistically masked DAGs (RPM-DAG). A discretized version of DPM-DAG (Rittel
& TTsciiatcchkk 2023) that we refer to as RPM-DAG, since the score-function gradient estimator is applied instead of pathwise Suftee-SoftSort gradient estimation,
ARCO-DAG. An autoregressive model over the causal order (Toth et al. 2024),
Mirlure models of RRMM-DGGs/ARCO-DAAA A categorical distribution over K order-based models as depicted Figure 3e (with K equals the basc model).
DPG. GFlowNet that generates DAGs by the sequential addition of odges starting from the emply graph (Deleu al. 2022).
In the supervised setting, we minimize the forward KT, divergence between the target distribution and the model distribution using the Adam Optimizer with decoupled weight decay (Loshchilov & Hutter (0019) over 1000 optimization steps. Further details and the used hyperparameters for each model are provided in
sect ion D.2
5.1 Markov equivalent graphs
The first experiment consists of the MEC represented by the chain graph A - B C motivated in Example I The target distribution assigns all three graphs of the equivalence class the samne probability, i.e., 3
In Table 2a we report the mean of the reverse KL. divergence and the total variation distance logether with their standard deviations error over 20 independent training and evaluation runs for the RPM-DAG, ARCO-DAG, and GFlowNet-DAG model. For the particle distributions we compute the optimal empirical values according to Eq. (13) [14) and [12)
The results align with the theoretical analysis that RPM-DAG fails to distribute the probability equally among all graphs of the MEC. The observed preference for the common cause ( is model bias arising from the fact that G() induces only a partial order that is compatible with two total ones.
For = 3, a particle distribution can perfectly match the target distribution. For fewer particles, this matching is not possible and the particle distribution is outperformed by ARCO-DAG and GFlowNet-DAG. However, our experimental result for the mixture of 3 RPM-DAG models indicates that these optimal theoretical values are not reached during model training. With an increasing number of components in the mixture models, the statistical divergence only decreases slowly, while their standard deviations increase initially, before they decrease again. The standard errors underline the validity of the reported mean values, since they are by a factor of /20 smaller than the reported standard deviations. The high variance is due to some outliers, illustrating that there is a regime in which the training is not stable.
In contrast, the results for ARCO-DAG demonstrate that can overcome the limitations of RPM-DAG by using a simple autoregressive model for the permutation weights. In comparison, to the GFlowNet-DAG model, yields considerably lower total variation distance, more stable, and is much faster in training and evaluation which can be attributed partly to the reduced number of model parameters (see Table 1l.
Note that the structure consisting of only three variables also appears in graphs with more variables and comes w.l.o.g.. In the case of a graph with C unconnected components, each consisting of a graph of the simple MEC, its MEC contains 3 graphs. Under assumed ideal conditions and a supervised setting, a particle method picks K of these graphs and assigns A as probability instead of 3-C Eq. (13) and (14) then imply high statistical divergences that are not competitive with parametrizod distributions over orders, e.g. RPM-DAG, that scale well for unconnected components.
10

---

Published in Transactions on Machine Learning Research (10/2025)
Table 2: Reverse KL divergence DRL, total variation distance Drv, and graph probabilities g(G) for different candidate models 9c and target distribution Pq- Empirical means and standard deviations are computed from 20 independent runs, except for the particle representations for which their analytical optimal values are computed using Eq. (13) and (14)
(a) Example MEC class X1 X2 X3.
QG PG
(K [[PG)
Drv(gallpo (Pg))llll)) 4
((1)
((Ga)
(Gaa)
0
0
0.3
0.3
0.3
1 graph particle
1.099
0.6
1.
0
2 graph partieles
0.406
0.3
0.5
0.5
0
graph particles
0.3
0.3
RPM-DAG
0.32673
+ 0.00004
0.34377
= 0.00003
0.53318
= 0.00001
0.161445
+ 0,00002
0.16145 + 0.00002
5-RPM-DAG
0.05
1 0.12
0.05
1 0.11
0.36
1. 0.07
0.31
1. 0.06
0.07
0.06
10-RPM-DAG
0.0013 + 0.0009
0.003
= 0.008
0.335
+ 0.007
0.333
+ 0.001
0.331
+ 0.009
ARCO-DAG
0.00370
0.00005
0.0038
+ 0.0003
0.3322
0.0007
0.3320
0.0007
0.3321
+ 0.0006
5-ARCO-DAG
0.0021
0.0006
0.003
- U.004
0,333
0.004
0.332
0.333
+ 0.002
IO-ARCO-DAG
0.0015
0.0006
0.0220
+ 0.0021
0.333
0.002
0.3327
0.0017
0.3324
+ 0.0012
GFlowNet-DAG
0.004 1 0.015
0.02
T 0.04
0.328
T 0.026
0.34
T 0.04
0.328
1 0.023
(b) Example 2
Coupled edges X1 X3 X2 and identifiable causal effect X1 X2.
QG
(ga(ge/Ipc))
Drviga ( 1
((1))
(G22)
(Gs) + ((.))
PG
0
0.3
0.6
0.1
graph particle
0.5108
0.4
U
graph particles
0.10544
0.5
U
3 graph particles
0.0513
0.05
0.3158
0.6316
0.0526
4 graph particles
0.3
0.6
0.1
RPM-DAG
0.333
+ 0.009
0.36
+ 0.01
0.126
+ 0.012
0.414
+ 0.023
0.44
+ 0.02
2-RPM-DAC
0.0024
0.0008
0.0046
+ 0.0047
0.2986
0.0006
0.5970
0.0046
0.1004
1 0.0011
ARCO-DAG
0.3288446
0.000002
0.355704
= 0.000012
0.122402
0.000015
0.42189
0.00003
0.45447
= 0.00002
2-ARCO-DAG
2.00126
0.00005
0.0021
1 0.0017
0.2999
1 0.0020
0.5993
0.0020
0.1004
1 0.0002
GFlowNet-DAG
0,0001
1 0.0003
0.004
f 0.004
0.30
0.05
0.599
1 0.004
0.099 T 0.004
5.2 Coupled edges
Dependent edges in the posterior in Example motivate the second experiment in which consider the corresponding graphs of Figure 4 with the following probabilities. The target distribution concentrates 60% of the probability mass on the graph G that misses X and X₂ as parents of Xs, further 30% are assigned to the true causal graph G(), To represent a coupling the two causes, the graphs (3) and (4)) that miss either of these parents get only a probability of 0.05 each.
In Table 2b the corresponding metrics for 20 runs are listed, following the same reporting as in section 5.1 Both RPM-DAG and ARCO-DAG, assign the graphs G and (4)) too high probability as they fail to account for the coupling of the edges X₁ + X3 and X₂ X3. Due to its autoregressive model over edges, GFlowNet-DAG can approximate the target posterior with very high accuracy in terms of individual graph probabilities as well as the two evaluated statistical divergences. By design, a particle distribution is not constrained by dependent edges. Since 90% the probability is concentrated on two graphs, a reasonable approximation can be obtained by modcling only two graphs in this particular example. However, our results show that both probabilistic mixture models, with two components of RPM- ARCO-DAG models each, outperform distribution with 3 graphs particles.
11

---

Published in Transactions on Machine Learning Research (10/2025)
Table 3: Target distribution concentrated around the MAP graph G.. Graphs with the same SHD, G GSHD, are assigned the same probability.
Table 4: Statistical divergences for different candidate models for the synthetic posterior distribution in Table 3 averaged over 20 runs.
SHD p(GsmD)
GSHD
p(G c GSHD)
(Ra(salip)) Drv(galliG)
0.15
0.15
1 graph particle
1.8971
0.8500
0.425
8
0.05313
10 graph particles
0.5395
0.4170
2
0.225
28
0.00804
25 graph particles
0.3516
0.2964
3
0.1
61
0.00164
50 graph particles 0.1969
0.1787
4
0.05
0.00053 94
100 graph particles
0.1042
0.0989
0.01581
111
0.00014
150 graph particles
0.0751
0.0723
0.01439
101
0.00014
250 graph particles
0.0426
0.0417
7
0.01068
75
0.00014
RPM-DAG
0.133 + 0.008
0.166 1 0.007
8
0.00627
44
0.00014
2-RPM-DAG
0,107 + 0.009
0.149 = 0.009
9
0,00242
17
0,00014
3-RPM-DAG
0.070 + 0.007
0.133 + 0.006
10
0.00043
3
0.00014
4-RPM-DAG
0.061 0.006
0.123 0.006
5-RPM-DAG
0.057 0.005
0.120 0.005
10-RPM-DAG
0.041 1.000
0.09 1 0.01
20-RPM-DAG
0.033 + 0.004
0.072 = 0.006
ARCO-DAG
0.087 0.009
0.146 + 0.007
2-ARCO-DAG
0.010
0.137 0.008
X1
X₂
X4
X3
3-ARCO-DAG
0
0.135 0.008
4-ARCO-DAG
0..09
0.130 = 0.007
3-ARCO-DAG
0.07 0.01
0.130 = 0.007
10-ARCO-DAG
0.058 0.007
0.114 1 0.007
Figure 5: Assumed true graph G₁ representing the MAP graph in the posterior modeled in Table 3
20- AD-ARCDDA
0.006
0.09 = 0.10
GelowNet-DAA
0.25 +
0.25 0.08
5.3 Concentration of posterior mass
For an identifiable FCM and a high number of samples, the posterior should be concentrated on graphs that show high likelihood. Tn the idealized setting of perfect regression of child its parents, the likelihood can be expected to peak for the true underlying graph G provided that parameter uncertainty regularization prevents superfluous edges that are not in Gt. This implies that high-scoring graphs are 'similar' to G2 implied by the FCM.
In the absence of an analytic posterior that motivates such similarity, we generate a synthetic target distri- bution around the assumed muzzimmmum-a-ppsleeiiin (MAP) graph G₁ depicted in FFiiue 5 that has positive support for all 543 possible DAGs with 4 nodes. We start by assigning probability masses to the groups of graphs that have structural Hamming distance (SHD) up to 4, i.c., up to 4 different entries in the adjacency matrix, and split the remaining probability mass equally among all graphs with higher SHD. For the five groups with SHD up to 4, we distribute their cumulated probability masses equally among all graphs within the respective group. The target distribution is summarized in [abll 3 which lists the probability masses for the groups of graphs with the same SHD alongside the individual ones.
We report in Table the mean values and standard deviations for the statistical divergences over 20 in- dependent runs for all candidate models, except the particle representations for which the optimal values are computed analytically using Eq. (13) and [14) The results highlight the expressivity of the RPM-DAG and its extension ARCO-DAG that both outperform particle distribution consisting of the 50 graphs with the highest probability in the target distribution that account for a cumulative probability of 86.3%. A single ARCO-DAG yields lower statistical divergences than a mixture of 2 RPM-DAG models in this setting, When comparing mixture models with K 2 3, the RPM-DAG mixtures show better (mean) results than ARCO-DAG. Although we trained in an idcalized very low-dimensional setting with samples drawn from the target distribution, the good performance of the GFlowNet-DAG model from section 3.1 and 5.2 cannot be observed in our considered setting. Its statistical divergences are much higher than for single RPM
12

---

Published in Transactions on Machine Learning Research (10/2025)
ARCO-DAG model. We conjecture that the gradients from a variety of different graphs do not yield sulficiently stable training signal that is necessary to tune the transformer network of GFlowNet-DAG.
In addition to the results for the expressiveness of the candidate distributions, we report the training times of the candidate models and considered mixtures them in section D.3 Notably, training the GFlowNet-DAC model takes 250 times longer than training the RPM- and ARCO-DAG.
6 Conclusion
Bayesian causal discovery promises uncertainty quantification in the prediction of the causal graph underly- ing the observed data. We reviewed several candidate models for distributions over DAGs and investigatod limitations of their expressivily. To minimize confounding factors such the influence of functional re- lationships the size of the training set, we considered an idealized supervised setting that can also be used to specify prior distributions. We showed in our theoretical analysis and experimental results that all considered candidate models, except for the autoregressive GFlowNet-DAG model, are theoretically not sufficiently expressive to match simple target distributions with coupled edges. While GFlowNet-DAAG is theoretically sufficiently expressive, we cannot confirm its expressiveness in a low-dimensional experiment where graph samples are drawn from the larget distribution with support over all possible DAGs. Since causal structure learning is typically an unsupervised problem, this poses a major limitation of this model.
Coupled edges, due to interaction effects, pose general challenge in causal structure learning and affect most Aigorithns-noo only Bayesian approaches. Our results suggest. that mixl ure of a small number of simple probabilistic models such as the RPM-DAG model may approximate distributions with coupled edges sufficienlly well in practical applications and outperform particle distributions with a moderate number of modeled graphs. Its exlension, ARCO-DAG, is more expressive w.r.t. probabilistic causal orders and shows competitive performance to GelowNet-DAGG while having a substantially lower number of parameters, and is more stable and by magnitude faster in training. When comparing RPM- and ARCO-DAG mixture models, the RPM-DAG mixtures outperform the ARCO-DAG mixtures for 10 mixture components consistently in all experiments. Therefore, we conjecture that mixture of RPM-DAG models is particularly suited to scale Bayesian causal discovery with higher numbers of variables.
References
Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong, BayesDAG: Gradient-basee posterior inference for causal discovery. In Advances in Neural Information Processing System, volume 36, 2023.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: learning DAGs via M-matrices and a log- determinant acyclicity characterization. In Advances in Newral Information Processing Systems, volume 35, 2022.
Bertrand Charpentior, Simon Kibler, and Stephan Günnemann. Diffcrentiable DAG sampling. In 10th International Conference on Learning Representations. Openreview, 2022.
Chris Cundy, Aditya Grover, and Stefano Ermon. BCD nots: Scalable variational approaches for Bayesian causal discovery. In Advances in Neural Information Processing Systems, volume 34, pp- 7095 7110, 2021.
Tristan Deleu, Anlónio Góis, Chris Chinenye Emezue, Mansi Rankawal, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. Bayesian structure learning with generative flow networks. In Conference on Uncer- tainty in Artificial Intelligence, volume 180 of Proceedings of Machine Learning Research, pp. 518-528. PMLR, 2022.
Nir Friedman and Daphne Koller. Being Bayesian about network structure. A Bayesian approach to structure discovery in Bayesian networks. Machine Learning, 50(1-2):95-12 2003.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10:524, 2019.
13

---

Published in Transactions on Machine Learning Research (10/2025)
Patrik O. Hoyer, Dominik Janzing, Joris M. Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, volume 21, pp. 689 696, 2008.
Alexander Immer, Christoph Schultheiss, Julia E. Vogt, Bernhard Schölkopf, Peter Bithlman,, and Alexan- der Marx. On the identifiability and estimation of causal location-scale noise models. In International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 14316- 14332. PMLR, 2023.
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. In 5th International Conference on Learning Representations. Openreview, 2017.
Daphne Koller and Nir Friedman. Probabilistic Graphical Models Principles and Techniques. MIT Press, 2009.
Augustine Kong. A note on importance sampling using standardized weights. University of Chicago, De- partment of Statistics, Technical Report, 348:14, 1992.
Po-Ling Loh and Peter Bühlmann. High-dimensional learning of linear causal networks via inverse covariance estimation. Journal of Machine Learning Research, 15(1):3065 3105, 2014.
Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. DiBS: Differentiable Bayesian struc- ture learning. In Advances in Neural Information Processing Systems, volume 34, pp. 24111 24123, 2021.
llya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In 7th International Conference on Learning Representations. OpenRevicw, 2019.
Tan Lundberg, Rebecca Johnson, and Brandon M Stewart. What is your estimand? Defining the target quantity connects statistical evidence to theory. American Sociological Review, 86(3):532-565, 2021.
Chris J. Maddison, Daniel Tarlow, and Tom Minka. A sampling. In Advances in Neural Information Processing Systems, volume 27, pp. 3086 3094, 2014.
Chris J. Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. In 5th International Conference on Learning Representations. Openreview, 2017.
Achille Nazaret, Justin Hong, Elham Azizi, and David Blei. Stable differentiable causal discovery. International Conference on Machine Learning. OpenReview, 2024.
OFIS Foundation Inc. The On-Line Encyclopedia of Integer Sequences, 2023. Published electronically at http://oeis.org orge//0022///
Judea Pearl. Causality- Cambridge University Press, 2nd edition, 2009.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causul inference: Foundations and learning algorithms. The MIT Press, 2017.
Robin L. Plackett. The analysis of permutations. Journal of the Royal Statistical Society Series C: Applied Statistics, 24(2):193-200 1975.
Sebastian Prillo and Julian Eisenschlos. Softsort: A continuous relaxation for the argsort operator. In International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 7793 7802. PMLR, 2020.
Simon Rittel and Sebastian Tschiatschek. Specifying prior beliefs over DAGs in deep Bayesian causal struc- ture learning. In European Conference on Artificial Intelligence, volume 372 of Frontiers in Artificial Intelligence and Applications, pp. 1962-1966 IOS Press, 2023.
Simon Rittel and Sebastian Tschiatschek. On differentiable Bayesian causal structure learning. In UAI Workshop Causal Inference, 2024.
14

---

Published in Transactions on Machine Learning Research (10/2025)
Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. The Annals of Slatislics, (8(3):1514-1536 2020.
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti J. Kerminen. A linear non- non-Gaussian acyclic model for causal discovery. In Journal of Machine Learning Research, volume 7, pp. 2003 20300 JML.R, 2006.
Peter Spirtes, Clark Glymour, and Richard Scheines. Cousation, Prediction and Search. Adaptive Compu- tation and Machine learning. MIT Press, 2nd edition, 2000.
Chandler Squires and Caroline Uhler. Causal structure learning: A combinatorial perspective. Foundations of Compulational Malhemalics, pp. 1 35, 2022.
Richard P. Stanley. Acyclic orientations of graphs. Discrete Mathematics, 5(2):171-178, 1973.
Mare Teyssier and Daphne Koller. Ordering-baae bassdd search: simple and effective algorithm for learning Bayesian networks. In Proceedings of the 21sl Conference Uncertainty Artificial Intelligence, pp- 548 549. AUAI Press, 2005.
Ryan Thompson, Edwin V. Bonilla, and Robert Kohn. Prodag: Projection-induced variational inference for directed acyclic graphs. orXiv:2405.155 15167, 2024.
Christian Toth, Christian Knoll, Franz Pernkopf, and Robert Peharz. Effective Bayesian causal inference via structural marginalisation and autoregressive orders. In ICML Workshop on Structured Probabilistic Inference & Generalive Modeling, 2024.
Benjic Wang, Matthew R Wicker, and Marta Kwiatkowska. Tractable uncertainty for structure learning. In International Conference on Machine Learning, pp. 23131-23152 2022.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 7154-7163. PMLR, 2019.
Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlincar causal model. In Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence, pp. 647 655. AUAI Press, 2009.
Xun Zheng, Bryon Aragam, Pradeop Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31, pp. 9492-95440 2018.
15

---

Published in Transactions on Machine Learning Research (10/2025)
A Probability mass functions over DAGs
For conciscness, we presented in section 3 only the generative models for DAGs that induce probability distribution. In the following, we provide their analytic probability mass functions.
A.1 Unconstrained, binary adjacency matrix with independent edges
Sampling ofl-diagonal adjacency matrix A with independent edges as in Eq. (3) does not guarantee that the sampled graph is a DAG. Rejecting all acyclic graphs implicitly defines distribution whose support is limited to g. Even for moderate values of individual edge probabilitics IA,, rejection sampling can be highly sample-inefieient, more cycles. To evaluate any given DAG G, renormalization of the probability of qa (G) by the constant Z is required which can be approximated by the relative frequency that sampled graph is a DAG:
particular for graphs with many nodes, i.e., higher values of D, that allow for
D II
8A (A)
aAu (Aij), i
(17)
1G GG) Z (A((G)
with == Σ [A (A) = EA [A AEA
(18)
However, using the unnormalized distributions QA allows to evaluate the odds of two DAGs without any sampling or computation of normalization constant Z.
A.2 Projection of samples from an unconstrained, continuous adjacency matrix
To avoid costly rejections of cyclic graphs in the generative process, a continuous adjacency matrix W drawn from the generative model defined in Eq. (4) can be projected on the space of acyclic, weighted adjacency matrices V {W W ((W)) = 0}:
V = (W) = arg min IW VI
(19)
Vev
As the corresponding weighted DAG V E V is still fully connected, L1 regularization can be applied to model sparser graphs. Finally, to obtain a binary DAG G Reg g from V, edge-wise thresholding denoted by the function is has to applied in a final step:
G = =Projn, (W) = t5 (Proj,,
(20)
Due to the implicit correspondence between W and G over the resulting projection Proja,s: the sample- efficient importance sampling as proposed in Appendix cannot be applied. This implies a higher cost of ovaluating the probability of single DAG G. Its relative frequency has to be approximated using Monte
Carlo sampling:
Tw (W)
II to., (Wy)), ifj
(21)
qG (G) = /
G Projn (W)) (w( dW Ew [G Proje, WW)
(22)
For further details on the projection, we refer the reader to Thompson al. 2024
A.3 Gibbs distribution
Eq. (5) defines a generative model where the probability for any cyclic graph is negligible small given a sufficiently high prefactor X. As in /eccionn A.I the odds of two graphs can be evaluated even without any sampling, while the caleulation of the probability of a given DAG, ga (G) requires normalization. In contrast.
16

---

Published in Transactions on Machine Learning Research (10/2025)
to Eq. (18) for high value of À samples of the generative model violate acyclicity only with very low probability and the normalization constant 2 is almost 1:
@G (G) 7 exp ( (A
with = Σ
exp (=M((NM))a((()))EE
- 1.
(23)
A-200
GE9
A.4 Particle distribution
A particle distribution for DAGs as defined in Eq. (6) number of DAGs {G()) Jk=1
distributes the probability mass among limited
K GG (G) = Σ
G = G(
qu (x(()))
(24)
k=1
A.5 Mixture distribution
A mixture distribution models K components of given base distribution, e.g,, RPM- or ARCO-DAG. The probability for a given DAG G then results from weighting by their individual probabilities Ic= (G):
qG GG
Σ 96* (k))
(25)
A.6 Probabilistically masked DAG
The probability of a graph G under the generative order-based model defined in Eq. (8) can be obtained by marginalization over the space of all permutations P over D nodes and corresponding upper-triangular matrices (m C A, i.e., the subset of edges that is not masked by the permutation matrix II:
GG (G)
Σ
(-MM o A Qu III)
II +()<(()
(A (At))
with M MIM
(26)
NEP.AEMA
Sampling without replacement from categorical distribution with some fixed weights (val2 is equivalent to drawing a permutation over (D) and known as the Plackett-Lace (PI. distribution (Plackett) 1975). At each sampling stage d, the categorical weights for the remaining variables that were not yet sampled are normalized:
D-1
D-1
4n (T) = II In(d) (=(d)
II v Exed (T()
(27)
d=1
Di
A differentiable, continuous relaxation of the discrete sampling of permutation matrices II can be obtained by pairing the Gumbel-Softtax trick (Maddison et al. 2014) with Softsort. (Prillo & Fisenschlos 2020)
A.7 Autoregressive model over all potential edges
Denoting the mapping between the sequence of edges S E S to a graph G by : S + {0, Dx the probability of sampling a graph G by the generative model defined in Eq. (10) equals:
4G (G) = Σ 9G.s (G,
Σ G 9(S) (S()
(28)
Ses
SES
where the probability of a sequence qs is the product of the autoregressive edge probabilities and probability of stopping the generative process:
E
)
[E-1
(S)) = II
(s. (S.(S:.-1)
II
y (0]5.--1)
(105)-1))te/ssee--- (1/S)
(29)
17

---

Published in Transactions on Machine Learning Research (10/2025)
The sequential process of GFlowNet itself defines DAG where all non-leaf nodes represent states of an unfinished causal graph G = g(..) Starting with the empty graph state the root node and probability 1, it splits the probability of a state representing a preliminary graph G among its children who are either other preliminary graph states with an additional edge or a terminal state, i.. a leave node representing the probability of the DAG G under the GelowNet model Deleu et al. (2022).
B Particle distributions
In the following we present a proof of Lemma 1 and derive the minimal values for the four statistical distances between a discrete target distribution that assigns some positive probability mass to more than K graphs and a distribution with only K particles. the trivial case where p distributes all probability only on K fewer graphs, a particle distribution q can model p perfectly without any approximation error and Lemma follows directly.
B.1 Minimal total variation distance
Proof. The total variation distance Drv between the candidate and target distribution, qG and pc, respec- tively, can be rewritten as:
Prv(gallra) = Σ |GG)
Σ
GG) +
g (((())) - -P(G()))
(30)
GEC
GESi
Σ ((G)))
In the RHS of Eq. (30) the sum over all graphs is split into two terms. The graphs that are not in the support of the particle distribution (G lead to a distance that is independent of the individual probabilities for the graphs. The second term measures their absolute deviations which would be minimized if p(G(k)) = (G(()) Vk [K]. However, due to the normalization constraint, the accumulated probability mass of the not modeled graphs has to be distributed additionally among the K graphs in QG- Consequently, we are interested in minimizing the divergence, (((()) Z ((G(k) yk [K] and thus
K
K
min Drv (PG) min 1 QEQ
1 ---Me (G k=1
Σ ((Gtt))-p(G(*)))
- (G k=1
-1-En(Guy
(31)
q
2
Consoquently, both sums in Eq. (31) amount preciscly to the loss of the accumulated probability mass for the graphs not modeled by 4a and the total variation distance is minimized when (G assigns non-zero probability the K graphs with the highest probability pc:
Vk [K]
(())
arg max PG
(32)
GEIX(G
Their corresponding probabilities (kk))) are not unique in general. The approximate probability distribu- tion qG with the minimal Drv to PG does not constrain how the remaining probability mass distributod among the K modeled graphs. The probability of single graph could account for the missing accumulated probability mass yielding only single biased value, but distorting the relative probabilities. This would limit the utility of such an approximated distribution heavily and motivates imposing the following constraint that is stated in Lemma as premise:
(G(())
p( ((G())
Vijje [
(33)
(GU))
p (G))
18

---

Published in Transactions on Machine Learning Research (10/2025)
Preserving the relative probabilities of PG in qG then uniquely defines dG by the normalized probabilities
((
(G Σj=ı
(34)
and yields an approximation that is for all K graphs overconfident w.r.t. to the approximated distribution PG, i.., (G(()) > ((G()) yk E [K] c [
B.2 Minimal Kullback-Lebble divergence
In contrast to the total variation distance, the Kullback-Leibler divergence not symmetric. Due to the assumed smaller support of the considered particle distribution (G compared to the target distribution over all DAGs PG only the reverse KL divergence, DKL(GG PPG), is well defined.
Proof. The constrained minimization of the KL divergence DKI(sallPG) can solved by using its corre- sponding Lagrangian function C and multiplier y
K
(( p p(G +7
C(gG. 7) Σ (G log
11 (((())
(35)
Setting its derivatives w.r.t. the probability of the particles to zero, reveals that (((()) is proportional to ((G()));
ac (i(G))) = log (G(())
+ +1-7=0
k(G
(((()) = =-- = p((())
(36)
As direct consequence, the normalization constraint leads to the same probabilities as in Eq. (34) that preserves the relative probabilitics pa in (G. Thus the KL divergence simplifics to:
(G)
K
((())) log (G (G
DKL(gGllpa)
Σ
(G)
log
Σ
(GG)
GEG
k=1
K
K
Σ
(G log Σ ((G
log (G
(37)
Note that the negative logarithm is minimized when its argument, the sum of the probability of K graphs in the target distribution Pc, is maximized. Hence, the constrained KL divergence is minimized, when qG
assigns non-zero probabilities to the K graphs with the highest probabilities in pc as in Eq. (12)
B.3 Minimal Hellinger & Bhattacharyya distance
Proof. The derivation of the minimal Hellinger distance Dn follows analogously to the minimal KI divergence in section B.22 The first steps is the derivation of the constraint of preserved relative probabilities as stated
in Eq. (34) from the minimization problem. Since the argument of the outer square root is non-negative and squaring is monotonic function R the distribution g E Q that minimizes the Hellinger distance is preserved when it is squared:
arg min Du ( = arg min Dn ( gEQ
(38)
geQ
19

---

Published in Transactions on Machine Learning Research (10/2025)
Furthermore,
1 Du (nG = Σ
(V[GG) VV(G))
2
GEG
K
(VaaGuA)) - V 2
1
Σ
(GG) + Σ
2
GEO
k=
1
- (V((())))((
2
K
= =1-
(G
(39)
This allows to rewrite the minimization problem as a maximization problem:
K
arg min DH Da'(aalDaeee QEQ
arg max g
a ((
(40)
Then the corresponding Lagrangian with multiplier y reads:
K C(g., 7) = V ((G())))(G()))+ +7 (
k
(41)
Taking its derivative w.r.t (G setting to zero and squaring both sides yields:
dc 1 V(G
-7 = 0
(y( 2 V(()())
(G 4-2 ((
(42)
The normalization constraint then again implies minimal squared Hellinger distance:
Eq. (34)
that can be inserted into Eq. (39) yielding the
K = 1- Σ
p(G (G
K (G
min D (gGllpG) QEQ
=
(43)
EEERit ((G())
A=1
The Bhattacharyya distance Dg is directly related to the squared Hellinger distance D3- Since the logarithm is a monotonic function, the previous result of the maximization problem can be rensed:
DB(gGlIPG) In Σ
Vq(G)p(GG) - In (1 DH (gG
(44)
K
min Dp(9cllpG) In Σ ((G()) t=1
(45)
Hence, the Hellinger distance and the Bhattacharyya distance are both minimized by selecting the graphs with the highest probability in pc as in Eq. (12)
20

---

Published in Transactions on Machine Learning Research (10/2025)
C Importance sampling for parametrized distributions with an auxiliary discrete structure
C.1 Construction of the proposal distribution
Implicit generative models for discrete structures such DAGs can define hierarchical sampling process with discrete auxiliary variables Y that partition the discrete sample space. The weights of discrete objects, e-g.; edges or permulations, that do not align with some larget structure can be set to 0. Since the relative probabilities of admissible substructures are preserved in this case, it is equal in probability to sampling form the unconstrained distribution and rejocting any sample that is not admissible under the target structure. W.l.o.g. we assume for the remainder of the discussion that the target structure a DAG G. Denoting the unnormalized weights for a combination (G, Y) of the candidate model qs by WG,Y,, the proposal distribution g becomes
Lo (G Y") (Y)
WG-.Y Y Ex G = (Y) "Y
(46)
and is well defined YG E g. It constrains the sample space to the target graph G and can be directly derived from the weights w of the candidate distribution Qs at any optimization step:
yΦ G", Y*) E W
WG Y
(47)
Ey
Order-based model For order-based models as described section 3.4 and depicted in Figure 3d the weights of the Bernoulli distribution for any edges that do not appear the target graph can be set to zero, i.e., J with G 0 dby O. In addition, the sampling of total order in Eq. (27) can be constrained to the (partial) order induced by the target DAG by setting the weights of all variables to 0 as long their parents have not been sampled as predecessor. Note that the analytic value of o (G', Y") requires summation of the probabilities of all admissible total orders, in our practical implementation we approximate the sum by finite number of sampled admissible total orders (lower bound). The more edges a DAG has, the fewer admissible permutations exists and the better the approximation by a fixed number of sampled random permutations. In the following we outline the procedure for RPM-DAG with a simple example.
Example 3 (Collider structure). Consider the setting with D 3 and assume the target graph G forms collider, i.e., X, + X3 X2. Then two total orders with X3 as the last. variable are admissible with G", i.., T (123) and IT2 (213).
Under the RPM-DAG model, the probability of G is then the sum of the probabilities of the two options that generate G*:
(G (G') Qur (TT) 9Az2 (0) (1)4A (1) + qo (T22) gAis (1)gaa (1)
(48)
= m1)
=q m²)
Now consider that the proposal distribution generates G by sampling πh Then the resulting joint proba- bility only differs in the probability of sampling πi, since the not masked entries of A are predetermined by the entries of G*:
q" GG" πT) = (044, (1)4Ass (1).
(49)
While on admits sampling of all possible permutations,
4n (T:)
(50)
112 2/ whs 1/2 4/3
th is constrained to ni and πτ. Hence, the normalization of the weights 2o changes according to the restricted set of admissible options:
+
2 (T()) =
(51)
1 1/2
21

---

Published in Transactions on Machine Learning Research (10/2025)
Sequenco-based model Sampling of the next edge in sequence described in section 3.5 and depicted in Figure 31 can be directly constrained to the edges that appear in the target graph. The ancestral masking that ensures acyclicity remains unaffected by importance sampling. The number of edges in a graphs then limits the length of the sequence to consider. In contrast to order-based models, fewer edges imply fewer different sequences that generate the same DAG.
For Example 3 we sketch the implementation for GFlowNet-DAG model. Let denote the two edges of G by G13 (X1 + Xx) and G23 == (X2 Xa). Then the probability of sampling G is the sum of both sequences S3,22 = (G13; G23) and S23,13 = (G23, G13) that generates G*
qG (G' qG.s GG" 513,25) + gG,s(G", S23,13)
(52)
Now assume that the sequence S = S13,23 is sampled by the proposal distribution q². Starting with the empty graph, the probability of stopping has to evaluate to false, i.e, gs(so = 0). After sampling G13, the process continues with fl = 1 and sampling G23, before it stops by f2 1. Note that the probabilities of discontinuation are the same for both, 9G and q:. since they fixed by G The proposal distribution only differs by renormalization of the probability for the first edge, since the second edge the last and there is no alternative anymore. As a consequence, the probability of sampling the sequence S to generate G under the proposal distribution equals
le (G", S) (5,(() (Si) 46 (0) qn (0) (8)))), (1|(S,,
(53)
C.2 Sample-efficiency the proposed importance sampling
In the following, we provide the proof of Lemina 2 by computing the effective sample size Kong (1992 :
Proof. Let's define QMC as the standard Monte Carlo estimator of (G())) based on Eq. (15)
QMC G
1 N
G" s(Y(n)))
(54)
and Òιs the importance sampling estimator based on Eq. (16) with (G+(a), (=(n)) samples from zz:
hss G (G") Σ rs (G+(i),
1 N N
yu(r)) where rG (G) = l0 (G" Y*) q GG', Y+)
(55)
is the likelihood ratio of the larget and proposal distribution that is well-defined by construction VG E g. Both estimators unbiased, but have different variances:
Nar (MC = Ex
G = g{ - Ex E [IG==
= -xxnga
GG" - -(()) 9 ()')²
(56)
NVar [iss( (G*) Br
Fip (G', y³ Eyr~ (To (G", Y²)
(
Σ (G", Y (G', Y*) - 2 (G*)2
==(")"==-]]]] - (G+)²
q
(57)
The derivation of Eq. (57) switches the reference distribution for the expectation value from o to q and is purely notational transtormation. Note that the summation still limited to (G', Y") which explains the added Iverson bracket. Since the parameters & are shared and go (G*) is computed by renormalization over
22

---

Published in Transactions on Machine Learning Research (10/2025)
the constrained space of Y'; have GG E G ra Y*) < 1 and, hence, it is trivial that the effected sample size has to be strictly greater than N:
(S0(4 = =N Var (arl4MG")) G")]
q0 (p()")
N> N.
(58)
Var hi G+)
9Φ (O*("))
@ (G))
Whenever there exists only one admissible Y* that generates G e.g. for RPM-DAG when G implies a total order or for the Gilow-Net-DAG when G* the empty graph, the likelihood ratio r(G", Y*) is constant and equals the very probability of sampling G
ro'"","" Ey G = g W
W Er Y
[us =
(59)
This directly implies then that Var (G)) 0 as a single sample of the proposed importance sampling estimator yields the true value. this case, the estimator Eg. (55) is then optimal w.r.t. to sample-eficieney. In the general case, the proposal distribution in Eq. (46) is not optimal, but sample-efiiciet according to
Eq. (58)
D Details on experimental setup
D.1 Model architectures
RPM-DAG Following the described parametrization of DPM-DAG, we model an off-diagonal adjacency matrix containing the logits of Bernoulli distributions for the directed edges and log-weights for the PL distributions. In contrast to DPM-DAG, we implemented discrete sampling of permutation and adjacency matrices and no continuous relaxations of it.
ARCO-DAG The auloregressive neural network of ARCO-DAG consists of a simple two layer perceptron with H = 30 hidden neurons and ReLU-activations and follows the official implementation on https: //github. ctttuu...om//hhiiiie http://rttooobbearrrrooo
The transformer architecture for the GFlowNet-Dag model follows the official implemen- tation provided on https:: httts:///itthut //github.. com/ caam//risseee trrstandelee// jaauuu/ja-aaaa dag gflownet, as default parameters it uses embedding size of HE 128, a key size of H 64 and Hr, -7 layers of Transformer blocks.
GFlowNet--Dag
D.2 Experimental parameters
In the experiments in section 5.1 and 5.2 evaluate the forward KL divergence in all three, respectfully four, graphs with positive probability mass. The probability of a graph under the candidate model is approximated by [16) using importance samples (IS). To speed up the calculations during training, the number of IS is chosen to be smaller than for the final evaluation of the trained model. The learning rate was chosen by a grid search over {1,5x 11551001110 10-1, 1x 1001 1.5 x 10 2, x 10 10-2,5 x 10 10-3, 9,1x 10 10-3, 5 x 10- 4, 10 *).. Training was performed over 1000 optimization steps to account for the instability of Net We document the used hyperparaaeterss in Tabll 5
Due to the low number of DAGs, first experiment (section 5.1p
allows to evaluate all 25 graphs within single forward pass during training and, hence, optimization step equals one epoch. To account for the higher number of possible sequences, we average the probability of sampling a target graph from the model distribution using 10 importance samples for GFlowNeet-DAG. For evaluation of each graph probability we drew 100 importance samples for all methods.
In the second experiment (soction 5.2), we had limit the number of samples from the target distribution for single optimization step. Due to the high number of trainable parameters we reduced the number of samples from the target distribution from 100 to 25 and used 5 instead of 10 importance samples for GFAG
23

---

Published in Transactions on Machine Learning Research (10/2025)
For the evaluation of each of the 543 DAGs, we decreased the number of importance samples from 100 to 20.
Table Hyperparameters.
(a) For the experiments
section 5.1 and 5.2
Graph model
Learning rate
# forward KL samples
# IS for training
# IS for evaluation
RPM-DAG
0.5
25
100
ARCO-DAG
0.5
25
100
GFlowNet-DAG
0,001
10
100
(b) the experiments
section 5.3
Graph model
Learning rate
# forward KL samples
# IS for training
# IS for ovaluation
RPM-DAG
0,05
100
10
20
ARCO-DAG
0.05
100
10
GFlowNet-DAs
D.3 Training time
In Table we provide qualitative comparison of measured training times of the different graph models for the experiment described in section 5.3 with the hyperparameters listed in Table 5) All models were evaluated over 20 independent runs, each consisting of 1000 optimization steps. The computations were conducted on a 11th Gen. (R)) Core 17-1165G7 processor with 2.80 GHz, 4 cores and 8 logical processing units paired with 32 GB DDR SDRAM. The reported average times for RPM- and ARCO-DAG scale linearly with the number of mixture components, while GlowNet-DAG takes more than 250 times longer than their base models. Its increased computation time is a result of the huge number of parameters of the transformer architecture and the sequential generation of the DAGs.
Table 6: Averaged training time and standard deviations over runs for the probabilistic graph models in the experiment of section 5.3 alongside the number of the models' parameters.
Craph model
# learnable parameters
Training time in min
RPM-DAG
100
0.152 0.005
5-RPM-DAG
505
0.761 + 0.018
10-RPM-DAG
1010
1.51 - 0.04
20-RPM-DAG
2020
3.13 1 0.16
ARCO-DAG
646
0.181 1: 0.004
5-ARCO-DAG
3235
1.855 0.085
10-ARCO-DAG
6470
1.79
+ 0.05
20-ARCO-DAG
12940
3.70
+ 0.03
GFoowNetDDAG
6 176 722
46.4 3.4
24


Reviews:
[{"review_id": "g5MSxUL6Eq", "reviewer": "Reviewer_jnvY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- The topic considered is new and interesting by focusing on expressivity of DAGs\n- Provided illustrative counter‑examples (Examples 1 & 2) that showcase (i) even mass across Markov‑equivalent graphs and (ii) edge‑level dependencies.\n- Lemma 1 derives closed‑form minima for total variation and KL distances\n- Proposed simple fix with empirical gains\n\n**Weaknesses**\n\n- Contribution is essentially a mixture models of distributions. Can you elaborate the novelty over existing literatures of mixture models? \n- There lack details for the proposed mixture method, its computational complexities, and its theoretical guarantees. I am also curious how their theoretical gain beyond baseline methods.\n- Lemma 1 / the main theory and the proposed method is disconnected. Sections 4.1 and 4.2 are also disconnected. Can you elaborate the   utility of results of Lemma 1? \n- The proposed importance-sampling for evaluating implicit models seems very straightforward and lack analysis of variance. Any assumption or condition for the denominator, i.e., the proposal distribution? What if the denominator goes to 0?\n- Assumptions need to be more clearly stated, not blended in the preliminaries of \"causal structure learning\". Do you have any other assumptions? What's \"assumed ideal conditions\" in Section 5?\n- Related work on universal or highly expressive DAG priors—such as DAG-VAE or flow-matching models—is not discussed.\n- While the paper includes results for a few mixture sizes (e.g., 2, 5, 10), it lacks systematic ablation or guidelines on how to choose $K$.\n- Real‑world utility remain untested. Can the author consider some real applications? Again, it is a bit abstract on how your work can be \"helpful to researchers and practitioners alike by demonstrating shortcomings of recently proposed distributions over DAGs\". Can you try to assess on downstream causal‑discovery accuracy or decision‑making impact?\n- The simulation studies seem to be limited to $D\\leq 4$. Can you try a larger scale?\n- Hyper‑parameter tuning procedure is under‑specified.\n- Several parts of the manuscript are hindered by grammatical errors and unclear phrasing. Figures could be better annotated, and the paper would benefit from substantial proofreading. E.g., \" (Example 1 or capture dependences of edgesin the graphs\" Please fix them.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper positions itself as the first systematic study of expressiveness for several recently‑proposed parametric distributions over DAGs. The authors considered two practical gaps: (i) assign equal mass within a Markov‑equivalence class and (ii) model edge‑dependencies; and introduced (a) a recipe for importance‑sampling evaluation and (b) simple mixture models that overcome the gaps.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please find my comments and requested changes in Weaknesses above.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "QVDmzVJA5d", "reviewer": "Reviewer_4dVp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\nThis works effectively summarizes different causal structure learning approaches, discussing pros&cons of each algorithm and comparing them with a set of synthetic experiments. Moreover, an extension for DPM-DAG and ARCO-DAG is presented.\n\n## Weakness\n\nThe authors clearly explain why they decided to use only synthetic data, and the provided explanation is definitely reasonable. Nevertheless, the presence of experiments with real world data would have allowed the reader to understand the behavior of the different models in a non-synthetic scenario.\n\nFurthermore, the authors did not explicitly discuss the time complexity of the different models.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors of this work discuss the expressivity of different parametrized distributions over DAGs for causal discovery. Specifically they focus their attention on learning the structure of a Functional Causal Model under the assumption of causal sufficiency.\n\nThe distribution considered by the authors are:\n\n- Independent edges\n- Graph particles\n- DPM-DAG\n- ARCO-DAG\n- GFlowNet-DAG\n\nFurthermore, the authors propose to use mixture models of DPM-DAG and ARCO-DAG.\n\nFinally, the paper presents a set of synthetic experiments to compare per performances of the\nmodels in terms of statistical divergence between the learned model and the real one.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1.  I suggest dedicating an entire section to the mixture model approach in order to make more clear the separation between well known models and the proposed mixture model.\n2. I suggest to briefly discuss the execution time of the different approaches.\n3. Do you think that an extension of this work to the classical Bayesian network with discrete variables would be possible/useful as future work?", "broader_impact_concerns": "No Concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "FD9QWPVDQk", "reviewer": "Reviewer_G1Wi", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nS1. \nThe paper provides an instructive summary of recent Bayesian causal structure learning approaches. The review is well-written and easy to follow.\n\nS2. The submission identifies three illustrative problems that can provide insight into the fundamental limitations of causal structure learning algorithms and uses them to analyze the behavior of existing methods.\n\nS3. \nThe paper discovers that using mixtures of RPM-DAG and ARCO-DAG instead of individual models allows them to overcome their shortcomings in expressiveness.\n\n**Weaknesses**\n\nW1. The problems used to evaluate the algorithms are small in terms of the number of variables: two have three variables, and one has four. Furthermore, there is only one instance of each problem. How well the conclusions would scale to cases with more variables is unclear. \n\nW2. The paper's organization is slightly confusing. For instance, there is no clear definition of expressiveness or summary of the method limitations, and an explanation of the proposed mixture models is fragmented.\n\n*Kindly note that causal structure discovery is out of my area of expertise, so I cannot confidently assess the relevance or novelty of the results. I checked the general logic and presentation of the preprint.*", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper reviews Bayesian causal structure learning algorithms to uncover their limitations. It analyzes the expressiveness of distributions over directed acyclic graphs obtained in these methods and verifies the findings experimentally. Specifically, the preprint considers five model types. It investigates how they perform on three example problems, which aim to assess if algorithms can learn correctly in certain basic cases, such as when one considers Markov equivalent graphs or when one variable depends on two variables simultaneously. Based on the results, the preprint concludes that mixtures of RPM-DAG and ARCO-DAG models are the most expressive.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Crucial**\n\nC1. Currently, while testing for the fundamental issues, the examples are small, with either 3 or 4 variables and only one instance of each problem. \nBefore suggesting changes, I would appreciate clarification of the authors' perspective, specifically on whether it would be feasible to add experiments with several graphs and more variables in each example and whether additional experiments would add value to the results.\nBased on my understanding, additional experiments can help better assess methods for which evaluation results are not deterministic, such as GFlowNet-DAG, and also reason more confidently about how well the methods can scale.\n\n**Minor**\n\nM1. Add a definition of expressiveness somewhere at the beginning of the paper, along with an explanation of the approach used to measure it. Expressiveness often appears in varied contexts with different definitions, and works use different methods to measure it. Hence, a clear definition would help with the understanding.\n\nM2. Improve the consistency of the discussion of the suggested mixture models.\n\n (a) Add a subsection defining proposed mixtures and reference it in the Contributions paragraph since there is no dedicated section introducing them, and the explanation is fragmented. \n\n (b) Specify the models in the mixtures in the last sentence of the Abstract, \"To overcome them, we propose mixture models of distributions over DAGs,\" since otherwise, it suggests that an entirely new model is introduced.\n\n (c) Rephrase the last sentence in the conclusion, \"Therefore, we conjecture that a mixture of RPM-DAG models is particularly suited to scale Bayesian causal discovery to higher numbers of variables.\" which is not well-supported by the current experimental setup, given that the experiments are for a single instance of each problem and only for a small number of variables.\n\nM3. Add a summary statement or a table summarizing the method limitations since, currently, explanations are distributed across multiple sections.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "jTcI3h04kV", "reviewer": "Reviewer_sSDn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper considers several approaches for representation that already exist in the literature, and conducts a comprehensive evaluation on a single synthetic problem set. \n\nThe authors use clear explanations and notation in their exposition. The concepts and complications of probabilistic DAG representation are clearly explained, as are the evaluation methods.\n\nWhile mostly complete, there are some suggestions below which should be considered by the authors.\n\n## Section 2, paragraph Bayesian causal discovery\n\nThe notation $q(\\mathcal{D}|\\mathbf{G})p(\\mathbf{G})$ is introduced in the Bayesian causal discovery paragraph. There are two things to note: (1) This notation for a generative model should be introduced earlier in the Probability and random variables paragraph. (2) $q(\\mathcal{D}|\\mathbf{G})p(\\mathbf{G})$ is the product of a generative $q$ and the prior $p(\\mathbf{G})$, however in the subsequent paragraph on Bayesian model error line 5, it states \"The model for the prior distribution $q(\\mathbf{G})$...\". This is an ambiguous switch of notation without clarification.\n\n## Section 3.1 Independent edges\n\nAnother recent work that proposes the projection of a cyclic $\\mathbf{A}$ onto a non-cyclic set as the constraint function $h(\\mathbf{G})$ is Thompson et al (2024).\n\n## Section 3.3 Probabilistic models over orders\n\nThe permutation-upper-triangular representation has been discussed previously in the literature but these citations are missing. For instance, distributions over permutations (Stern 1990), and recent (preprint) work exploring distributions over this representation (Bonilla et al 2024, Davies et al 2025). In particular, an all-at-once fully-autoregressive representation over $\\mathbf{\\Pi},\\mathbf{U}$ was explored in Davies et al 2025 where an extension to the MADE (Germain et al, 2015) architecture was used to represent a fixed-size distribution over this tuple. \n\nEquation 8 in section 3.3 repeats much of equation 7, section 3.2. The split of the introduction and discussion permutation-upper-triangular representation over two sections should be revisited in order to reduce ambiguity. Especially since, as stated in line 10 of Section 3.2, the particle representation in general can be defined as a probabilistic representation.\n\n## Sections 4/5 Evaluation/Experiments\n\nRelying solely on the structural Hamming distance for evaluation of gthe concentration of probability mass (Table 3) can be misleading. For a more comprehensive evaluation, can the authors please also provide results using other standard scores such as F1 and/or AUROC? It might also be worth referencing or adding a discussion on why use of a single score such as SHD can be misleading.\n\nIt appears that the synthetic data used for evaluation is generated from a four-node DAG. Hence, the scalability of each probabilistic representation is not assessed in higher-node scenarios, a topic of recent interest in the research community. Perhaps this restriction could be highlighted to avoid ambiguity.\n\n## References:\n\nThompson, R., Bonilla, E. V., & Kohn, R. (2024). ProDAG: Projection-Induced Variational Inference for Directed Acyclic Graphs. arXiv preprint arXiv:2405.15167.\n\nHal Stern. Models for distributions on permutations. Journal of the American Statistical Association, 85(410):558–564, 1990.\n\nBonilla, E. V., Elinas, P., Zhao, H., Filippone, M., Kitsios, V., & O'Kane, T. (2024). Variational dag estimation via state augmentation with stochastic permutations. arXiv preprint arXiv:2402.02644.\n\nDavies, L., Mackinlay, D., Oliveira, R., & Sisson, S. A. (2025). Amortized variational transdimensional inference. arXiv preprint arXiv:2506.04749.\n\nGermain, M., Gregor, K., Murray, I., & Larochelle, H. (2015, June). Made: Masked autoencoder for distribution estimation. In International conference on machine learning (pp. 881-889). PMLR.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present a detailed review of probabilistic approaches for representation of directed acyclic graphs. In their experiments, they compare an assortment of recently contributed probabilistic representations and assess them via forward KL divergence and total variational distance. Due to the small enumeration of graphs in the evaluation (where the four node setup yields 543 possible DAGs) a complete evaluation over the entire enumeration is possible. The focus of the evaluation is to assess the expressiveness of each probabilistic representation considered.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "In addition to the notes in the Strengths and Weaknesses section above, the following changes should be considered:\n\nSection 2 Preliminaries, paragraph 1 probability and random variables. Line 6, the sentence \"We mainly apply this concise notation in the main text\" does not read well: \"mainly\" is superfluous.\n\nSection 2 Paragraph: Preliminaries, functional causal models. As this is a review, a figure would be well-placed to support the notation and concepts.\n\nSection 2 Bayesian causal discovery: the marginal likelihood arises by \"averaging over the (conditional) prior distribution\" is ambiguous. Changing \"averaging\" to \"taking the expectation over\" may help.", "broader_impact_concerns": "No ethical or broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)