ProDAG: Projected Variational Inference for Directed Acyclic Graphs
Authors: Ryan Thompson
Edwin V. Bonilla
Robert Kohn
Venue: NeurIPS 2025 poster ()
Keywords: structure learning
directed acyclic graph
variational inference
Primary area: probabilistic_methods

Abstract: Directed acyclic graph (DAG) learning is a central task in structure discovery and causal inference. Although the field has witnessed remarkable advances over the past few years, it remains statistically and computationally challenging to learn a single (point estimate) DAG from data, let alone provide uncertainty quantification. We address the difficult task of quantifying graph uncertainty by developing a Bayesian variational inference framework based on novel, provably valid distributions that have support directly on the space of sparse DAGs. These distributions, which we use to define our prior and variational posterior, are induced by a projection operation that maps an arbitrary continuous distribution onto the space of sparse weighted acyclic adjacency matrices. While this projection is combinatorial, it can be solved efficiently using recent continuous reformulations of acyclicity constraints. We empirically demonstrate that our method, ProDAG, can outperform state-of-the-art alternatives in both accuracy and uncertainty quantification.

ProDAG: Projected Variational Inference for Directed Acyclic Graphs
Ryan Thompson" University of Technology Sydney
Edwin V. Bonilla CSIRO's Data61
Robert Kohn University of New South Wales
Abstract
Directed acyclic graph (DAG) learning is a central task in structure discovery and causal inference. Although the field has witnessed remarkable advances over the past few years, it remains statistically and computationally challenging to learn a single (point estimate) DAG from data, let alone provide uncertainty quantification. We address the difficult task of quantifying graph uncertainty by developing a Bayesian variational inference framework based novel, provably valid distribu- tions that have support directly on the space of sparse DAGs. These distributions, which
to define prior and variational posterior, induced by pro- jection operation that maps an arbitrary continuous distribution onto the space sparse weighted acyclic adjacency matrices. While this projection is combinatorial, it can be solved efficiently using recent continuous reformulations of acyclicity constraints. We empirically demonstrate that our method, ProDAG, can outperform state-of-the tht-thh-rrt art alternatives in both accuracy and uncertainty quantification.
1 Introduction
directed acyclic graph (DAG) has directed edges between nodes and no directed cycles, allowing it to represent complex conditional independencies compactly. DAGs play central role in causal inference (Spirtes et al., 2001: Pearl, 2009) and consequently have found applications spanning broad range of domains such as psychology (Foster, 2010), economics (Imbens, 2020), and epidemiology (Tennant et al., 2021). The literature on DAGs, dating back at least to the 1980s (Lauritzen and Spiegelhalter, 1988; Pearl, 1988), has recently witnessed resurgence of interest due to a computational breakthrough by Zheng et al. (2018) that formulated a continuous characterization for the discrete notion of acyclicity. This breakthrough opened the door for scalable first-order algorithms that can learn DAG structures from datasets containing many more variables than previously assumed feasible. We refer the reader to Vowels et al. (2022) and Kitson et al. (2023) for surveys of this work.
Unlike frequentist approaches that only learn point estimate DAGs (Zheng et al., 2018: 8018; Bello et al., 2022; Andrews et al., 2023), Bayesian approaches learn entire posterior distributions over DAGs. The posterior quantifies graph uncertainty while acknowledging that the true DAG might be identifiable only up to its Markov equivalence class.? Uncertainty quantification also benefits downstream tasks such as treatment effect estimation (Geffner al., 2022) and experimental design (Annadani et al., 2023b). Most Bayesian approaches use Gibbs-type prior based on continuous acyclicity penalty (Annadani et al., 2021; Lorch al., 2021; Geffner et al., 2022) or model DAGs in an augmented space where acyclicity is readily satisfied (Charpentier et al., 2022; Annadani al., 2023a; Bonilla et al., 2024). While the former group (Gibbs-type priors) do not guarantee an exact posterior over DAGs, the latter group (augmented spaces) require relaxations or inference over discrete structures.
"Corresponding author. Email: ryan. thompson- edw.ee.uuuuuuuuuuuueee... au
'The Markov equivalence class the set of DAGs that encode the conditional independencies.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Wus
Wiz
Wize
W
prnty
0
D
W ~ P W (ro,(W) DAG F ball
Was
Figure I: Illustration of ProDAG's projected distributions. Samples (blue dots) from unconstrained continuous space (blue ellipse) are projected onto the nearest acyclic matrix within an & e-conntraieeed region (orange diamonds). Projected samples (orange dots) satisfy acyclicity and sparsity constraints. Theoretically, we show that for any continuous W IP the projection pro,, (W) is unique and measurable. This result implies that the projected distribution a valid distribution over DAGs.
We propose ProDAG, Bayesian method for learning DAGs that departs from the approaches above. Leveraging the power of variational inference, ProDAG models the prior and variational posterior using new distributions that natively have support the space of DAGs. These DAG distributions are induced by mapping a continuous distribution IP onto the space of DAGs via a minimal distance projection. Specifically, consider DAG represented by the weighted adjaceney matrix W R2*x with entry Wyk nonzero if and only if directed edge exists from node j to node k. We define distribution over such matrices W, and hence over DAGs, via the data generating process
W P, X A~m, W = pro, (W),
(1)
where À 2 0 is parameter governing the graph's sparsity (number of edges) with prior and
1 |IW - W
pro =
arg min
WEDAG, [Wlle SA
is the projection of the continuous matrix W RPxP onto the set of l, constrained acyclic matrices. Figure visually illustrates the process. We prove this projection is measurable mapping with unique solution, yielding valid probability distributions. A key quality of these distributions is that they place probability mass on exact zeros, necessary for sampling valid, sparsely connected DAGs.
Inspired by recently developed continuous DAG characterizations (c.g., Zheng et al., 2018), we reformulate the combinatorial optimization problem (2) as an equivalent continuous optimization problem that can be solved in parallel on a GPU. We use this formulation in variational framework by forming the prior and variational posterior using the new distributions, enabling scalable Bayesian DAG learning. Our framework accommodates DAGs represented as linear and nonlinear structural equation models (SEMs), making our approach suitable for problems and data of varying complexity. Experiments on synthetic and real data demonstrate that our approach often delivers more accurate inference than existing methods for DAG learning. Our toolkit for ProDAG is available on Gi tHub.
In summary, our core contributions are:
1. A new Bayesian model that directly samples from provably valid DAG distributions, ensuring exact acyclicity and eliminating the need for discrete approximations or relaxations;
2. A scalable variational inference framework that leverages GPU-accelerated projections, enabling efficient and accurate uncertainty quantification for linear and nonlinear DAGs;
3. An extensive suite empirical evaluations that demonstrate state-of-the-ar results across wide variety of synthetic datasets, and validated on a real biological dataset; and
4. A fser-friendly, well-documented, and open-source Julia implementation, made publicly available to promote adoption and reproducibility among researchers and practitioners.

---

2
Related work
Existing approaches for Bayesian inference on DAGs broadly fall into two groups: (1) regularization- based methods that include a continuous acyclicity penalty via Gibbs-type prior in the variational
objective (Annadani et al., 2021; Lorch et al., 2021; Geffner et al., 2022), and (2) auteeauugentttion methods that, e.g.. place prior over the graph's 's topological ordering (Viinikka et al., 2020; Cundy et al., 2021; Charpentier et al., 2022; Bonilla et al., 2024; Toth et al., 2024). The node-potential augmentation of Annadani al. (2023a) is related to the latter category. See also Wang et al. (2022) and Deleu et al. (2023) for approaches based on probabilistic circuits and generative flow networks. Our approach differs by designing prior and approximate posterior directly over DAGs without modeling the (discrete) topological ordering or using a penalty that does not place mass on zeros.
More broadly, our paper complements the literature on computational methods for point estimate DAG learning. The seminal paper of Zheng et al. (2018), which introduced the NOTEARS method, cleverly exploited the correspondence between the number of cycles in matrix and the trace of that matrix's power. A similar connection is now known to hold for other continuous functions and is exploited by methods such DAGMA (Bello et al., 2022), which we use in this work, as well those of Ng et al. (2020), Yu et al. (2021), and Gillot and Parviainen (2022). These approaches that use continuous characterizations of acyclicity are appealing since they scale to large graphs and can be augmented with discrete combinatorics (Manzour et al., 2021; Deng et al., 2023) if required. Though not pursued here, our projection onto DAG space also amenable to discrete techniques.
Finally, although distributions over DAGs via projections a new idea, others have used projections to construct distributions over non-grad non-graphical structures. Xu and Duan (2023) designed distributions over sparse regression coefficients by projecting a vector drawn from an arbitrary continuous distribution onto the & ball. They illustrated several computational advantages of their distributions priors in a Hamiltonian Monte Carlo (HMC) framework over classical spike-and-slab priors and prove that the resulting posterior enjoys the minimax concentration rate. Xu et al. (2023) extended this work to general proximal mappings with HMC, which includes the convex & projection as a special case. While our approach projects matrices onto the L ball, we also impose acyclicity, giving rise to nonconvex combinatorial optimization problem that does not constitute a proximal mapping.
3 Projected distributions
3.1 Deseription
To recap the introduction, we induce a distribution over weighted acyclic adjacency matrices W as W= = pro, (W), where pro, is defined in (2), W ~ P is a continuous matrix, and À ~ π bounds the & norm of W. Since acyelic graphs have self-loops (i.e., W always has zero diagonal), can take IP to be any continuous distribution over the off-diagonal elements of W (e.g., a multivariate Gaussian), with the diagonal elements fixed zero. The matrix W can thus be interpreted as the weighted adjacency matrix of a simple directed graph, with the projection removing all cycles.
Besides removing cycles, the projection maps W. onto the E ball of size A. This part of the projection is necessary for sampling parsimonious graphs. Without this constraint, the projection would return the densest possible acyclic graph with (p" ppppp// p)/2 edges. The advantage of constraining 4 ball, compared with modeling W as a Laplace distribution, is that the l ball gives probability mass exact zeros beyond that already given by the acyclicity constraint. It is known that Laplace distributions do not concentrate exactly on zero except at their mode (Park and Casella, 2008).
to the
On a conceptual level, the resulting distribution p(W) can be understood components and exact zeros.* 4 To elaborate, first observe that the nonzeros in W may be indexed by a sparsity pattern S E {0, 1 1 (i.e., an unweighted acyclic adjacency matrix). The nonzeros are then Ws. where Ws is the restriction of W to S. Because the & constraint shrinks the nonzeros towards zero, the continuous Ws is transformation fx of Ws, i..., Ws fxWs). where fx is implicit in
mix of continuous
'The distribution can be Dirac for distribution over graphs of fixed sparsity (same number of edges) an exponential for a distribution over graphs of varying sparsity. "This mixture means that the projected distribution is multimodal, as different regions of the continuous distribution map distinct acyclic graphs (except in degenerate point-mass cases).
3

---

the projection. This mixture interpretation suggests an interesting representation of p pWW):
p(W) Σ P(W,, S.)p(S)) - Σ P(fA((Ws.) (.),(S)),
where St...., Sm is the finite collection of all acyclic adjacency matrices. We emphasize that this representation is purely conceptual and therefore we do not define or use the corresponding distributions over S. Ws. or Ws. However, it highlights that p(W) forms a distribution over acyelic adjacency matrices S without explicitly having model the discrete, high-dimensional probability mass function pS)). Moreover, as we later show, the projection can be reformulated as a continuous optimization problem, meaning one sample from p(W) without combinatorial complexity.
3.2 Properties
One might ask if p(W) a valid distribution. Theorem 1 addresses this question in the affirmative by showing that the projection is almost surely unique and measurable. The proof is in Appendix A.
Theorem 1. Let W be endowed with a continuous probability measure. Then it holds:
1. Projection (2) is unique almost surely with respect to the measure on W; and
2. Projection (2) is measurable with respect to the measure on W.
Theorem 1 states that, for any continuous W ~ P. the function pro, (W) is uniquely defined and measurable. Uniqueness here means that for almost every W, the projection has a single minimizer; different values of W may still map to the same W, but no single W admits multiple optimal projections. This important result implies that projected distribution is valid distribution since measurable function y applied to a random variable = produces well-defined random variable g(2).
Scalable projections
4.1 Continuous reformulation
To facilitate scalability, we reformulate the combinatorial projection (2) as an equivalent continuous optimization problem. This recharacterization follows recent work (Zheng et al., 2018; Bello et al., 2022) that established a useful and important equivalence relationship for certain functions h(W)::
k(W)= = 0 + W E DAG.
Bello et al. (2022) showed that for all matrices W E W = (W E RPXP: p(W o W) < 1}. where the spectral radius and the Hadamard product, the function h(W) log det([ W W) satisfies this property. This function, which constitutes the main innovation of DAGMA, attains the value zero if and only there are no cycles in W. Thus, its level set zero exactly corresponds the entire set of DAGs. Though it is a nonconvex function, ((W) is useful because it is both continuous and differentiable in W. We thus exploit this equivalence of constraints and rewrite projection (2) as
pro (w) WEW, - [Wile, SA
arg min
| -
As diseuss next, it is possible (i) solve this projection using first-order methods and (ii) evaluate its gradients analytically. These attributes pave the way for scalable, gradient-base posterior learning.
4.2 Algorithms
To solve projection (3), we employ two-step procedure that produces a solution the intersection of sets. Step 1 projects the matrix W onto the set of acyclic matrices. We call the result W. Step projects W. onto the & ball, yielding the final output W that remains acyelic but is further sparse.
For the acyclicity projection in Step 1, we use a fath-following algorithm similar Bello et al. (2022).
This algorithm reformulates the problem by shifting the acyclicity constraint into the objective:
min f (W; W) I - w11 + h(W). Wew
(4)

---

Here, 2 0 is taken along a sequence μ()) > v1(2) >... > p.) The wolh-following algorithm
solves (4) at μ = μ ((+1)) by using the solution at ( as an initialization point (see Appendix B). The limiting solution as - is guaranteed to satisfy the acyclicity constraint hu (WW = 0 (Lemma 6, Bello et al., 2022). Crucially, the optimization problem (4) is amenable to gradient descent. The gradient of f readily follows from the gradient of h. given by Vh(W) 2(1 W w)- W.
For the 2 projection in Step 2. we employ a matrix version of the & projection algorithm developed by Duchi et al. (2008). This non non-iterative algorithm soft-thresholds the elements of the weighted adjaceney matrix W that is produced by Step 1. This second step results in the final matrix W, solution to the original projection (3). For completeness, we provide the algorithm in Appendix B.
4.3 Gradients
In many Bayesian learning frameworks, including our variational framework, it is necessary to compute gradients. In our case, this computation is complicated by the presence of the projection since we need its gradients with respect to W and X. Though it is technically possible to differentiate through our algorithms using automatic differentiation, it is practically infeasible due to the iterative nature of the path-following algorithm. As it turns out, it is possible to analytically evaluate the gradients via the implicit function theorem. This technique involves differentiating the Karush- Kuhn-Tucker (KK')) conditions of the projection that implicitly define W function of W and A. Intuitively, the KKT conditions express how optimality of W is preserved under small perturbations of W or X, and the implicit function theorem formalizes this dependence to yield local expressions for how W varies. The gradients from this technique are in Appendix C alongside their derivation.
5 Posterior learning
5.1 Prior and variational posterior
Suppose we have data X E Raxp with n observations on p variables I E RR. Let p(W) be a projected prior distribution over the DAG W. The conditional likelihood p(X W) is given by the linear SEM x = W x + E, where E ~ E). Our task is to infer the posterior distribution P(W X). which is analytically intractable. Though our distribution admits convenient algorithms, remains taxing to learn the posterior via traditional approaches due to the high-dimensionality of the posterior space. We thus approximate the posterior using variational inference, leading to ProDAG.
Variational inference (e.g.. Blei et al., 2017) optimizes the parameters 0 of an approximate posterior qe (W) such that it is nearest to the true posterior p(W) X) in Kullback-Leibler (KL) divergence. Like the prior, the approximate posterior is a projected distribution. However, as we discuss later, performing variational inference W alone is hard because the projection is non-invertible (multiple W project the same W). Hence, we focus on learning the joint posterior over W and W. The joint posterior does not suffer from non-invertibility because for any W the associated W is also observed.
Let ( W X) be the true joint posterior and let go W) be the variational joint posterior. Since W is deterministic function of W. the variational joint posterior is the product of the variational marginal posterior qo (W) and Dirac delta function: qa (W, W) go (W))(W pro (W)). The marginal qo (W) can be continuous distribution such as multivariate Gaussian. In that case, the parameters 0 contain means and covariances. Choosing qo (W) is the same as choosing F in (1). For simplicity and clarity of exposition, we treat the sparsity parameter À as fixed, i.e., π in (1) is a Dirac.
5.2 Evidence lower bound
As in standard variational inference, we proceed by maximizing the evidence lower bound (ELBO):
EL) F. g. Plog((X W, W)] [T (W, W) p(W, W)).
(5)
The first term on the right- side of (5) is the expectation of the log-likelihood under the variational posterior. Since the log- likelioodoood depends only on W through W, this term simplifies as
yn (W. w) hiw.ww/logg [logp(X W, W)] = (log W)).
5

---

To estimate this expectation, we sample DAGs from go (W) and evaluate the log- likelihhodd on these samples. To sample from go (W), we need only sample from go (W, W) and discard W. The second term on the right-hand side of (5) is the KL divergence between the variational posterior qo (W, W) and the prior p(W, W) p(W) Lro (W)). Since W is a deterministic function of W, it is not too difficult to see that the KI. divergence between the joint posterior and prior is equal to the KL divergence between the marginal posterior and prior on W. Hence, the second term also simplifies as KL[go(W, W) II p(W, W)] = KL[ (W) |I p(W)]
Provided that (W) and p p(W) are chosen from distributional families with closed-form density functions such as multivariate Gaussian distributions, this term is amenable to analytical evaluation. The gradients of Gaussians and other (more complex) distributions for W are readily accommodated in the ELBO using the well-known reparameterization trick (see, e.g., Kingma and Welling, 2014).
The challenge mentioned earlier of learning the marginal posterior of W arises due to the KL divergence appearing in (5). If we instead focus on the marginal posterior, the divergence term in (5) is between go W)) and p(W). With some work, one can show that the divergence between these terms involves the divergence between the (unknown) prior and posterior conditional distributions of W given W (since W is non-invertible in W). In any case, this issue is avoided in the joint posterior.
5.3 Algorithm
In light of the above remarks, we maximize the ELBO using the gradient descent-based Algorithm 1.
Algorithm 1 ProDAG
Input: Initialization 0, data X E R" learning rate 1) >> O, no. of DAG samples L E N while Not converged do
Sample W(() ~ (W) and set W() prox (W) for I....,L
Compute ELBO = El logp(X X (W)) - KL (W) || p(W)
Update 0 < 0 + vVV (OF
end while
Output: Optimized parameters 0
The most taxing part of Algorithm is the step that samples the W matrices for the expected log- likelihood. Specifically, at each iteration, the algorithm projects a sample of L matrices W onto the constrained space. This projection, formulated with the continuous acyclicity constraint, involves the inversion of p p matrices, which has cubic complexity. Fortunately, these can be done in parallel on a GPU using a batched CUDA implementation of the projection algorithms in Appendix B. Meanwhile, the gradients (required in the update step for the variational parameters ()) have quadratic complexity in p (see Appendix C). Hence, the overall complexity of training ProDAG is O(p²).
Though cubic complexity is not ideal, most state-of-the-ar art frequentist and Bayesian approaches to DAG learning also have cubic complexity (e.g., Lorch et al., 2021; Bello et al., 2022; Annadani et al., 2023a). Moreover, cubic complexity is the worst case complexity. We investigate the expected case by plotting the average training time over a range of p in Figure 2. This figure also reports the average time perform a parallel projection of 100 matrices, since this is the main source of complexity. Though certainly not scaling linearly, the times are reasonable for many applications, with training taking under minute for p 20 variables and slightly more than hour for p 100 variables.
6
Nonlinear DAGs
6.1 Structural equation model
The discussion so far assumes the SEM is linear. However, our approach accommodates nonlinear SEMs of the form x = f(x) + where ʃ IRP - RR is nonlinear function with an acyelic Jacobian matrix. Acyelicity of the Jacobian matrix is necessary for a DAG because 0fa (x)) /Ox; * 0 indicates an edge exists from node j to node k. standard choice of the function I is a neural network structured as f(x) (fi fo(r)) T. where fx( is a feedforward sub-network for node k.

---

Training time, (sccords))
Projection time (secardii)
3020
Number variaalls
Number of variables
Figure 2: Computation times in seconds with a sample size n 100. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
fp- Let w = (wi,...,,wp wip) be a p-tuple of neural network weights corresponding to f1..... Denote
by wht E Rxxd the restriction of the kth network's weights wi to the first hidden layer, where is the number of neurons in that layer. If the jth row of wf. denoted w/x E R is a zero vector (i.e., |k /|2 = 0), then f. (z)/0x; 0 and no edge exists from node to node k. Hence, we can force the Jacobian of (1) to be acyclic by constraining the p x p matrix ((w) [Ila)[ be acyclic.
6.2 Posterior learning
Towards nonlinear DAG, we construct distribution over all weights w of the neural network f(z): IP, nr, p
where P is continuous distribution over the network weights. This data generating process yields a fully Bayesian neural network with every layer having prior over its weights. The projection returns the weights w nearest to ß such that the resulting adjacency matrix W (w) is acyclic and sparse:
P
pro (
arg min W(W) (W(()))= ([Wlll SX
Σ 1 wlll
(6)
The matrix W (w) forming the constraints in (6) is a function only of the first hidden layer weights. Hence, the projection only modifies those weights and leaves the remaining layer weights untouched.
Though the number of weights on the input layer is more numerous than the number of weights in the linear SEM case by a factor of d, the projection reduces to simpler problem that is no more complex than the linear case. Proposition presents this result. The proof is Appendix Proposition 1. Let is be a tuple of neural network weights and let a/y be the vector of weights for the jth input to the first hidden layer of the kth output. Define the matrix W elementwise as Wyk || ll2 and let W be the projection of W in (3). Then projection (6) is solved by taking
Wyk
w W
Wik
for j. 1, .......p and leaving the remaining layers' weights unmodified.
Proposition states that we need only project the matrix W constructed from w to obtain the weights for acyclic neural network fx)) Consequently, the nonlinear
requires new algorithms.
Performing posterior inference on the neural network weights w within our variational framework requires no major modifications. The only change required is to re-express the ELBO in terms of w:
EL) Eg [log p(X w)] KL (w) || p())),
where the variational parameters 0 now parameterize a distribution of the same dimension as w.
7 Experiments
7.1 Baselines and metries
Our method ProDAG is benchmarked against several state-of-the-art methods: DAGMA (Bello et al., 2022), DiBS and DiBS+ (Lorch et al., 2021), BayesDAG (Annadani al., 2023a), and BOSS (Andrews
7

---

et al., 2023). DAGMA is a frequentist version of our approach that delivers a point estimate DAG. DiBS, DiBS+, and BayesDAG are alternative Bayesian approaches that learn posterior distributions over DAGs. BOSS is frequentist comparator that relies on discrete (greedy search) methods rather than continuous acyclicity characterization. Appendix E. describes hyperparameters and implementations.
We consider several metrics that characterize different qualities of the learned posterior distributions. As measure of uncertainty quantification, we report the Brier score under the posterior. As measures of structure recovery, we report the expected structural Hamming distance (SHD) and expected FI score under the posterior. Finally, as measure of discriminative power, we report the AUROC under the posterior. Lower values of Brier score and expected SHD are better, while higher values of expected F1 score and AUROC are better. All four metries compare the learned posterior to the ground truth DAG. Since DAGMA and BOSS return single point estimate graph, the metrics are evaluated by taking its posterior as Dirac delta distribution that places all mass on the learned graph. For all other methods, the metrics are evaluated on a sample of graphs from the learned posterior.
7.2 Linear synthetic data
Figures 3 and 4 report results on datasets simulated from linear Erdós Erdös-Rényi DAGs with p = 20 and p 100 nodes (s 40 and 200 edges). The sample size varies from 10 to 1000.
PrODAC DRS RayesDAG
DAGMA DISS DDSS
100
1.00
75
ExpectedF1 score
0.75
Bieer score
Expected SHD
AUROC
0.30
120
124
1244
104
109
124
70%
70m4
10
12'
104
1024
104
102
7244
1020
10m0
124
Sample
Sample
Sample
Sample
Figure 3: Performance on synthetic datasets generated from linear Erdós Erdös-Rényi DAGs with p = 20 nodes, s = 40 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are mcasured over 10 independently and identically generated datasets.
Rényi
ProDAS
BayesDAG
DAGMA DuBS+ ROKS
1060
1.00
830
El score
0.75
brief score
620
Exprect ed
AUROC
88
0.50
420
220
0.25
0.6
437
105
1077
101
107
107
1000
1000
100
12'4
1007
1000
5377
100
1077
1377
1000
1077
120
Sample Size
Sample
Sample
Sample
Figure 4: Performance on synthetic datasets generated from linear Erdös-Rényi nodes, 8 = 200 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
Rényi DAGs with 100
Appendix F details the full setup. ProDAG delivers comparatively strong performance across the full suite sample sizes, especially regarding uncertainty quantification. For small sample sizes, where there is little hope of recovering the true graph, ProDAG is as competitive other Bayesian methods. Meanwhile, for moderate and large sample sizes, ProDAG overtakes the Bayesian methods across most metrics. ProDAG performs approximately well as the frequentist DAGMA and BOSS for
8

---

the largest sample size, where little uncertainty remains regarding the true graph. Except for their discriminative power (AUROC), DiBS and BayesDAG do not significantly improve with the sample size. Lorch et al. (2021) studied DiBS on sample size of 100, and their results roughly match ours at n = 100. Annadani al. (2023a) evaluated BayesDAG on smaller graphs of size p =
7.3 Nonlinear synthetic data
Figures 5 and 6 report results on datasets simulated from nonlinear DAGs. Following the setup in
DACMAA EayssDAG
PIDDA DES
Expected Ft score
Biner score
28
M
AURDC
Expected
1-|
10
0.25
Sample sire
Sample Sizp
Sample sire
Sample size
Figure Performance on synthetic datasets generated from nonlinear Erdós Rényi DAGs with p 10 nodes, 20 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
DegesDAG
Expected F1 score
Bieer Score
50
AUROC
Expected
06
120
1244
1.50
1053
70
002
124
105
7227
100
13"
100
1054
100
105
12.1
102
1997
125
Sample
Sample size
Sample
Sample size
Figure 6: Performance synthetic datasets generated from nonlinear Erdōs-Rényi DAGs with p = 20 nodes, s = 40 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
Zheng et al. (2020), each variable in the graph is generated from a single hidden-layer neural network whose inputs are that variable's parents. The learning task is more challenging than in the linear setting (more than 10x parameters in each graph), so we study graphs with 10 and p 20 nodes (s 20 and 40 edges). BOSS cannot model nonlinear DAGs, so is excluded. DiBS and DiBS+ exceed available memory when p 20, so are excluded in that setting. ProDAG delivers strong performance and improves over the frequentist and Bayesian baselines most of the time. Again, particularly excels at uncertainty quantification. The results for DiBS approximately match those Lorch et al. (2021), where its expected SHD is similar to that of the null graph. Annadani et al. (2023a) studied BayesDAG here, ProDAG does not require a sample of similarly large size to produce high-quality inference.
sample size of 5000, much bigger than the largest sample size
7.4 Real data
The flow cytometry data of Sachs et al. (2005) is biological dataset designed to aid the discovery of protein signaling networks. The dataset contains n = 7466 human cell measurements on p = 11 phosphoproteins and phospholipids. To our knowledge, it is the only real dataset with expert

---

consensus graph available, so it is widely used for real-world evaluation of DAG learning approaches. Table reports results on the flow cytometry dataset. with the consensus graph consisting of 18 directed edges taken as the ground truth. ProDAG outperforms all approaches across all metrics. It.
Table 1:
Performance on the flow cytometry dataset of Sachs et al. (2005). The averages and standard emros are measured over 10 splits of the data. The best value of each metric is indicated in bold.
Brier score
Exp. SHD
Exp. F1 score (%)
AUROC (%)
ProDAG
16 0.2
0.2
24 17 + 1.4
0.6
1 1.6
DAGMA
36 + 2.1
29 1 1.6
49 1.3
DiBS
16 - 0.2
22 0.1
16 16
55 = 1.6
DiBS+
26 = 1.1
2
15 +
51 =1=1 1.4
BayesDAG
25 2,5
59 1.6
BOSS
34 0.7
28 + 0.5
15 + 1.3
19 0.8
achieves the lowest Brier score, tied with DiBS, and surpasses DiBS+ and BayesDAG in uncertainty quantification. Additionally, ProDAG demonstrates superior performance in expected SHD, expected Fl score, and AUROC, reinforcing its strong overall capability in recovering the consensus graph structure. Meanwhile, ProDAG uniformly improves on DAGMA, highlighting the Bayesian benefits of ProDAG over DAGMA, since both are identical in their continuous characterization of acyclicity,
7.5 Additional experiments
suite of additional experiments are included in Appendix G. In partieular, Appendix G.1 presents results from further synthetic experiments covering (i) denser graphs, (ii) scale-free graphs, (iii) non-Gaussian noise distributions, and (iv) heteroscedastic noise distributions. Appendix G.2 reports results on semi-synthetic datasets generated from a large real-world medical diagnostic network. Appendix G.3 provides comparisons with a Markov chain Monte Carlo (MCMC) baseline. Finally, Appendix G.4 compares the run times of ProDAG with those of its primary Bayesian competitors.
Concluding remarks
DAGs are complex structures due to their acyclic nature, making Bayesian inference formidable task. ProDAG is a novel approach that employs new projected distributions that have support directly on the space of DAGs. We study key properties of these distributions and show that they facilitate state-of-the-arr Bayesian inference for linear and nonlinear DAGs in scalable variational framework. While ProDAG is computationally costlier than frequentist methods, is well-suited to settings where reliable posterior uncertainty critical (e.g., decision making), where limited data make sparsity valuable (e.g.. biology), or where strict acyclicity and interpretability are essential (e.g., policy).
As with other approaches to structure discovery using purely observational data, identifiability of the DAG structure in ProDAG depends on the correspondence between the assumed likelihood and the true data generating process. When the likelihood correctly specifies the noise distribution implied by the underlying SEM, standard identifiability results apply under assumptions such as causal sufficiency, the Causal Markov condition, and faithfulness. For instance, in a linear SEM with homoscedastic Gaussian noise, the DAG becomes fully identifiable (Peters and Bühlmann, 2014). In this setting, which experiments cover, ProDAG's posterior increasingly concentrates correct DAG as the sample size grows. Conversely, the likelihood is misspecified, different graphs may induce indistinguishable likelihoods, leading to loss of identifiability. In such cases, ProDAG represents the remaining structural ambiguity through the spread of its posterior distribution.
the
Recent work has extended continuous acyclicity approaches to interventional data (Brouillard et al., 2020; Faria et al., 2022; Lopez et al., 2022). Though our focus is observational data, ProDAG could be adapted for interventional settings by replacing the observational likelihood with one appropriate for interventions while leaving the projected distributions unchanged. Developing this extension and assessing ProDAG in interventional contexts offers natural and important next step for research.
10

---

References
Bryan Andrews, Joseph Ramsey, Ruben Sanchez-Romero, Jazmin Camehong, and Erich Kummerfeld. Fast scalable and accurate discovery of DAGs using the best order score search and grow-shrink trees. In Advances in Neural Information Processing Systems, volume 36, pages 63945-63956 2023.
Yashas Annadani. Jonas Rothfuss, Alexandre Lacoste, Nino Scherrer, Anirudh Goyal, Yoshua Bengio, and Stefan Bauer. Variational causal networks: Approximate Bayesian inference over causal structures. 2021. arXiv: 2106.07635.
Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong. BayesDAG: Gradient- posterior inference for causal discovery. In Advances in Neural Information Processing Systems, volume 36, pages 1738-1763, 2023a.
Yashas Annadani, Panagiotis Tigas, Desi R Ivanova, Andrew Jesson, Yarin Gal, Adam Foster, and Stefan Bauer. Differentiable multi-target causal Bayesian experimental design. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 34263-34279, 2023b.
Kevin Bello, Bryon Aragam, and Pradcep Ravikumar. DAGMA: Learning DAGs via M-matrices and a degedettrminat acyclicity characterization. In Advances in Neural Information Processing Systems, volume 35, pages 8226-8239, 2022.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the Americun Statistical Association, 112:859-8777 2017.
Edwin V Bonilla, Pantelis Elinas, He Zhao, Maurizio Filippone, Vassili Kitsios, and Terry O' O'Kane. Variational DAG estimation via state augmentation with stochastic permutations. 2024. arXiv: 2402.02644.
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexan- dre Drouin. Differentiable causal discovery from interventional data. In Advances in Neural Information Processing Systems, volume 33, pages 21865-21877, 2020.
Bertrand Charpentier, Simon Kibler, and Stephan Günnemann. Differentiable DAG sampling. In International Conference Learning Representations, 2022.
Chris Cundy, Aditya Grover, and Stefano Ermon. BCD Nets: Scalable variational approaches for Bayesian causal discovery. In Advances in Neural Information Processing Systems, volume 34, pages 7095-7110, 2021.
Tristan Deleu, Mizu Nishikawa-Toomey, Jithendaraa Subramanian, Nikolay Malkin, Laurent Charlin, and Yoshua Bengio. Joint Bayesian inference of graphical structure and parameters with a single generative flow network. In Advances in Neural Information Processing Systems, volume 36, pages 31204-31231, 2023.
Chang Deng, Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. Optimizing NOTEARS objectives via topological swaps. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 7563-7595, 2023.
John Duchi, Shai Shalev-Shwartz, Yoram Singer, and Tushar Chandra. Efficient projections onto the balball for learning in high dimensions. In Proceedings of the 25th International Conference on Machine Learning, pages 272-279, 2008
Gonçalo R A Faria, André F Martinsf, and Mário A T Figueiredo. Differentiable causal discovery under latent interventions. In Proceedings of the 1st Conference on Causal Learning and Reasoning, volume 177, pages 253-274, 2022.
E Michael Foster. Causal inference and developmental psychology. Developmental Psychology, 46: 1454-1480, 2010.
Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Agrin Hilmkil, Joel Jennings, Nick Pawlowski, Miltiadis Allamanis, and Cheng Zhang. Deep end-to-end causal inference. NeurIPS 2022 Workshop on Causal Machine Learning for Real-World Impact, 2022.

---

Pierre Gillot and Pekka Parviainen. Learning large DAGs by combining continuous optimization and feedback are set heuristics. In Proceedings of the AAAl Conference on Artificial Intelligence. volume 36, pages 6713-6720, 2022.
Guido W Imbens. Potential outcome and directed acyclic graph approaches to causality: Relevance for empirical practice in economics. Journal of Economic Literature, 58:1 119-1179, 2020.
Michael Innes, Elliot Saba, Keno Fischer, Dhairya Gandhi, Marco Concetto Rudilosso, Neethu Mariya Joy, Tejan Karmali, Avik Pal, and Viral B Shah. Fashionable modelling with Flux. In Workshop on Systems for ML and Open Source Software at NeurIPS 2018, 2018.
Diederik P Kingma and Jimmy Lei Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
Diederik P Kingma and Max Welling. Auto-encoding variational Bayes. In International Conference Learning Representations, 2014.
Neville Kenneth Kitson, Anthony C Constantinou, Zhigao Guo, Yang Liu, and Kiattikun Chobtham. A survey of Bayesian network structure learning. Artificial Intelligence Review, 56:8721-8814, 2023.
S L Lauritzen and D J Spiegelhalter. Local computations with probabilities graphical struc- tures and their application to expert systems. Journal of the Royal Statistical Society: Series (Methodological), 50:157-224, 1988.
Romain Lopez, Chrrsthitnaann Hütter, Jonathan K Pritchard, and Aviv Regev. Large-scale differen- tiable causal discovery of factor graphs. In Advances Neural Information Processing Systems, volume 35, pages 19290-19303, 2022.
Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. DiBS: structure learning. In Advances in Neural Information Processing Systems, volume 34, pages 24111-24123, 2021.
Differentiable Bayesian
Hasan Manzour, Simge Küçükkavuzz Hao-Hsiang Wu, and Ali Shojaie. Integer programming for learning directed acyelic graphs from continuous data. INFORMS Journal on Optimization, 46-73, 2021.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. In Advances in Neural Information Processing Systems, volume 33, pages 17954,17954, 2020.
Trevor Park and George Casella. The Bayesian lasso. Journal of the American Statistical Association, 103:681-686, 2008.
Judea Pearl. Probabilistic Reasoning in Intelligent Systems. Morgan Kaufmann, San Francisco, USA, 1988.
Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, New York, USA, 2nd edition, 2009.
J Peters and P Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101:219-2281 12021992222222 2014.
Karen Sachs, Omar Perez, Dana Pr'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308:523-529, 2005.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. MIT Press, New York, NY, USA, 2001.
Peter W G Tennant, Eleanor Murray, Kellyn F Amold, Laurie Berrie, Matthew P Fox, Sarah C Gadd, Wendy J Harrison, Claire Keeble, Lynsie Ranker, Johannes Textor, Georgia Tomova, Mark S Gilthorpe, and George T H Ellison. Use of directed acyelic graphs (DAGs) to identify confounders in applied health research: Review and recommendations. International Journal of Epidemiology, 50:620-632 2021.
12

---

Ryan Thompson, Edwin V Bonilla, and Robert Kohn. Contextual directed acyclic graphs. In Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, volume 238, pages 2872-2880, 2024.
Christian Toth, Christian Knoll, Franz Pernkopf, and Robert Peharz. Effective Bayesian causal inference via structural marginalisation and autoregressive orders. In Structured Probabilistic Inference & Generative Modeling Workshop of ICML 2024, 2024.
Jussi Viinikka, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. Towards scalable Bayesian learning of causal DAGs. In Advances in Neural Information Processing Systems, volume 33, pages 6584-6594, 2020.
Matthew Vowels, Necati Cihan Camgoz, and Richard Bowden. D'ya like DAGs? structure learning and causal discovery. ACM Computing Surveys, 55:1-36, 2022.
survey on
Benjie Wang, Matthew Wicker, and Marta Kwiatkowska. Tractable uncertainty for structure learning. In Proceedings of the 39th International Conference on Machine Learning, volume 162, pages 23131-23150, 2022.
Maoran Xu and Leo L Duan. Bayesian inference with the 11 -all prior: Solving combinatorial prob- lems with exact zeros. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 85:1538-1560 2023.
Maoran Xu, Hua Zhou, Yujie Hu, and Leo L Duan. Bayesian inference using the proximal mapping: Uncertainty quantification under varying dimensionality. Journal of the American Statistical Association, 2023.
Yue Yu, Tian Gao, Naiyu Yin, and Qiang DAGs with no curl: An efficient DAG structure learning approach. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 12156-12166, 2021.
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Sociery: Series B (Statistical Methodology), 68:49-65 67, 2006.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. DAGs with NO TEARS: Continu- ous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. Learning sparse nonparametric DAGs. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statisties, volume 108, pages 3414-3425, 2020.
13

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction accurately reflect the remainder of the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed in Appendix H.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Assumptions where relevant are stated. Proofs are provided in the appendices.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Core details for reproducing the experimental results are included in Appen- dices and F. Code for exactly reproducing the experimental results is released.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code, data, and instructions for reproducing the experiments is available.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Experimental settings are detailed in Appendix E.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: One standard errors of the mean are reported for all experiments.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Computational resources described in Appendix E.
9. Code of ethics
14

---

Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https:/ http:://werrrp...pppoooooahcaaaa hctps//ecpeeeeeiiiccccceeodaaaa.
Question:
Answer: [Yes]
Justification: The paper conforms with the NeurIPS Code of Ethics.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Societal impacts are discussed in Appendix I.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses such risks.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All existing assets are cited in the paper. Version numbers and URLs are provided. The license name is included for each asset. See Appendix E.
13. New assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [Yes]
Justification: package ProDAG is released with accompanying documentation.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: The paper does not involve research with human subjects.
16. Declaration of LLM usage
Question:
Does the paper describe the usage of LLMS if it is important, original, or non-standard component of the core methods this research? Note that if the L.LM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development does not involve L.MSs.
15

---

A Proof of Theorem 1
Proof. We proceed by proving each claim of the theorem in turn.
Uniqueness We begin by proving the first claim that projection (2) is unique almost surely. We
assume hereafter that À 0 (if X = 0 uniqueness holds trivially). First, observe that the weighted
adjacency matrix W can be expressed in terms of an unweighted (binary) adjacency matrix S e {0, 1 p and a continuous matrix V RPx::
min WEDAG, [Wllle
|
WI?
min SE SEDAG VER Sx
fs(V; W) = 1|W - S O
where @ denotes the Hadamard product and we adopt the convention that Vyk: = 0 Syt = 0. We can write the minimization on the hinht-hand side as min-min problem:
min SE(0.1) SEDAG VERPER, IVVll
fs(V; W)
min SE(0.1)xx> VERP*y SEDAG
min Js(V:W)
[[VIl SA
Consider the inner optimization problem on the right-hand side. For any W, let (W) denote optimal solution this problem, and define the set of W with non-unique solutions as
As == (W E R"*P| fV (W); W) = fV (W); W)) VT(WW * Vi (W)).
Since (V; W) is strictly convex in V and the & ball is convex, the optimal solution V*(W) is unique for all W. Hence, the set As is empty. Now, define the set of W with non non-unique solutions across any two different adjacency matrices S1 and S2
Ss = (W E RRXx/ Js, (VW
Ss, (V; (W); W), VT(WW) * VV))
Because Vit Vit) + (W), there at least one coordinate (j, k) on which they differ. Equality of the two quadratic objectives then forces W to lie on a single affine hyperplane in RPxp, Such hyperplane has measure zero when W is endowed with a continuous probability measure. Moreover, because there are finitely many pairs S2) with S A S2, the union of As over all pairs also has measure zero. We thus conclude that projection (2) is unique almost surely with respect W.
Measurability
To prove the second claim that projection (2) is measurable, define the functions
gA (W, .) == arg min
W SS VI 2
V VERPXP V 1A, SA
where again adopt the convention that Vyk 0 if Syk
0, and
hx hw
arg min se(0.1 SEDAG
W 5 (W, S)
where, if the minimizer not unique, we choose the lexicographically smallest It follows that
pro, (W) = (W (W)).
We now demonstrate that each of these functions are measurable in W. First, for fixed S, observe that the function 9x (W, S) is the minimizing argument of strietly convex optimization problem (with compact feasible set) and W + 1/2[[W S O VII is continuous in W. Hence, by Berge's maximum theorem, 9x (W. :S) continuous in W. Continuous functions are measurable with respect to their argument, and therefore ga (W, S) is measurable W. Second, since the set of binary adjacency matrices is finite, the function h2 (W) takes the pointwise minimum finitely many measurable functions and returns the adjacency matrix indexed by that minimum (breaking any ties per the lexicographic order). A pointwise minimum of measurable functions is itself measurable. Lastly, since gx ( and hx(W) are both measurable in W, their composition gx gx (W, h (W) must also be measurable in W. We thus conclude that projection (2) is measurable with respect to W.
16

---

B Projection algorithms
Algorithm 2 provides the routine for the projection of the matrix W onto the set of acyelic matrices. As discussed in the main text, the algorithm minimizes penalized objective function along
Algorithm 2 Solver for acyclicity projection
Input: Continuous W TRPx,, initialization W()) E W, and scaling coefficients {u(0) for t = 0, ,....,,'' , T do
Initialize W at (())) and set
((t-1) arg min fa W) Wew
end for
Output: Acyclic matrix W- = W(T)
decreasing sequence of scaling coefficients 1( > > (T At each stage, the solution W(()) serves warm start for the next problem at (((+1), This continuation scheme produces sequence of unconstrained matrices [W")) that converges μ 0 to point satisfying the constraint h(W) 0. In our implementation,
set (()) 1, p(++1) (t) /2, T 10, and (0 0 W.
Algorithm 3 projects the acyclic matrix W from Algorithm 2 onto the & ball of size X. The algorithm
Algorithm 3 Solver for & projection
Input: Acyclic matrix W E RPxP and sparsity parameter À 2 0
Sort vec(W) V1l > V22 1.... 2..22 Mpxpl
Set
Jmax + max j: /:11
Σ Vvl - -2) 15
Juna (lvyl A) /Jmax
0
Compute W E RexP elementwise as
Wikk (Wjk) max([Wje] 0,
Output: & constrained acyclic matrix W
sorts the entries of W by absolute value, computes threshold 0 from the largest entries consistent with the & constraint, and then applies soft-thresholding to shrink them toward zero (setting some exactly zero) while preserving signs. The resulting matrix W satisfies [[]]e, S À by construction.
C Gradients
Proposition presents the gradients of the minimizer in (3) with respect to the continuous matrix W and the sparsity parameter X. The proof the end of this section.
Proposition 2. Let W (Wzz) RpxP be the output of projection (3) applied to the matrix W (W;x) RexP with sparsity parameter > Let A {(j, k) Wyk + 0} be the set of active edge weights. Then, if the & constraint is binding, the gradients of Wik with respect to War and X are
sign(w (g sign(W;;)
Ow,k Ni gr
ggr jk
if(,k) A
A
otherwise
17

---

and
sign(W5k)
Owyk 0x
((,K)E
A
otherwise,
(3). The KKT conditions for stationarity and where 0% :- k[(ik) (g. r] and [A] is the cardinality of A. Otherwise, the constraint is non-binding, the gradients are
Jwyk ÖWier
((,,))
otherwise
and
Ow,k
0.
0X
Proposition 2 indicates that the Jacobian is sparse and has nonzeros only in positions corresponding to the surviving elements in W. In the binding case, where the & constraint has a shrinkage effect, the value of the nonzeros depends on the size of the active set A (the number of nonzero edges). In the non-binding case, where the { constraint imparts no shrinkage, all nonzero components equal one. Importantly, Proposition indicates that only the projection output W is needed evaluate the gradients; no additional computations are required. This result is useful for posterior learning since it means that a backward pass through the model demands virtually no further overhead cost.
Thompson et al. (2024) presented a result related to Proposition in which they provided the gradients of DAG projection for contextual learning. Their gradients differ slightly from those here because their sparsity constraint different. The gradients relating to were not considered in their work.
Proof. Let W' denote an optimal solution to projection complementary slackness are
(1
ll w + (IWW"l A) hh((w))
(8)
and
v¹ (I lle. - =0
(9)
where v* and 17* are the optimal values of the KKT multipliers. The optimal solution W* is implicitly a function of the problem data W and X, i.e., W = W* (W, A), as are the optimal KKT multipliers v* = <W, A) and 17" (W,, A). Denote the set of nonzeros by A.
Non-binding & constraint Suppose the & constraint is non-binding. Evaluating the subderivative on the left-hand side of the stationarity condition (8) gives
Wfk Wsk 0,
where use that D* and (l(W))/aws 0. Differentiating this expression with respect
Wer for ki,k)= A, gives
Dwž 59 NWgr
For those kj.k) & A. wfk = and hence
Owi 0. Digr
Moreover, because 2* = 0, we have
Ow
0x
for all (j, k).
18

---

Binding & constraint: active set Suppose the & constraint is binding and consider the gradients
for Wis in the active For these (j. k) E A, it follows from the stationarity condition (8) that
Wiz Wik + 1/ sign(w);) = O,
(10)
where we again use that Ol(W)) = 0. The complementary slackness condition (9) for (i, k) E A gives
Σ
|wjal -
(11)
(AA)EA
We proceed by first deriving the gradients of wik with respect to Wor- Differentiating (10) with respect Wer gives
Ow jk Niggr
Ov
v² θ sign(w3<) 0,
sign(W}()
Oviggr
from which it follows
Ov sigu(W]k)) Oiigr
(12)
DWgr
Next, differentiating (11) with respect to War gives
Ow Dwgr gr
Σ sign(w){)
(13)
(JR))AA
Substituting (12) into (13) and solving for Ov* 10 yields
Ov-
wign(w")))
sign(war)
(14)
DWer
AA
Finally, substituting (14) back into (12) gives
sign (en (Wy.)) sign(w;3)
8
Owir
A
We turn now to the gradients of wiR with respect to A. Differentiating (10) with respect to X gives
Owik ax
Ov+
sign (wgn(W5h) + vz
sign(w)() O,
ax
ax
from which it follows
Ow jk
On*
sign(w
(15)
01
ax
Next, differentiating (11) with respect to À gives
Owz
Σ
= 0.
(16)
0x
(.A)EA
Substituting (15) into (16) and solving for dv* /0x yields
Ov*
(17)
0x
A
Finally, substituting (17) back into (15) gives
Öw
sign(w),)
0x
A
Binding & constraint: inactive set Suppose the & constraint is binding and consider the gradients for wix in the inactive set. These (j. k) & A play no role in either of the KKT conditions, and hence
have no gradient with respect to Wik or A:
Ōw jk
0
and
Ow
0.
ax
19

---

D Proof of Proposition
Proof. The acyclicity constraint forces ||a/s ll2 for some j and k, meaning the whole vector "ij 0. Likewise, the sparsity constraint, which represents group lasso regularizer (Yuan and Lin, 2006), is also known to set some w/e zero. Those surviving vectors that not set to shrunk towards zero by a factor |/ l / The combined effect of these two constraints gives
1 l2
w/th
which equals zero if || 1illl2 = Now, observe that the objective function can be expressed solely in terms of the norms of each vector:
Σ j.k
will3
|1o0 gkk ll
k.k
Iwefi ll
12
Σ
W/
1
L ll2
3.k
2
|wefia ll
Σ
| 11
1
| ll
3,
(1
ll2
1
3.k
Taking Wyk- = || and wyk- = |W yields the claim of the proposition.
E
Implementation
ProDAG uses a multivariate Gaussian prior W. with independent elements, each having mean zero and variance one. The variational posterior W is also taken multivariate Gaussian with independent elements, where the means and variances are initialized at their prior values. We use the Adam optimizer (Kingma and Ba, 2015) with a learning rate of 0. minimize the variational objective function. The posterior variances are held positive during optimization using softplus transform. At each Adam iteration, 100 samples of W are drawn from the posterior to estimate the objective function. To project these sampled W and obtain W, we solve the projection using gradient descent with learning rates of 1/p and 0.25/p for the linear and nonlinear settings, respectively. Similar to DAGMA, the resulting W can contain small but nonzero elements that violate acyclicity. These small nonzeros occur because gradient descent seldom produces machine-precision zeros. To that end, we follow standard practice and set small elements of W to zero with a threshold of 0. At inference time, this threshold can be set such that is guaranteed acyelic (see, e.g., Ng et al., 2020).
Our implementation of ProDAG uses the machine learning library Flux (Innes et al., 2018) 8) and per- forms projections using parallel GPU (batched CUDA) implementations of the projection algorithms.
For the benchmark methods, we use the respective authors' open-source Python implementations:
DAGMA: haaatttttttr.ggggtcoooco////e/ea Apache License 2.0;
DiBS and DiBS+: http:://gggggggttooo...////////s v1.3.3, MIT License;
BayesDAG:
http:///gtaaooooo//tttrootttttpe v0.1.0, MIT Li-
cense; and
BOSS: http:://gghaa..c//mmmmappp//llat v7.6.6, GNU General Public License v2.0.
We use 50 particles for DiBS's particle variational inference. More particles exceed available memory.
20

---

The sparsity parameter of each method is tuned over grid of 10 values between A and Amax,5
ProDAG: > min 0 and Amax is the average & ball from the learned posterior with À = oc;
DAGMA: \ A 10-3 and A = 1;
DiBS and DiBS+: do include sparsity parameter;
BayesDAG: X A 10 and Amax 103; and
BOSS: does not include a sparsity parameter.
The above grids are chosen so that the sparsity level matches the true sparsity level for some value between A and A Amax We choose the specific value of À using a separate validation set of size [0.1n]. We cannot extract the weighted adjacency matrices or neural networks from BayesDAG (only the binary adjacency matrices), we choose its À that it is closest the true sparsity level.
In the nonlinear setting, the neural networks for ProDAG and DiBS consist of a single hidden layer with 10 neurons that use ReLU activation functions. DAGMA employs the same architecture but the neurons use sigmoid activations as its implementation does not support custom activation functions. The network for BayesDAG consists of two hidden layers, each with 128 neurons that use ReLU activations. BayesDAG's 's implementation does not support custom numbers of layers or neurons.
The experiments are run Linux workstation with an AMD Ryzen Threadripper PRO 5995WX CPU, 256GB RAM, and NVIDIA GeForee RTX 4090 GPUs. Each method allocated either single GPU single CPU core, with the experiments run in parallel across GPUs/CPU cores.
F Simulation design
The synthetic datasets generated from linear DAGs are constructed as follows. We first randomly sample an graph are then oriented according to a randomly generated topological ordering. Each directed edge in the graph is then weighted by sampling weights uniformly on [-0.7, -0.3] U [0.3, 0.7], yielding weighted adjacency matrix W. We then draw iid noise En and generate the variables as
or scale-free graph with nodes and edges. The edges of this undirected
I₁ W + Ess i i= n.
For synthetic datasets generated from nonlinear DAGs, we likewise sample Erdós Rényi scale- free graph and orient it according to a random topological ordering. We then generate a feedforward neural network fa with a single hidden layer of 10 neurons for k 1, :,....P. p. The weights of fa sampled uniformly from [-0.7, - 0.3) U [0.3, 0.7]. with hidden layer weights of node k's non-parents set to zero. The activation is ReLUU The noise & -
Enp is drawn iid and the variables taken as
K; = f(x;) + Es. i= 1=1,.... 1,...,,,
where (fi)) (f.) fp()).
The noise distribution is specified follows. For the Gaussian experiments, we draw p independent components from a Normal(O, 1) distribution. For the exponential experiments, the components are drawn from an Exponential(1) distribution, and for the Gumbel experiments from a Gumbel(0. 1) distribution. For heteroscedastic Gaussian experiments, we draw Eıj (the jth component of Ex) from Normal of), where oj is sampled from Uniform(2/33 4/3) and held fixed across observations.
The Brier score, used as a metric in the experiments, is defined as
Brier := Σ(1 * 0 ?
J.k
where 1 [w * 0] is the indicator that the true DAG contains a directed edge from node node k and Pizk is the posterior probability that this edge exists, as estimated from posterior samples.
'The parameter X in ProDAG refers to the é ball size, whereas À in DAGMA and BayesDAG refers to the coefficient on the & penalty.
21

---

G Additional experiments
G.1 Alternative simulation designs
We investigate the robustness of ProDAG to several alternative simulation designs. Figures and report results on dense graphs containing two and three times many edges as nodes. Figures and 10 report results on scale-free graphs whose degree exponent set to two and three. Figures 11 and 12 report results from Gumbel noise distributions, while Figures 13 and 14 report results from exponential noise distributions. Finally, Figure 15 report results from Gaussian noise distributions with heteroscedastic variance. A Gaussian likelihood is used consistently across settings. Under all simulation designs, the relative performance of ProDAG roughly similar to that the main experiments, indicating that ProDAG remains reliable across broad range of realistic scenarios.
G.2 Semi-synthetic datasets
We generate semi-synthetic datasets using the MUNIN (subnetwork #1) DAG, a large medical diagnos- tic network from the Bayesian Network Repository (https: //www//wwwwww bnlearn com/bnreppn wwwpbbeeeiiittoooo torry. The graph, which contains p = 186 nodes and = 273 edges, is used to generate semi semi-synthetic datasets by sampling edge weights and noise in the same fashion as the fully synthetic experiments. Table 2 reports results with sample size n = 100, showing that ProDAG performs strongly with the lowest Brier score and expected SHD. While DAGMA scores well on expected F1 score and AUROC,
Table 2: Performance on semi-synthetic datasets of size = 100 generated from the MUNIN graph from the Bayesian Network Repository with p 186 nodes, $ 273 edges, and Gaussian noise. The averages and standard eтrors are measured over 10 independently and identically generated datasets. The best value of each metric is indicated in bold.
Brier
Exp. SHD
Exp. F1 score (%)
AUROC (%)
ProDAG
146 + 4.0 305 65.8
201 4.9
37 + 2.0
87 = 0.6
DAGMA
274 64.9
67 4.1
0.5
DiBS
DiBS+
BayesDAG
279 = 0.3
557 1 0.9
1.0
53 = 0.7
BOSS
ProDAG provides superior calibration and structure recovery. To our knowledge, this experiment represents one of the largest Bayesian DAG evaluations reported to date. Aside from ProDAG, the only methods successfully run are DAGMA and BayesDAG. The other baselines encounter difficulties: DiBS and DiBS+ produce nonsensical results indicative of convergence issues, while BOSS crashes.
G.3 MCMC baseline
We compare ProDAG with Gadget (Viinikka al., 2020), earlier Bayesian approach that learns posteriors over DAGs via MCMC. We focus on linear Erdōs-Rényi DAGS with p = 20 nodes and 40 edges, Gadget does not scale larger setting. Table 3 reports results with sample size n 100, showing that Gadget substantially underperforms ProDAG. This underperformance is not
Table 3: Performance on synthetic datasets of size n 100 generated from linear Erdós Erdös-Rényi DAGs with p 20 nodes, 40 edges, and Gaussian noise. The averages and standard errors are measured over independently and identically generated datasets. The best value each metric indicated in bold.
Brier score
Exp. SHD
Exp. F1 score (%)
AUROC (%)
ProDAG Gadget
13 + 1.8
14 + 1.5
78 + 2.7
91 + 1.6
52 1 1.8
37 1 0.5
17 1.5
65 = 3.0
due to Gadget's configuration, evaluated many settings and report the best-performing
22

---

Expected F1 score
0.55
80
Biier Score
Expected
AUROC
0.90
0.35
030
101
Sample
Sample
Sample
Sample
Figure 7: Performance synthetic datasets generated from linear Cedös-Rényi DAGs with p 20 nodes, s = 60 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) measured over 10 independently and identically generated datasets.
PrODAC: DRS DAyesAAG
DAGMM DBBS+ BOSS
120
1.00
102
ExpectedF1 score
0.75
Biner score
AUROC
88
Expected
0.90
06
120
123
125
12:44
7044
109
124
70m
105
12'4
1021
1024
100
1234
7000
129
Sample size
Sample size
Sample size
Semple size
Figure 8: Performance on synthetic datasets generated from linear Erdös Frdös-Rényi DAGs with p = 20 nodes, 80 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
DBS BayenDAG
DAGMM 0BBC- 8OSS
Expected F1 score
Biner score
AUROC
Expected
00
Sample size
Sample size
Sample size
Sample size
Figure 9: Performance on synthetic datasets generated from linear scale-free (degree exponent 2) DAGs with = 20 nodes, s 40 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
23

---

00
Expected F1 score
0.75
Brier Score
AUROC
08
0.30
0.25
030
1014
Sample
Sample
Sample
Sample
Figure 10: Performance synthetic datasets generated from linear scale-free (degree exponent 3) DAGs with = 20 nodes, = 40 edges, and Gaussian noise. The averages (solid points) and standard errors (error bars) measured over 10 independently and identically generated datasets.
PrODAG: DRS DAyesAAG
DAGMA DBBS+ BOSS
100
1.00
75
ExpectedF1 score
0.75
Biner score
Expected
AUROC
88
0.90
0.
120
123
12:4
7014
109
124
7099
70m
123
72'4
102
1024
100
12.4
7000
129
Sample size
Sample size
Sample
Semple size
Figure 11: Performance on synthetic datasets generated from linear Rényi DAGs with p = 20 nodes, 40 edges, and Gumbel noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
DayesDAG
DAGMA DJBS+ BOSS
Expected F1 score
Biter score
AUROC
8.7
350
0.00
Sample size
Sa Sample size
Sample size
Sample size
Figure Performance on synthetic datasets generated from linear Erdös Erdös-Rényi nodes, s = 200 edges, and Gumbel noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
Rényi DAGs with p 100
24

---

Expected F1 score
Brier Score
AUROC
88
0.30
25
0.25
030
10
100
Sample
Samplc size
Sample
Sample
Figure 13: Performance synthetic datasets generated from linear Erdós-Rees DAGs with p 20 nodes, = 40 edges, and exponential noise. The averages (solid points) and standard errors (error bars) measured over 10 independently and identically generated datasets.
PrODAC a DISS DAGAADAS
DAGMA DISS+ BO5S
1000
1.00
750
750
ExpectedF1 score
0.75
Biner score
Expected
AUROC
590
0.90
0.5
1244
1000
7231
7014
104
1044
1020
1220
1250
72"4
104
1004
7091
100
72.4
7090
129
Sample size
Sample
Sample size
Semple size
Figure 14: Performance on synthetic datasets generated from linear Rrdös-Rényi DAGs with 100 nodes, 200 edges, and exponential noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
DBS
BayenDAG
DAGMM D.BC- 8OSS
Expected F1 score
Biner Score
AUROC
--
Expected
0:95
Sample size
Sample size
Sample size
Sample size
Figure 15: Performance on synthetic datasets generated from linear Erdös-Rémt Rényi DAGs with p = 20 nodes, = 40 edges, and heteroscedastic Gaussian noise. The averages (solid points) and standard errors (error bars) are measured over 10 independently and identically generated datasets.
25

---

G.4 Timings
We compare the run times of ProDAG with those of the Bayesian baselines in our main experiments. Table 4 reports timings on (p = 20 nodes, 40 edges) and nonlinear (p 10 nodes, = 20 edges) Erdös-Rény: DAGs with Gaussian noise. The timings show that ProDAG is faster than BayesDAG
Table Run times in seconds with a sample size n 100 for linear (p = 20 nodes) and nonlinear (p 10 nodes) settings. The averages and standard errors are measured over 10 independently and identically generated datasets. The best value in each setting is indicated in bold.
Linear
Nonlinear
ProDAG DiBS
30 + 7.0
61 + 3.8
8
+ 0.4
BayesDAG
57 + 1.1
122 + 5.5
and slower than DiBS. However, ProDAG's superior performance in uncertainty quantification and structure recovery makes a strong choice, even with the higher computational cost over DiBS.
H
Limitations
Though approach has several favorable properties, it naturally has limitations, the main highlight here. As with most Bayesian approaches that employ variational inference, ours does not guarantee that the learned variational posterior approximates the posterior well. The possibility of a poor variational approximation is further complicated by the nonconvexity of the acyclicity constraint, giving rise nonconvex optimization problem that can be solved only a stationary point. We remark, however, that nonconvexity is a property of most existing DAG learning methods.
I Societal impacts
DAG learning tools play key role in causal inference, where the goal is to establish causal rela- tionships between random variables. Our approach benefits from quantifying uncertainty in this process, which can help guide downstream decision-making However, as with any DAG learning tool applied to observational data, cannot causally interpret the DAGs learned by approach without making empirically unverifiable assumptions such as causal sufficiency. Therefore, it is crucial to involve subject-matter experts capable of determining the validity of any such assumptions.
26


Reviews:
[{"review_id": "qE3F0kZ2Cs", "reviewer": "Reviewer_e12N", "summary": "This paper introduces ProDAG, a Bayesian structure learning method that employs variational inference to infer an approximate posterior over DAGs given observational data. Unlike prior variational approaches which either (i) enforce acyclicity in the variational family or (ii) modify the prior to penalize cyclicity, ProDAG takes the novel approach of choosing an unconstrained variational family that induces a distribution over DAGs through the use of a many-to-one projection operation. It is shown that this defines a valid distribution over DAGs (weight matrices), and how to propagate gradients through the projection operation. Empirical results show that ProDAG outperforms existing variational Bayesian structure learning methods on standard metrics, particularly as the sample size increases.", "questions": "- Intuitively, why is it possible to derive the gradients of the optimal solution (Apx C) in closed form but not the optimal solution itself? Given that the iterative procedure only produces an approximation solution to the projection, is there a significant discrepancy between the exact autograd gradients for the approximate solution and the gradients of the optimal solution?\n- Why did the authors choose a multivariate Gaussian variational posterior, given that the true posteriors can be multimodal? (unless the projection step induces multimodality). It would be interesting to see some examples in low dimensions to visualize the effect of the projection step.\n- Did the authors try learning a variational posterior with a Laplace-style prior instead of the hard sparsity projection? If so and it performed worse, are there some qualitative insights into why the hard projection step should be preferred?\n- How did the authors tune the sparsity hyperparameter $\\lambda$ on the validation data? How sensitive is the performance of the method to the choice of $\\lambda$ (i.e. would it be reasonable to simply choose a single $\\lambda$ value and use it across all experiments)?", "limitations": "Yes, limitations are addressed in the appendix.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Overall, this is a high quality paper introducing an interesting approach for Bayesian structure learning. The hard projection step onto weight matrices with a particular sparsity is a novel idea and is shown to lead to strong empirical performance. The idea is backed up with solid theoretical work showing the measurability of the projection mapping and derivation of the ELBO and gradients. The clarity of the paper is also excellent with the key ideas being described well in the main paper with full details presented in the Appendix. \n\nRegarding weaknesses, while the empirical evaluation shows strong performance, I would like to have seen some more insight into *why* ProDAG performs better than existing approaches (see questions below). Additionally, the empirical evaluation compares to variational approaches such as BayesDAG and DiBS that employ augmented spaces. However, for completeness there should also be comparison to other types of Bayesian methods, in particular (i) permutation-based variational methods such as BCDNets; and (ii) MCMC-based approaches such as Gadget [1] which may exhibit qualitatively different behavior (e.g. consistency of MCMC). The related work section should also mention some of these MCMC-based approaches. Finally, a seeming limitation is that ProDAG requires the sparsity hyperparameter $\\lambda$ to be set in advance, and due to the hard projection constraint this may inevitably lead to errors even in the large-data limit (i.e. if one sets $\\lambda < d(d-1)/2$ then one could never learn the complete graph even if it is the true graph). \n\n\n[1] Towards Scalable Bayesian Learning of Causal DAGs. Viinikka et al. (2020). NeurIPS", "was_revised": true, "final_justification": "The authors have fully addressed my concerns regarding experimental evaluation, and I support acceptance of the paper.", "rebuttal": "We appreciate your time and attention in reviewing our work.\n\n## Weaknesses\n\n> Regarding weaknesses, while the empirical evaluation shows strong performance, I would like to have seen some more insight into why ProDAG performs better than existing approaches (see questions below).\n\n`ProDAG` places both the prior and variational posterior directly on DAG space, enforcing exact acyclicity through projection. Competing methods (i) operate in relaxed spaces with soft penalties that cannot guarantee acyclicity, or (ii) introduce discrete variables (e.g., permutations) that complicate optimization. See the discussion in Section 2. By pushing irrelevant edges to exact zeros while keeping the ELBO differentiable, `ProDAG` yields sharper uncertainty estimates and lower SHD scores.\n\n> Additionally, the empirical evaluation compares to variational approaches such as BayesDAG and DiBS that employ augmented spaces. However, for completeness there should also be comparison to other types of Bayesian methods, in particular (i) permutation-based variational methods such as BCDNets; and (ii) MCMC-based approaches such as Gadget [1] which may exhibit qualitatively different behavior (e.g. consistency of MCMC). The related work section should also mention some of these MCMC-based approaches.\n\nFor permutation-related variational methods such as `BCDNets`, we note that we have, in fact, compared to `BayesDAG` which is representative of the state-of-the-art in this class of methods and has been shown in its original paper to exceed `BCDNets`'s performance. See Section 3.2 of the `BayesDAG` paper for the equivalence to permutation approaches and Section 6 for the comparisons against `BCDNets`. For MCMC-based methods, we have now added comparisons against `Gadget` with $p=20$ and $n=100$. **Table R1** shows it substantially underperforms `ProDAG`. This underperformance is not due to `Gadget`’s configuration, as we evaluated many settings and report the best-performing one. We will expand the related work section to include a discussion of `Gadget` and MCMC-based methods.\n\n**Table R1** Comparisons against `Gadget`.\n\n|          | Brier score | Exp. SHD | Exp. F1 score (%) | AUROC (%) |\n|:---------|:------------|:---------|:------------------|:----------|\n| `ProDAG` | 13 ± 1.8    | 14 ± 1.5 | 78 ± 2.7          | 91 ± 1.6  |\n| `Gadget` | 52 ± 1.8    | 37 ± 0.5 | 17 ± 1.5          | 65 ± 0.5  |\n\n> Finally, a seeming limitation is that ProDAG requires the sparsity hyperparameter $\\lambda$ to be set in advance, and due to the hard projection constraint this may inevitably lead to errors even in the large-data limit (i.e. if one sets $\\lambda<d(d-1)/2$ then one could never learn the complete graph even if it is the true graph). \n\nWe treat $\\lambda$ purely as a regularisation knob and apply a grid search, so the tuning procedure (see our response below) selects the $\\lambda$ best supported by the data. The grid search includes $\\lambda=\\infty$ so the true graph, however dense, is never ruled out. Only an intentionally small $\\lambda$ would exclude it.\n\n## Questions\n\n> Intuitively, why is it possible to derive the gradients of the optimal solution (Apx C) in closed form but not the optimal solution itself? Given that the iterative procedure only produces an approximation solution to the projection, is there a significant discrepancy between the exact autograd gradients for the approximate solution and the gradients of the optimal solution?\n\nThe projection has no closed-form solution because it solves a high-dimensional non-convex problem, but once a stationary point $W^\\star$ is attained its KKT conditions are algebraic relations that can be implicitly differentiated, giving the compact gradient formulas in Proposition 2 (in short, local sensitives are easier to obtain than global minimizers). The gradient formulas are evaluated at the $W^\\star$ returned by the projection, so they exactly coincide with the gradient autograd would produce at the same point, i.e., there is no discrepancy.\n\n> Why did the authors choose a multivariate Gaussian variational posterior, given that the true posteriors can be multimodal? (unless the projection step induces multimodality). It would be interesting to see some examples in low dimensions to visualize the effect of the projection step.\n\nWe choose a Gaussian in the unconstrained space $\\tilde{W}$ for differentiable reparameterization, closed-form KL, and efficient optimization. We stress, however, that the projection $\\tilde{W}\\mapsto W$ does indeed induce a non-Gaussian, *multimodal posterior* over DAGs. We will clarify this point explicitly in the revised text. We are not allowed to include visualizations in the rebuttal, but Figure 1 in the paper roughly illustrates the effect of the projection.\n\n> Did the authors try learning a variational posterior with a Laplace-style prior instead of the hard sparsity projection? If so and it performed worse, are there some qualitative insights into why the hard projection step should be preferred?\n\nYes, during prototyping we evaluated the DAG projection under a Laplace prior and saw noticeably higher SHD scores because Laplace only shrinks edges toward zero and does not set them exactly to zero, so the resulting graphs were too dense. In contrast, the $\\ell_1$ ball prior assigns positive probabilities on exact zeros and thus yields crisper sparsity patterns and better edge recovery. See Xu & Duan (2023, arXiv:2006.01340) for a detailed comparison of $\\ell_1$ ball versus Laplace. \n\n> How did the authors tune the sparsity hyperparameter $\\lambda$ on the validation data? How sensitive is the performance of the method to the choice of $\\lambda$ (i.e. would it be reasonable to simply choose a single $\\lambda$ value and use it across all experiments)?\n\nWe sweep $\\lambda$ over ten values from $\\lambda_\\min=0$ to $\\lambda_\\max$ ($\\lambda_\\max$ is the average $\\ell_1$-norm of $W$ from the learned posterior when $\\lambda=\\infty$) and pick the value that maximizes the expected likelihood on a held-out 10% validation set (see Appendix E). Because $\\lambda$ directly controls graph sparsity, performance drops if it is set too low (underconnected) or too high (overconnected), so exploring a few points on a grid is preferable to a single value. We apply the same grid‑search to every method with a sparsity parameter."}, {"review_id": "reR2s9Szh4", "reviewer": "Reviewer_VF9w", "summary": "This paper utilize Bayesian variational inference framework to construct a valid distribution on the space of sparse DAGs for the causal discovery task. The proposed method ProDAG models the prior and variational posterior of potential DAGs. For this purpose the authors use the continuous reformulations of acyclicity constraints and perform gradient descent based on that.", "questions": "Below I provide my questions.\n\n* How do the author set the value for $\\lambda$.\n* The authors said that the projection (2) is unique. Why can there not be two different $\\tilde{W}$ with equal distance from $W$?\n* The authors mentioned that step 1 projects the matrix $\\tilde{W}$ onto the set of acyclic matrices through equation 4. How does the algorithm have access to $W$ to calculate the subtraction here?", "limitations": "yes.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Below I provide my comments.\n\n# Strengths\nFigure 1 is very helpful to understand the algorithm concepts. The idea of using variational inference for causal discovery appears novel.\n\n# Weaknesses\n\n* The following works also considers a continuous constraint optimization problem\nThe authors should compare their methods to the following works qualitatively and use some of them as benchmarks.\n\n1. Brouillard, Philippe, et al. \"Differentiable causal discovery from interventional data.\" Advances in Neural Information Processing Systems 33 (2020): 21865-21877.\n2. Bhattacharya, Rohit, et al. \"Differentiable causal discovery under unmeasured confounding.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2021.\n3. Faria, Gonçalo Rui Alves, Andre Martins, and Mário AT Figueiredo. \"Differentiable causal discovery under latent interventions.\" Conference on Causal Learning and Reasoning. PMLR, 2022.\n4. He, Yue, et al. \"Daring: Differentiable causal discovery with residual independence.\" Proceedings of the 27th ACM SIGKDD conference on knowledge discovery & data mining. 2021.\n5. Lopez, Romain, et al. \"Large-scale differentiable causal discovery of factor graphs.\" Advances in Neural Information Processing Systems 35 (2022): 19290-19303.\n\n* The authors should provide more details on the “path-following algorithm”.\n* The authors should provide some intuition on “implicit function theorem.” and “Karush-Kuhn- Tucker (KKT) conditions of the projection” in the main paper.\n* The role of $l_1$ ball is not clear. The author should discuss how large or small $l_1$ ball affects the algorithm performance. Also it is not clear how $l_1$ is chosen.\n* It is not clear how the authors obtain $q_{\\theta}(\\tilde{W})$ with variational inference in their algorithm. I would suggest the authors discussing this in detail.\n* The authors should show more experimental results with real-world data. They might check https://www.bnlearn.com/bnrepository/ for this purpose.\n* Figure 2 might be replaced with a table and the space can be used for describing other concepts in detail.", "was_revised": true, "final_justification": "This work has some novel contributions and the authors also added new experiments according to my suggestion. Thus, I keep my positive score.", "rebuttal": "Thank you for your thoughtful review and engagement with the paper.\n\n## Weaknesses\n\n> The following works also considers a continuous constraint optimization problem The authors should compare their methods to the following works qualitatively and use some of them as benchmarks.\n\nThank you for these references. We agree they relate to our work but target different settings and all are non-Bayesian. Brouillard et al. (2020), Faria et al. (2022) and Lopez et al. (2022) rely on interventional data, whereas `ProDAG` assumes purely observational data. Bhattacharya et al. (2021) also uses observational data but learns acyclic directed mixed graphs (graphs with bidirected edges), not DAGs. He et al. (2021) adds a residual‑independence regularizer to frequentist `NOTEARS`‑type models, but it is unclear how to embed this within `ProDAG`’s Bayesian projection prior. We will cite the papers and discuss the extension of `ProDAG` to interventional data as future work.\n\n> The authors should provide more details on the “path-following algorithm”.\n\nWe agree that the path-following algorithm (presented in Appendix B) would benefit from additional explanatory text to complement the algorithmic details. We will revise the appendix to include a more descriptive explanation of the algorithm’s steps.\n\n> The authors should provide some intuition on “implicit function theorem.” and “Karush-Kuhn- Tucker (KKT) conditions of the projection” in the main paper.\n\nWe agree. While the full details were placed in Appendix C due to space limitations, we will include descriptions in the main text in the final version.\n\n> The role of $l_1$ ball is not clear. The author should discuss how large or small $l_1$ ball affects the algorithm performance. Also it is not clear how $l_1$ is chosen.\n\nThe $\\ell_1$ ball constraint is included to sample sparse graphs; a larger $\\ell_1$ ball mean more edges while a smaller $\\ell_1$ ball means fewer edges. In our numerical experience, the algorithm's speed is not particularly sensitive to the size of the ball, though smaller balls tend to converge slightly faster than larger balls. The size $\\lambda$ of the ball is chosen using a validation set (see Appendix E).\n\n> It is not clear how the authors obtain $q_\\theta(\\tilde{W})$ with variational inference in their algorithm. I would suggest the authors discussing this in detail.\n\nThe distribution $q_\\theta(\\tilde{W})$ is the variational family chosen by the user (see the last paragraph of Section 5.1) such as a multivariate Gaussian, and its parameters $\\theta$ are optimized via gradient descent to maximize the ELBO. The projected distribution over DAGs, $q_\\theta(W)$, is implicitly defined by applying the projection to samples from $q_\\theta(\\tilde{W})$.\n\n> The authors should show more experimental results with real-world data. They might check bnrepository for this purpose.\n\nThank you for the reference to this great repository. We have now added **Table R2**, which reports results on the Andes network ($p=223$) from the repository, drawing $n=100$ samples from its DAG with our usual data generation pipeline. This large, real‑world‑structured benchmark shows that `ProDAG` performs strongly on very high-dimensional DAGs, achieving the lowest Brier score and SHD. While `DAGMA` scores well on F1 and AUROC, `ProDAG` provides superior calibration and structure recovery.\n\n**Table R2** Experiments on real-world Andes DAG with $p=223$ nodes and $s=338$ edges. We report 8 runs here due to time constraints and will extend to 10 in the final version, but results are already clearly differentiated.\n\n|            | Brier score | Exp. SHD   | Exp. F1 score (%) | AUROC (%) |\n|:-----------|:------------|:-----------|:------------------|:----------|\n| `ProDAG`   | 143 ± 3.4   | 198 ± 4.7  | 39 ± 2.0          | 87 ± 0.8  |\n| `DAGMA`    | 337 ± 78.7  | 306 ± 77.8 | 65 ± 4.9          | 89 ± 0.6  |\n\nTo our knowledge, this experiment represents one of the largest Bayesian DAG evaluations reported to date. Aside from `ProDAG`, `DAGMA` was the only method to successfully run on this large dataset. The other Bayesian baselines encountered difficulties: `DiBS` and `DiBS+` produced nonsensical results (Brier scores and SHDs in the tens of thousands), indicative of convergence issues, while `BayesDAG` ran out of memory (24 GB available). `BOSS` encountered crashes.\n\n> Figure 2 might be replaced with a table and the space can be used for describing other concepts in detail.\n\nIf space is a concern in the final version, we will move Figure 2 to the appendix to make room for additional explanations in the main text.\n\n## Questions\n\n> How do the author set the value for $\\lambda$.\n\nThe parameter $\\lambda$ is chosen by maximizing the expected likelihood on a validation set. See Appendix E for details.\n\n> The authors said that the projection (2) is unique. Why can there not be two different $\\tilde{W}$ with equal distance from $W$?\n\nUniqueness means that for almost every draw of $\\tilde{W}$ the projection (2) has one minimizer, so $\\operatorname{pro}_\\lambda(\\tilde{W})$ is single‑valued for each fixed $\\tilde{W}$. This result does not imply injectivity: different $\\tilde{W}$ can still project to the same $W$, but no single $\\tilde{W}$ has more than one optimal projection. We will clarify this point in the text.\n\n> The authors mentioned that step 1 projects the matrix $\\tilde{W}$ onto the set of acyclic matrices through equation 4. How does the algorithm have access to $W$ to calculate the subtraction here?\n\nThe $W$ in (4) appearing inside $f_\\mu(W;\\tilde{W})=\\frac{\\mu}{2}\\|\\tilde{W}-W\\|_F^2+h(W)$ is the current optimization variable. The optimizer repeatedly updates this variable while it computes the difference $\\tilde{W}-W$ until $f$ is minimized. After convergence we reuse the same symbol to denote the minimizer, so the duplicate notation is purely stylistic and we never assume prior knowledge of the solution."}, {"review_id": "8CrBHABRXw", "reviewer": "Reviewer_agaZ", "summary": "This paper introduces a method for learning DAGs that addresses the challenging problem of uncertainty quantification in causal structure discovery. The main innovation is the development of \"projected distributions\" that have support directly on the space of sparse DAGs by projecting samples from continuous distributions onto the set of acyclic matrices using a combination of acyclicity and sparsity constraints. The authors prove these projected distributions are mathematically valid, develop efficient algorithms for the required projections, and demonstrate through a range of experiments that their method performs well in both accuracy and uncertainty quantification.", "questions": "1. Can you demonstrate the method on larger networks (p=200-500), or propose any algorithmic improvements that could \nsignificantly increase the scaleability of the approach?\n\n2. Can you determine which components are responsible for the performance improvements? Some ablations could help with this, such as  comparing against a simpler baselines with the same acyclicity constraints without projected distributions. \n\n3. Can you provide recommendations about the settings in which your approach should be chosen over existing methods, given the modest improvements and higher computational cost? \n\n4. Can you provide more analyses as to when the variational approximation works well? This could include diagnostics or an experiment comparing against exact posteriors on a small synthetic problem.", "limitations": "The limitations are discussed above. I do not see any potential negative societal impact of the work.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "This is a strong paper that makes solid theoretical and practical contributions to Bayesian DAG learning. The core idea of projected distributions addresses an important problem, as most existing Bayesian DAG methods either do not guarantee exact acyclicity or require complex discrete inference procedures. The method appears technically sound and shows improvements over baselines, though gains are often small. The projected distribution idea could likely be applied beyond DAGs to other constrained structure learning problems, increasing the significance of the work.\n\nStrengths:\n\nQuality: The technical execution is very good, with rigorous mathematical foundations and proofs, and comprehensive experimental evaluation in both linear and nonlinear SEMs with multiple metrics and baselines. The variational inference framework is well-designed and the extension to nonlinear DAGs via neural networks is natural, though the method still inherits standard limitations like its cubic complexity and approximate posteriors from variational inference.\n\nClarity: The paper is well-written, with the projected distribution concept presented intuitively and the theoretical framework clearly explained. The novel ideas are integrated nicely with known results.\n\nSignificance: This addresses an important problem in causal inference with practical implications for uncertainty quantification in downstream tasks. The method shows strong uncertainty quantification compared to frequentist approaches and performs well for large sample sizes where other Bayesian methods fall off, with particularly good results on the real protein signaling dataset.\n\nOriginality: The key innovation is to define a novel distribution P(W) over the space of DAG adjacency matrices, which can be used in a variational framework. The approach seems innovative compared to previous approaches, especially this contribution of defining this distribution with continuous support. The entire framework also combines the new ideas with known results/tools nicely. \n\nWeaknesses:\n\nQuality: The method inherits standard limitations including cubic complexity, which limits its scalability (experiments only go to p=100). The use of variational inference does not guarantee good approximation to the true posterior. The non-convexity of the acyclicity constraint means optimization can only reach stationary points rather than global optima.\n\nClarity: While generally well-written, some important experimental details are deferred to appendices rather than being in the main text, and the combination of multiple existing techniques (though novel in combination) makes it hard to work out which components drive the performance improvements.\n\nSignificance: The performance improvements over baselines are consistent but small, and the scalability limitations to networks with around 100 variables may limit practical applicability to larger real-world problems where the method's benefits would be most valuable.\n\nOriginality: While the overall framework is novel, the individual technical components such as the acyclicity constraints, L1-ball projections,  and variational inference) are all existing techniques, so the contribution is primarily in their creative combination rather than fundamental algorithmic innovations. Nonetheless there are some significant innovations in the work.", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to engage with our work and for your thoughtful feedback.\n\n## Weaknesses\n\n> The method inherits standard limitations including cubic complexity, which limits its scalability (experiments only go to p=100). The use of variational inference does not guarantee good approximation to the true posterior. The non-convexity of the acyclicity constraint means optimization can only reach stationary points rather than global optima.\n\nIt is correct that non-convexity means we can only attain a local optimum, though this property is true of all differentiable DAG learning methods (e.g., `DAGMA`). Regarding scalability, we now include larger experiments (see our response below).\n\n> While generally well-written, some important experimental details are deferred to appendices rather than being in the main text, and the combination of multiple existing techniques (though novel in combination) makes it hard to work out which components drive the performance improvements.\n\nWith the extra page allowed in the camera ready version, we will move additional experimental details (e.g., data generation and hyperparameter tuning) from the appendix to the main text. Regarding the source of performance gains, please see our response below concerning component attribution and ablations.\n\n> The performance improvements over baselines are consistent but small, and the scalability limitations to networks with around 100 variables may limit practical applicability to larger real-world problems where the method's benefits would be most valuable.\n\nWe agree that practical applicability depends on scalability. As noted in our response below, we now include results on the Andes network with $p=223$ in **Table R2** (see response to Reviewer VF9w), which demonstrate that `ProDAG` scales beyond $p=100$ while maintaining its performance advantages.\n\n> While the overall framework is novel, the individual technical components such as the acyclicity constraints, L1-ball projections, and variational inference) are all existing techniques, so the contribution is primarily in their creative combination rather than fundamental algorithmic innovations. Nonetheless there are some significant innovations in the work.\n\nThe individual tools (acylclicity characterizations, $\\ell_1$ balls) are indeed known, but their composition to define a projected prior is, as you say, new. No prior work samples exact sparse DAGs by projecting continuous draws and then instantiates that construction in a variational framework. This step is the key innovation that drives `ProDAG`'s gains.\n\n## Questions\n\n> Can you demonstrate the method on larger networks (p=200-500), or propose any algorithmic improvements that could significantly increase the scaleability of the approach?\n\nWe now include experiments on the $p=223$ Andes network from bnlearn (**Table R2**). On this larger, real‑world‑structured benchmark, `ProDAG` continues to delver excellent performance and represents one of the largest Bayesian DAG evaluations reported to date.\n\n> Can you determine which components are responsible for the performance improvements? Some ablations could help with this, such as comparing against a simpler baselines with the same acyclicity constraints without projected distributions.\n\nWe believe that `DAGMA`, a key baseline in our experiments, reflects the ablation you propose: it enforces the same acyclicity constraints but omits projected distributions. `DAGMA` can also be interpreted as a MAP version of `ProDAG`. The consistent improvements over `DAGMA` underscore the value of modeling full posteriors directly over DAGs, made possible by our projected variational framework. In practice, these gains trace to two components: (i) projected DAG distributions that place mass on sparse, exactly acyclic structures, and (ii) variational inference that leverages this distributional support to propagate uncertainty during learning.\n\n> Can you provide recommendations about the settings in which your approach should be chosen over existing methods, given the modest improvements and higher computational cost?\n\nChoose `ProDAG` when one of these situations arise:\n\n1. You need well-calibrated posterior uncertainty (e.g., decision making, active experimentation);\n2. The number of observations is limited, so precise sparsity and exact zeros aid recovery; or\n3. Exact acyclicity and hard sparsity matter for interpretability (e.g., policy, biology).\n\n> Can you provide more analyses as to when the variational approximation works well? This could include diagnostics or an experiment comparing against exact posteriors on a small synthetic problem.\n\nWhile exact posterior comparisons via MCMC are nontrivial, we report Brier score, SHD, F1 score, and AUROC as direct assessments of the variational approximation’s structure learning capacity. To further stress-test this approximation, we added new higher-dimensional experiments on nonlinear and non-Gaussian data (**Tables R3–R5**), and the Andes network. In all cases, `ProDAG` maintains strong performance, demonstrating robustness of the approximation across a range of challenging settings."}, {"review_id": "GiDKsplTzM", "reviewer": "Reviewer_7A4S", "summary": "This paper proposes a novel Bayesian causal discovery method ProDAG for learning directed acyclic graphs (DAG) with uncertainty quantification. The prior on the adjacency matrix $W$ of the DAGs is induced by projecting $\\tilde W\\in\\mathbb R^{p\\times p}, \\tilde W\\sim\\mathbb P$ onto its closest neighbor in the set of adjacency matrices. The author showed that the projection is unique and thus the prior is well-defined. Furthermore, the author introduced a variational inference method to learn the underlying DAG based on this prior. Experiments have shown that ProDAG performs competitively against the state-of-the-art (Bayesian) DAG learning methods.", "questions": "1. Do you assume knowledge of the true noise distribution for experiments on non-Gaussian data?\n2. Is ProDAG robust if the log-likelihood function is misspecified? For example, using a Gaussian likelihood for data generated by non-Gaussian noise distributions, or an equal variance likelihood for a heterogeneous non-Gaussian noise distribution?\n3. How sensitive are the results to the choice of $\\mathbb P$?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n1. The proposed method is novel and shows good performance for various noise distributions and graph structures.\n2. The prior distribution is well-motivated and justified.\n3. The projection procedure in the prior evaluation is efficient. Moreover, the author provides explicit gradient estimation for this projection.\n\n**Weaknesses**\n1. Identifiability. The author did not discuss what kind of DAGs the proposed model can identify. Since the VI procedure requires knowledge of the likelihood, it seems that the noise distribution has to be known to guarantee identifiability, which is a limitation that needs to be declared. \n2. Some experiments are not sufficient. For example, the experiment for nonlinear data is limited to $d=10$, and those for non-Gaussian noise distributions are limited to $d=20$. Experiments on higher dimensions, e.g., $d=100$ would be helpful.\n3. Time comparison with other Bayesian methods is not provided.\n\nI will be happy to increase the score if my concerns can be addressed.", "was_revised": true, "final_justification": "The author has clarified the experimental setting and provided more extensive experiments which addressed most of my concerns. My remaining concern is whether the method is robust under heterogeneous noise distribution. However, given its solid methodology, I find it not a reason to reject the paper. That is, even if it is not the case (but the author should mention this), the paper still provide a solid variational inference method for causal discovery in several scenarios. Therefore, I raise my score to 4.", "rebuttal": "We are grateful for your engagement and the constructive feedback on our submission.\n\n## Weaknesses\n\n> Identifiability. The author did not discuss what kind of DAGs the proposed model can identify. Since the VI procedure requires knowledge of the likelihood, it seems that the noise distribution has to be known to guarantee identifiability, which is a limitation that needs to be declared. \n\nThank you for highlighting this point. Indeed, identifiability depends on both the data-generating SEM (including the noise distribution) and the user-specified likelihood. For instance:\n\n- If the likelihood is correctly specified (e.g., Gaussian likelihood applied to Gaussian noise), the true DAG (or its equivalence class) is theoretically identifiable. In our experiments (e.g., Figure 3), we observe this behavior as the posterior increasingly concentrates on the true structure as $n$ grows; and\n- If the likelihood is misspecified (e.g., Gaussian likelihood applied to non-Gaussian noise), full identifiability may no longer hold. In this scenario, `ProDAG` will reflect the resulting uncertainty in its posterior.\n\nThus, `ProDAG` captures and reports the inherent uncertainty determined by the combination of the true SEM and the chosen likelihood. We will explicitly clarify this subtlety in the revised manuscript.\n\n> Some experiments are not sufficient. For example, the experiment for nonlinear data is limited to $d=10$, and those for non-Gaussian noise distributions are limited to $d=20$. Experiments on higher dimensions, e.g., $d=100$ would be helpful.\n\nWe added the following experiments with more nodes $p$ ($d$ in your notation) to address your concerns:\n\n1. Nonlinear data with $p=20$ (double the previous dimension) in **Table R3**;\n2. Gumbel noise data with $p=100$ in **Table R4**; and\n3. Exponential noise data with $p=100$ in **Table R5**.\n\nDue to NeurIPS image restrictions this year, we report a snapshot of the experiments at $n=100$.\n\n**Table R3** Larger nonlinear experiment with $p=20$.\n\n|            | Brier score | Exp. SHD   | Exp. F1 score (%) | AUROC (%) |\n|:-----------|:------------|:-----------|:------------------|:----------|\n| `ProDAG`   | 25 ± 1.5    | 49 ± 3.7   | 41 ± 2.1          | 81 ± 2.8  |\n| `DAGMA`    | 56 ± 8.1    | 55 ± 7.9   | 32 ± 3.9          | 62 ± 2.0  |\n| `BayesDAG` | 34 ± 0.3    | 67 ± 1.6   | 14 ± 0.5          | 67 ± 2.0  |\n\n\n**Table R4** Larger Gumbel noise experiment with $p=100$.\n\n|            | Brier score | Exp. SHD   | Exp. F1 score (%) | AUROC (%) |\n|:-----------|:------------|:-----------|:------------------|:----------|\n| `ProDAG`   | 81 ± 4.3    | 115 ± 4.6  | 63 ± 1.4          | 92 ± 0.8  |\n| `DAGMA`    | 131 ± 4.0   | 116 ± 3.3  | 72 ± 0.8          | 91 ± 0.4  |\n| `DiBS`     | 108 ± 1.6   | 204 ± 2.0  | 38 ± 0.5          | 97 ± 0.3  |\n| `DiBS+`    | 221 ± 4.2   | 196 ± 3.3  | 41 ± 1.0          | 70 ± 0.8  |\n| `BayesDAG` | 198 ± 0.2   | 282 ± 0.9  | 2 ± 0.1           | 55 ± 0.7  |\n\n**Table R5** Larger exponential noise experiment with $p=100$.\n\n|            | Brier score | Exp. SHD   | Exp. F1 score (%) | AUROC (%) |\n|:-----------|:------------|:-----------|:------------------|:----------|\n| `ProDAG`   | 93 ± 3.4    | 123 ± 3.2  | 55 ± 1.7          | 87 ± 0.9  |\n| `DAGMA`    | 215 ± 6.6   | 197 ± 6.1  | 61 ± 0.8          | 91 ± 0.3  |\n| `DiBS`     | 129 ± 0.8   | 182 ± 0.9  | 34 ± 0.3          | 93 ± 0.3  |\n| `DiBS+`    | 188 ± 2.4   | 172 ± 2.1  | 40 ± 0.6          | 68 ± 1.1  |\n| `BayesDAG` | 199 ± 0.1   | 286 ± 1.2  | 2 ± 0.0           | 54 ± 0.5  |\n\nThe results align with our earlier findings, confirming that `ProDAG` scales to larger, more complex settings. Note that `DiBS`/`DiBS+` ran out of memory in the big nonlinear experiment so are not included. Likewise, `BOSS` is not reported in Table R3 because it does not accommodate nonlinear models, and Tables R4 and R5 because it encountered errors. \n\n> Time comparison with other Bayesian methods is not provided.\n\nWe added a time comparison in **Table R6** for linear and nonlinear DAGs (corresponding to the experiments in Figures 3 and 5 of the paper). These timings show that `ProDAG` is faster than `BayesDAG` and slower than `DiBS`. However, `ProDAG`'s superior performance in accuracy and uncertainty quantification makes it a strong choice, even with the higher computational cost compared to `DiBS`.\n\n**Table R6** Run times (seconds) against Bayesian baselines.\n\n|              | Linear | Nonlinear |\n|:-------------|:-------|:----------|\n| `ProDAG`     | 30 ± 7 | 61 ± 4    |\n| `DiBS/DiBS+` | 8 ± 0  | 40 ± 0    |\n| `BayesDAG`   | 57 ± 1 | 122 ± 5   |\n\n## Questions\n\n> Do you assume knowledge of the true noise distribution for experiments on non-Gaussian data?\n\nNo, the model is always fit with a Gaussian likelihood $x=W^\\top x+\\epsilon,\\epsilon\\sim N(0,\\Sigma)$ independent of the true noise distribution in the simulated data.\n\n> Is ProDAG robust if the log-likelihood function is misspecified? For example, using a Gaussian likelihood for data generated by non-Gaussian noise distributions, or an equal variance likelihood for a heterogeneous non-Gaussian noise distribution?\n\nYes, as in our non-Gaussian experiments, when we deliberately apply the Gaussian likelihood to data generated with Gumbel or exponential noise, `ProDAG`'s performance is roughly similar to the Gaussian noise case, indicating it is resilient to moderate likelihood misspecification.\n\n> How sensitive are the results to the choice of $\\mathbb{P}$?\n\nAny continuous prior on $\\tilde{W}$ is pushed through the same projection $\\tilde{W}\\mapsto W$, so the induced prior is largely governed by the $\\ell_1$ ball and acyclicity constraints, not by the exact density of $\\mathbb{P}$. We therefore use an isotropic Gaussian, which is quite general. Having said that, some priors such as Laplace, which shrink towards zero, are not well-suited as they can combine with the $\\ell_1$-ball constraint to induce double shrinkage."}]

Author rebuttal:

Meta-review: This paper addresses the problem of learning a distribution over directed acyclic graphs by reformulating the discrete posterior over graphs into a relaxed distribution with a projection operator, which is in turn approximated with a variational posterior. Reviewers agreed this work is of high quality and makes an important contribution. I ask the authors to make sure to honor their commitment to include further experiments with heteroscedastic noise in the final version.
Decision: Accept (poster)