Paper deep dive
SVI-DAG: A Structured Variational Inference Approach to Bayesian Causal Discovery
Shrenik Zinage
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Bayesian causal discovery seeks to determine the posterior distribution of causal theories, which are interpreted as directed acyclic graphs (DAGs) that explain the observed data. The resulting posterior allows systematic reasoning regarding epistemic uncertainty within these theories. Nonetheless, finding such graphs is difficult due to identifiability problems and limited observational data. Furthermore, precisely approximating posterior over graphs is challenging given vast range of potential DAGs. Recent Bayesian approaches have addressed some of these challenges, yet they remain limited as they fail to encode dependencies between edges, and lack principled ways to incorporate domain knowledge as inductive biases during the search process. To overcome these limitations, we propose SVI-DAG, a structured variational inference approach to Bayesian causal discovery using observational data and prior beliefs that uses normalizing flows to model dependencies between edges, supporting expressive and multimodal posterior learning over DAGs. To mitigate mode seeking behaviour in evidence lower bound optimization and promote mode coverage, we use stein variational gradient descent to update the node potentials using a kernel in acyclicity space. We evaluate SVI-DAG against 5 state-of-the-art Bayesian DAG learning methods and demonstrate superior performance in uncertainty quantification while remaining competitive in terms of structural accuracy.
Tags
Links
- Source: https://arxiv.org/abs/2608.04930v1
- Canonical: https://arxiv.org/abs/2608.04930v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
80,715 characters extracted from source content.
Expand or collapse full text
SVI-DAG: A Structured Variational Inference Approach to Bayesian Causal Discovery Shrenik Zinage ∗ Massachusetts Institute of Technology shrenik@mit.edu Abstract Bayesian causal discovery seeks to determine the posterior distribution of causal theories, which are interpreted as directed acyclic graphs (DAGs) that explain the observed data. The resulting posterior allows systematic reasoning regarding epistemic uncertainty within these theories. Nonetheless, finding such graphs is difficult due to identifiability problems and limited observational data. Furthermore, precisely approximating posterior over graphs is challenging given vast range of potential DAGs. Recent Bayesian approaches have addressed some of these challenges, yet they remain limited as they fail to encode dependencies between edges, and lack principled ways to incorporate domain knowledge as inductive biases during the search process. To overcome these limitations, we propose SVI-DAG, a structured variational inference approach to Bayesian causal discovery using observational data and prior beliefs that uses normalizing flows to model dependencies between edges, supporting expressive and multimodal posterior learning over DAGs. To mitigate mode seeking behaviour in evidence lower bound optimization and promote mode coverage, we use stein variational gradient descent to update the node potentials using a kernel in acyclicity space. We evaluateSVI-DAGagainst 5 state-of-the-art Bayesian DAG learning methods and demonstrate superior performance in uncertainty quantification while remaining competitive in terms of structural accuracy. The code is available athttps: //github.com/shrenikvz/SVI-DAG. 1 Introduction Understanding causal relationships from data (Sch ̈ olkopf et al., 2021) remains a central problem across science and engineering. Causal discovery—the task of inferring which variables directly influence others—allows researchers to move beyond mere associations and reason about the consequences of interventions in complex systems. Such a capability holds the potential to accelerate progress across a broad range of domains, including epidemiology (Tennant et al., 2021), climate science (Runge et al., 2019), economics (Imbens, 2020), psychology (Foster, 2010), and beyond. Causal relationships are commonly represented as directed acyclic graphs (DAGs) (Peters et al., 2017), where nodes correspond to variables and directed edges encode direct causal influences. Learning the DAG from data, however, constitutes a formidable combinatorial problem (Heinze-Deml et al., 2018), as the number of admissible DAGs grows super-exponentially with the number of variables and with purely observational data, the true DAG may be identifiable only up to its Markov equivalence class (MEC) 2 . This fundamental identifiability limitation motivates a Bayesian treatment (Friedman and Koller, 2003), wherein the goal is to approximate the posterior distribution over all DAGs that explain the data. Recent work on safe and trustworthy AI (Bengio et al., 2026) further argues that epistemic ∗ This work was conducted while the author was a PhD student at Purdue University, independently of the author’s dissertation research. 2 The Markov equivalence class is the set of DAGs that encode the same conditional independencies. Preprint. arXiv:2608.04930v1 [cs.LG] 5 Aug 2026 uncertainty over causal theories is essential for building non-agentic Scientist AI (Bengio et al., 2025) systems that provide reliable predictions and can serve as guardrails against overconfident decision making. Existing Bayesian approaches to DAG learning have made substantial progress, yet several limitations remain. A predominant assumption is that edges are modeled as independent random variables, thereby disregarding the rich dependencies that exist between edges. Failing to capture these dependencies leads to difficulty in capturing the multimodal structure of posterior. Moreover, existing Bayesian methods largely lack principled mechanisms for incorporating domain knowledge, which limits their capacity to guide the search toward structurally plausible regions of hypothesis space. Our Approach. To address these shortcomings, we introduceSVI-DAG, a structured variational inference (SVI) approach to Bayesian causal discovery. Our core contributions are as follows: •We propose a differentiable Bayesian approach to causal discovery that encodes dependen- cies between edge logits using conditional normalizing flows to support expressive and multimodal learning over DAGs. • Our approach provides a principled way to incorporate domain knowledge as an inductive bias using a new type of prior that adapts based on how strong or weak the prior beliefs are using a Beta-Bernoulli distribution. We further show that this construction induces a tractable Logistic-Beta prior over the edge logits. •To mitigate mode seeking behavior in evidence lower bound (ELBO) optimization, we use stein variational gradient descent (SVGD) to update node potentials using kernel in acyclicity space. This improves mode coverage and quantification of epistemic uncertainty. 2 Related Works Given the vast literature on causal discovery (Squires and Uhler, 2023), we focus exclusively on approaches that pose causal discovery as an optimization problem (Vowels et al., 2022). Causal Discovery. A foundational shift in causal discovery began with NOTEARS (Zheng et al., 2018), which reformulated combinatorial DAG search as a continuous optimization problem. This idea was extended to nonlinear settings using graph neural networks (Yu et al., 2019), gradient based neural DAG learning (Lachapelle et al., 2019), and graph autoencoder architectures (Ng et al., 2019). Further improvements include DAGs with No Curl (Yu et al., 2021) which improves efficiency via a curl free condition, and methods using acyclicity such as DAGMA (Bello et al., 2022). Complementary approaches introduced masked gradients (Ng et al., 2022) or analyzed interaction of sparsity and acyclicity penalties (Ng et al., 2020) while NoDAGs-Flow (Sethuraman et al., 2023) relaxed the acyclicity assumption entirely, learning nonlinear cyclic causal structures via normalizing flows. On combinatorial side, BOSS (Andrews et al., 2023) demonstrated that efficient order score search can rival continuous relaxation methods in speed and accuracy, while Ban et al. (2024) proposed differentiable structure learning over partial topological orderings that flexibly interpolates between unconstrained and fully ordered search. Bridging constraint based and score based paradigms, Zhou et al. (2025) reformulated conditional independence testing differentiably, allowing end-to-end gradient based learning with statistical guarantees of constraint based methods. These methods demonstrate strong empirical performance but typically return point estimates and do not quantify uncertainty over graph structures. Bayesian Causal Discovery. For Bayesian approaches, early differentiable methods include D- VAE (Zhang et al., 2019), which encodes DAGs in a latent space via variational autoencoders. VCN (Annadani et al., 2021) and BCD Nets (Cundy et al., 2021) introduced variational inference (VI) formulations over graph structures, while DiBS (Lorch et al., 2021) proposed a differentiable particle based framework. Concurrent lines of work explored differentiable DAG sampling (Charpentier et al., 2022), tractable marginal uncertainty computation (Wang et al., 2022), and amortized inference for both static (Lorch et al., 2022) and temporal (L ̈ owe et al., 2022) settings. GFlowNets emerged as a powerful alternative, allowing diverse posterior sampling over structures (Deleu et al., 2022) and joint inference over graphs and parameters (Nishikawa-Toomey et al., 2022; Deleu et al., 2023). More recently, BayesDAG (Annadani et al., 2023) introduced gradient based posterior inference with 2 improved scalability, ProDAG (Thompson et al., 2024) proposed projected VI to enforce acyclicity, and meta learning formulations (Dhir et al., 2024) have been explored to generalize across tasks. 3 Methodology 3.1 Preliminaries and Notations 3.1.1 Causal Graph and Structural Causal Model LetG = (V,E)be a DAG on the variable setV = X 1 ,...,X m whereX i represents a random variable. LetA∈0, 1 m×m be the adjacency matrix encoding the causal relationships among the mvariables with entriesA ij = 1if there is a directed edgeX j → X i inGand 0 otherwise. Let X i ⊆ Rdenote the sample space for nodeX i , and define the joint sample spaceX =X 1 ×·×X m . LetD =x (n) N n=1 be the dataset, where each observation is a tuplex (n) = (x (n) 1 ,...,x (n) m )∈X withx (n) i ∈X i . A structural causal model (SCM) states that each variable is generated as a function of its parents plus some exogenous noise. We denote byX Pa(i) the set of parental nodes ofX i so that there is an edge from X j ∈ X Pa(i) to X i in DAGG. Formally, for each node i∈1,...,m, X i = f i X Pa(i) ;θ i + ε i , wheref i (·;θ i )is a deterministic function with parametersθ i and andε i ’s are independent noise variables with strictly positive densities w.r.t Lebesgue measure. If noise variables are Gaussian and functionsf i are not linear or constant, then additive noise model (ANM) is structurally identifiable (Peters et al., 2014; Hoyer et al., 2008). Causality typically assumes structural assignments do not form cycles and they induce a DAG (Pearl, 2009). 3.1.2 Bayesian Structure Learning Our goal is to approximate the posterior distribution over all possible DAGsAand associated parametersθ = θ 1 ,...,θ m in a neural network based SCM (NCM). Specifically, we target the posterior over structures and parametersp(A,θ |D) ∝ p(D | A,θ)p(θ | A)p(A)with a prior over parameters p(θ | A) and graphs p(A) (Friedman and Koller, 2003). The likelihood is: p(D | A,θ) = N Y n=1 m Y i=1 p x (n) i | x (n) Pa(i) ;θ i , where the dependency onAis implicit in the parent index setsPa(i). Concretely, assuming the noise is Gaussian,p(X i | X Pa(i) ,θ i ) = N X i | f i (X Pa(i) ;θ i ),σ 2 i .Eachθ i determines conditional distribution of X i given its parents respecting the DAG structure A. 3.2 Our Approach In this section, we detail each component of theSVI-DAGalgorithm which is specifically designed to address some of the critical challenges in causal structure learning. 3.2.1 Domain Informed Prior over Adjacency Matrix We start by defining a prior over the DAGs. A common choice isB ij ∼ Bernoulli(p ij ).The parametersp ij encode prior beliefs about edge presence. To improve robustness to misspecified prior beliefs, we introduce a hierarchical prior overB. For each ordered pair(i,j)withi̸= j, we introduceπ ij ∈ (0, 1)such thatp(B ij | π ij ) = π B ij ij (1−π ij ) 1−B ij . Instead of fixingπ ij to the prior valuep ij , we place a hyper prior onπ ij using a Beta distributionπ ij ∼ Beta α ij ,β ij . We choose hyperparameters so that Beta-density mode is p ij given by α ij = ν ij p ij + 1, β ij = ν ij (1− p ij ) + 1, ν ij = ν min + κ|p ij − 0.5| η for i̸= j, whereν ij > 0is a concentration parameter that is strictly increasing in|p ij − 0.5|whenκ > 0and ν min ,κ, andηare constants. This choice assigns higher concentration to nominal beliefs farther from 0.5 (see Figure 1). Asν ij → ∞, the marginal Beta–Bernoulli prior over each free edgeB ij 3 converges to aBernoulli(p ij )prior (see Proposition B.1 in appendix B). LetΠ = π ij i̸=j . We sample edges as B ij ∼ Bernoulli(π ij ), which gives the joint prior as p(B, Π)∝ Y i̸=j Bernoulli(B ij | π ij ) Beta π ij | ν ij p ij + 1, ν ij (1− p ij ) + 1 . 00.10.20.30.40.50.60.70.80.91 0 5 10 15 20 p Density ν = ν min + κ|p− 0.5| η , (ν min = 1,κ = 500,η = 2) p = 0.05, ν = 102.2p = 0.1, ν = 81.0 p = 0.3, ν = 21.0p = 0.5, ν = 1.0 p = 0.7, ν = 21.0p = 0.9, ν = 81.0 p = 0.95, ν = 102.2 Figure 1: Adaptive domain informed prior over adjacency matrix Marginalizing overπ ij recovers a Beta-Bernoulli model overB ij . Since causal graphs are acyclic, we enforce acyclicity in the prior by construction (Annadani et al., 2023). We introduce a binary free edge matrixB ∈ 0, 1 m×m with zero diagonal and a vector of node poten- tialsr = (r 1 ,...,r m ) T ∈ R m . Let’s define acyclic mask and induced adjacency by M (r) = P (r)LP (r) T , A = B⊙ M (r), where⊙denotes the Hadamard product. For any(B,r), Ais a DAG (see Theorem B.3 in appendix). Moreover, for any DAGA, there always exists a corresponding pair (B,r)such thatA = B⊙ M (r)(Annadani et al., 2023; Yu et al., 2021). We therefore place joint prior on the free edges and potentialsp(B, Π,r) ∝ p(B, Π)N (r | 0,σ 2 r I m ). This construction enforces acyclicity, so no additional acyclicity indicator or penalty is required. The permutationP (r)is well defined when entries ofrare pairwise distinct. Sincep(r)is Gaussian, this condition holds almost surely. By analytically marginalizing over Π (see Derivation A.1 in appendix A), we have p(B,r) = Y i̸=j B(B ij + α ij , β ij + 1− B ij ) B(α ij ,β ij ) N r | 0,σ 2 r I m , whereBdenotes the Beta function. Under this formulation,p ij is the nominal mode of the latent edge probabilityπ ij , not the exact finiteν ij marginal probability ofB ij . The latter isα ij /(α ij +β ij ) which approaches p ij as ν ij →∞. 3.2.2 Prior over NCM Parameters conditioned on the Adjacency Matrix Next, conditional onA, we place a prior on NCM parametersθ =θ 1 ,...,θ m . We consider a prior p(θ | A) = N θ | μ p (A), diag σ 2 p (A) = d θ Y k=1 N θ k | μ p,k (A), σ 2 p,k (A) , whered θ is the dimensionality andμ p (A)∈ R d θ andσ p (A)∈ (0,∞) d θ are deterministic functions ofA. In the simplest case, these functions are constant (for example, zero mean and fixed variance), but the notation permits dependence on A. The joint prior factorizes as p(A,θ) = p(θ | A) p(A). 3.2.3 ELBO without Normalizing Flow By marginalizingΠin the prior to getp(B,r)(see Derivation A.1 in appendix A) and rearranging terms (see Derivation A.2 in appendix A), we have the ELBO as ELBO(φ) = E q φ (B,r) h E q φ (θ|A) [logp(D | A,θ)] | z expected log likelihood −D KL (q φ (θ | A)∥p(θ | A)) | z KL divergence for θ i − D KL (q φ (B,r)∥p(B,r)) |z KL divergence for B and r . where q φ (θ | A) and q φ (B,r) denotes the conditional and structural guides. Remark 3.1. Note that the above ELBO is not equivalent to the ELBO defined usingq φ (A)under deterministic mapA = B⊙ M (r)because the pushforward from(B,r)toAis many-to-one. If q φ (A)andp(A)are defined as the pushforwards ofq φ (B,r)andp(B,r), then Proposition B.2 impliesD KL q φ (B,r)∥p(B,r) ≥ D KL q φ (A)∥p(A) . Consequently, an ELBO written directly in terms ofAis at least as tight as the ELBO in(B,r)space. However, the induced priorp(A)is generally intractable, whereasp(B,r)admits a simple factorized form. We therefore work with Eq. A.1 for computational tractability, even though ELBO in A space would be tighter. 4 3.3 Soft Relaxation for Adjacency Matrix and Acyclicity SinceBis discrete, we cannot differentiate with respect toB ij . We therefore introduce a continuous relaxation forBvia the Gumbel-Softmax trick (Jang et al., 2016; Maddison et al., 2016; Ng et al., 2022). For each(i,j)withi ̸= j, we introduce a learnable edge logitγ ij ∈ Rand a temperature T > 0. We draw logistic noiseR ij iid ∼ Logistic(0, 1)implemented asR ij = logU ij − log(1− U ij ), U ij iid ∼ Uniform(0, 1). Equivalently,R ij can be generated as the difference of two i.i.d. Gumbel(0, 1)variables. We use the same noise to define both a hard Bernoulli sample and its differentiable relaxation: B ij = 1γ ij + R ij > 0, ̃ B ij = σ (γ ij + R ij )/T , σ(x) = 1/(1 + e −x ),(1) whereB i = ̃ B i = 0. ThenB ij ∼ Bernoulli(σ(γ ij )), while ̃ B ij ∈ (0, 1)converges pointwise toB ij asT → 0, except on the probability-zero eventγ ij + R ij = 0(Ng et al., 2022). For the order, letP (r)be the hard permutation obtained by sortingr. We also define a backward surrogate through the Sinkhorn operatorS(Adams and Zemel, 2011). Forτ > 0, letS 0 (r) = r o T , whereo = (m,m− 1,..., 1) T , and applyingK S alternating row and column normalizations to exp(S 0 (r)/τ ) we have P τ,K S (r) =S K S (exp(S 0 (r)/τ )), M τ,K S (r) = P τ,K S (r)LP τ,K S (r) T . Double stochasticity is guaranteed in theK S →∞limit. More importantly, even an exactly doubly stochastic nonpermutation matrix generally makesM τ,K S (r) dense, so this relaxed order mask is not generally acyclic and can have positive entries in both directions. Lemma B.4 (see appendix B) states thatP τ,K S (r) → P (r) andM τ,K S (r) → M (r) asτ → 0andK S → ∞. We distinguish the hard DAG from its soft backward surrogate as A = B⊙ M (r), ̃ A = ̃ B⊙ M τ,K S (r)⊙ (1 m×m − I m ), A ST = ̃ A + sg(A− ̃ A),(2) whereM (r) = P (r)LP (r) T andsgis stop-gradient. ThusA ST = Ain the forward pass and every evaluated graph is a binary DAG by Theorem B.3. The diagonal factor removes self-loops from ̃ A, but does not make the finite-temperature relaxation acyclic. 3.4 Forward Model −6−4−20246 0 0.2 0.4 0.6 0.8 1 1.2 γ p ( γ ) ν = ν min + κ|p− 0.5| η , (ν min = 1, κ = 500, η = 2) p = 0.05, ν = 102.2p = 0.1, ν = 81.0 p = 0.3, ν = 21.0p = 0.5, ν = 1.0 p = 0.7, ν = 21.0p = 0.9, ν = 81.0 p = 0.95, ν = 102.2 logit(p) Figure 2: Induced prior over edge logits We now define forward model using Bayesian MLP for nodeithat uses the parentX Pa(i) as input and outputs a prediction ˆ X i . Formally, ˆ X i = f i X Pa(i) ; θ i .To incorporate adjacency, letX = (X 1 ,...,X m ) T and let (A ST ) i denote thei th row of the straight-through matrix in Eq. 2. We define g i (x;θ i ) = f i (A ST ) i ⊙ x;θ i . SinceA ST = Ain the forward pass, the evaluated input zeros every nonparent exactly and all likelihood values cor- respond to a hard DAG. During backpropagation, deriva- tives are routed through the generally cyclic relaxation ̃ A. These straight through derivatives are biased surrogate gradients. 3.5 ELBO with Normalizing Flow In many applications, the assumption thatq φ (B)factorizes independently across edges can be restrictive. SVI (Rezende and Mohamed, 2015) addresses this limitation by constructing a more expressiveq φ (B). We define a standard base distributionp 0 (z)for edge latent variables. We introduce an invertible conditional mapT φ : R m(m−1) × R m → R m(m−1) , parameterized byφ. The transformation givesγas a function ofzandrgiven byγ = T φ (z;r). By rearrangement (see Derivation A.3 in appendix A), we have the revised ELBO as ELBO(φ,ψ) = E r∼q ψ (r) E z∼p 0 E B∼p(·|T φ (z;r)) h E q φ (θ|A) logp(D | A,θ) |z expected log likelihood −D KL q φ (θ | A)∥p(θ | A) |z KL divergence for θ + log det ∂T φ (z;r) ∂z − logp 0 (z) + logp T φ (z;r) i − D KL q ψ (r)∥p(r) | z KL divergence for potentials . (3) 5 Under Derivation A.4 (see appendix A) and the modeling choice that model and guide share the hard conditionalp(B | γ), the domain informed prior now enters through prior over edge logitsp(γ)(see Figure 2) as p(γ) = Y i̸=j 1 B(α ij ,β ij ) σ(γ ij ) α ij 1− σ(γ ij ) β ij , σ(x) = 1/(1 + e −x ). Remark 3.2. Had we instead placed a Bernoulli priorB ij ∼ Bernoulli(p ij )on each edge, the edge probability would be deterministic (i.e.,π ij = p ij ). Equivalently,p(π ij ) = δ(π ij −p ij ), and the logit transformation would inducep(γ ij ) = δ(γ ij − logit(p ij )).Hence, a fixed Bernoulli prior induces a degenerate prior overγrather than an absolutely continuous density. This is not suitable to our formulation, sinceq φ (γ | r)is continuous—so the KL divergenceD KL (q φ (γ | r)∥p(γ))to this degenerate prior diverges unless the flow collapses to the same point mass. The Beta prior over π ij avoids this issue and gives the Logistic-Beta density. We further note that although the prior overπ ij is parameterized to have modep ij , this mode is not preserved under the nonlinear logit transformation. In particular,logit(p ij )is the image of the mode inπ ij space, whereas the mode of the Logistic-Beta density overγ ij islog(α ij /β ij ). Consequently, the peaks ofp(γ)need not align exactly with logit(p ij ), as seen in Figure 2. A single sample Monte Carlo (MC) estimator of revised ELBO is unbiased and practical in high dimension. The straight through backward pass is generally a biased gradient estimator. If both q φ (θ | A)andp(θ | A)are diagonal Gaussians with strictly positive component scales, their KL divergence has the closed form D KL q φ (θ | A)∥p(θ | A) = 1 2 P d θ k=1 " σ 2 φ,k (A) σ 2 p,k (A) + μ φ,k (A)−μ p,k (A) 2 σ 2 p,k (A) − 1 + 2 log σ p,k (A) σ φ,k (A) # .(4) 3.6Updating Node Potentials with SVGD using a Kernel defined in Relaxed Acyclicity Space Note that, even with a dependent guide, standard VI maximizes the ELBO, which is equivalent to minimizing the reverse KL divergenceD KL (q φ ∥p). For multimodal posteriors, this reverse KL objective can be mode-seeking (Blei et al., 2017; Zhang et al., 2018). In particular, it encourages the approximation to concentrate on a single mode, largely to avoid allocating mass to low-density posterior regions. Consequently, although a mixture-of-Gaussians guide forrcould, in principle, represent multiple modes, standard ELBO optimization would still not guarantee recovery of all posterior modes. To promote mode coverage, we replace the parametric guide overrwith a nonpara- metric particle approximation and update the particles using SVGD (Liu and Wang, 2016). If we condition the revised ELBO (see Equation 3) on a fixed r, we have L(φ;r) = E z∼p 0 E B∼p(·|T φ (z;r)) h E q φ (θ|A) logp(D | A,θ) |z expected log likelihood −D KL q φ (θ | A)∥p(θ | A) |z KL divergence for θ + log det ∂T φ (z;r) ∂z − logp 0 (z) + logp T φ (z;r) i . (5) This constitutes a genuine lower bound on the conditional log evidenceL(φ;r) ≤ logp(D | r). More precisely, defining∆ φ (r) = D KL (q φ (γ,B,θ | r)∥p(γ,B,θ |D,r)), the standard variational identity giveslogp(D | r)−L(φ;r) = ∆ φ (r)≥ 0. To introduce repulsive diversity pressure over node orders without assigning a density to a finite set of particles, we construct a proper target directly from the fixed r conditional bound. Specifically, for a fixed φ, we define Z φ = Z p(r) expL(φ;r)dr, ̃p φ (r |D) = Z −1 φ p(r) expL(φ;r). If0 < p(D) <∞andr 7→L(φ;r)is measurable, the nonnegativity∆ φ (r)≥ 0implies0 < Z φ ≤ p(D). Hence, ̃p φ (r | D)constitutes a proper, lower-bound-induced target that relates to the exact potential posterior by ̃p φ (r |D) = p(r |D) exp−∆ φ (r) E p(r|D) [exp−∆ φ (r)] . Thus, ̃p φ (r | D) = p(r | D)holds if and only if∆ φ (r)is constant forp(r | D)almost everyr including, in particular, the special case of an exact conditional guide. Assume thatr 7→L(φ;r)is 6 differentiable on the interior of each order chamber. This requires the flow, likelihood, and parameter guide to be differentiable there and an integrable function that allows differentiation under the continuous expectations. The derivative of a sampled hard threshold is not this exact derivative. It would require an analytic finite sum or an unbiased discrete score estimator. Under these assumptions, away from tie hyperplanes, the score is∇ r log ̃p φ (r | D) = ∇ r [L(φ;r) + logp(r)]. The hard permutation can makeLand ̃p φ (r | D)discontinuous across tie hyperplanes, even though these hyperplanes carry prior measure zero. We maintain a finite setr (k) K part k=1 purely as optimization variables. Holding these particles fixed, the exact conditional objective average that we optimize overφisK −1 part P K part k=1 L(φ;r (k) ) . Notably, this is not a global ELBO for an empirical guide over r. Let’s define the continuously differentiable, positive-definite Gaussian kernelk M (ψ,ψ ′ )and letψ(r) = vec(M τ,K S (r)). Its pullbackk r (r,r ′ ) = k M (ψ(r),ψ(r ′ ))is, in turn, a continuously differentiable, positive-definite kernel on potential space. Applying the standard smooth-target SVGD formula to this score gives the particle update rule r (k) ← r (k) + η r 1 K part K part X j=1 h k r (r (j) ,r (k) )∇ r (j) log ̃p φ (r |D) +∇ r (j) k r (r (j) ,r (k) ) i .(6) Here,∇ r (j) k r differentiates throughψ(r (j) ). The second term in the bracket is the repulsive force that helps in mode coverage. 3.7 Reparameterization Trick for Estimating Gradients In order to calculate gradients, letf i be anL-layer MLP with parametersθthat contain weights and biasesW (l) ,b (l) L l=1 such thatW (l) ∈ R d (l) out ×d (l) in andb (l) ∈ R d (l) out . We approximate the conditional posterior over network parameters with a Gaussian guide. In software the hypernetworkH φ receives A ST since its forward value is exactlyA, this evaluates the same density asq φ (θ | A)while routing a surrogate derivative through ̃ A: q φ (θ | A ST ) =N θ | μ φ (A ST ), diag σ 2 φ (A ST ) = d θ Y k=1 N θ k | μ φ,k (A ST ), σ 2 φ,k (A ST ) , where(μ φ ,ρ φ ) = H φ (A ST )andσ φ = softplus(ρ φ ) = log(1 + exp(ρ φ )), which enforces positivity. The conditioning is therefore explicit and differentiable. Changing A ST changes the full variational posterior over θ through H φ . θ = μ φ (A ST ) + σ φ (A ST )⊙ ε, ε∼N (0,I d θ ), which preserves differentiability through bothH φ and the predictorsf i . ̃ B is reparameterized using the Gumbel-Sigmoid relaxation as shown in Eq. 1. We refer the reader to Appendix A.6 for the pseudo code and computational complexity of SVI-DAG algorithm. 4 Results and Discussion We evaluateSVI-DAGagainst five state-of-the-art (SOTA) Bayesian DAG learning methods on linear synthetic data, nonlinear synthetic data and real data using brier score, expected structural hamming distance (SHD), expected F1 score, and area under the receiver operating characteristic curve (AUROC). The rationale underlying the selection of these metrics is provided in appendix C. Datasets.To isolate effect of domain informed prior, we use a two node graph with data generated from ANM with Gaussian noise, where small graph size makes role of prior tractable to analyze. For synthetic experiments at scale, we generate linear and nonlinear structural equation models over Erd ̋ os–R ́ enyi (ER) DAGs with 25 and 50 nodes. For real data, we evaluate on Sachs dataset. Baselines. We focus our evaluation on Bayesian methods to DAG learning and compare with ProDAG(Thompson et al., 2024),BayesDAG(Annadani et al., 2023),DDS(Charpentier et al., 2022), BCDNets (Cundy et al., 2021), andDIBS(Lorch et al., 2021). Appendix E describes the hyperparam- eters used for benchmarking these algorithms. 7 4.1 Effect of Domain Informed Prior Linear (DAG not identifiable from data) Nonlinear (DAG identifiable from data) Incorr. Noninf. Corr. 84.715.3 40.142.517.4 83.216.8 13.386.4 80.314.5 95.2 0255075100 0255075100 Posterior mass (%) X 1 → X 2 (true) X 1 ← X 2 X 1 X 2 Incorr. Noninf. Corr. Figure 3: Effect of prior on 2 node graph Figure 3 illustrates effect of prior on a 2 node graph using 1000 samples, chosen to balance influence of likelihood and prior. Af- ter optimization, we draw10000hard posterior DAG samples A = B⊙ M (r),rsampled uniformly from the SVGD particles,γ fromq(γ | r), hard Bernoulli edgesB = 1γ + R > 0and the hard permutation maskM (r) = P (r)LP (r) ⊤ , so every sample is a binary DAG and report the proportion of samples assigned to x 1 → x 2 ,x 1 ← x 2 , andx 1 ⊥ x 2 . With ground truthx 1 → x 2 , linear Gaussian setting (non-identifiable) shows that prior dominates orientation: an incorrect prior (p = 0.01 correct, 0.99 incorrect) con- centrates mass on the wrong direction (84.7%), a noninformative prior (p = 0.5) gives balanced posterior, and correct prior (p = 0.99 correct, 0.01 incorrect) shifts mass to true direction (83.2%). In nonlinear Gaussian setting (identifiable), prior acts primarily as regularizer: the incorrect prior suppresses both edges, concentrating 86.4% onx 1 ⊥ x 2 , while a correct prior sharpens posterior and assigns highest mass to x 1 → x 2 (95.2%). 4.2 Linear Synthetic Data 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0.04 0.06 0.08 0.10 Brier score 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 40 60 80 100 Expected SHD 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 20 40 60 80 Expected F 1 score (%) 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 70 80 90 100 AUROC (%) SVI-DAGProDAGBayesDAGDiBSDDSBCD Nets Figure 4: Results on synthetic data generated from linear ER DAGs withp = 25nodes ands = 40 edges under Gaussian noise. Each point denotes sample mean, and accompanying error bars indicate standard error, both computed across 5 independent draws of the dataset. 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0.02 0.04 0.06 Brier score 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 100 200 300 Expected SHD 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 20 40 60 80 Expected F 1 score (%) 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 50 60 70 80 90 100 AUROC (%) SVI-DAGProDAGBayesDAGDiBSDDSBCD Nets Figure 5: Results on synthetic data generated from linear ER DAGs withp = 50nodes ands = 80 edges under Gaussian noise. Each point denotes sample mean, and accompanying error bars indicate standard error, both computed across 5 independent draws of the dataset. Figures 4 and 5 presents results on datasets generated from linear ER DAGs withp = 25,s = 40and p = 50,s = 80across increasing sample sizes. Since the underlying structure is identifiable only up to its MEC, all metrics are evaluated at CPDAG level.SVI-DAGachieves lower Brier scores and higher AUROC across most of the sample sizes, indicating improved uncertainty quantification and better separation, while remaining competitive with baselines in terms of E-SHD, and E-F1 score. 4.3 Nonlinear Synthetic Data Figures 6 and 7 presents results on datasets generated from nonlinear ER DAGs withp = 25,s = 40 andp = 50,s = 80across increasing sample sizes. In this setting, the underlying DAG is identifiable, and all metrics are therefore evaluated at the DAG level. We can observe thatSVI-DAGattains a 8 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0.04 0.05 0.06 0.07 Brier score 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 40 60 80 100 Expected SHD 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0 20 40 Expected F 1 score (%) 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 60 70 80 90 AUROC (%) SVI-DAGProDAGBayesDAGDiBSDDSBCD Nets Figure 6: Results on synthetic data generated from nonlinear ER DAGs withp = 25nodes and s = 40edges under Gaussian noise. Each point denotes the sample mean, and the accompanying error bars indicate the standard error, both computed across 5 independent draws of the dataset. 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0.02 0.03 0.04 Brier score 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 100 200 300 400 Expected SHD 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 0 10 20 30 40 50 Expected F 1 score (%) 10 2.0 10 2.5 10 3.0 10 3.5 10 4.0 Sample size 50 60 70 80 90 AUROC (%) SVI-DAGProDAGBayesDAGDiBSDDSBCD Nets Figure 7: Results on synthetic data generated from nonlinear ER DAGs withp = 50nodes and s = 80edges under Gaussian noise. Each point denotes sample mean, and accompanying error bars indicate standard error, both computed across 5 independent draws of the dataset. lower Brier score and higher AUROC across most sample sizes, while also staying competitive with other algorithms in this setting as well. 4.4 Real Data Table 1: Performance on Sachs dataset. The average and standard errors are measured over 10-fold cross-validation splits of data. The best value of each DAG metric is indicated in bold. Brier score (↓)Exp. SHD (↓)Exp. F1 score (↑)AUROC (↑) SVI-DAG0.185± 0.00725.70± 0.9228.95± 1.7358.16± 2.27 ProDAG0.159± 0.00319.68± 0.3514.32± 1.5452.95± 1.81 BayesDAG0.174± 0.00828.66± 0.7826.75± 2.0862.23± 1.93 DDS0.168± 0.00424.62± 0.5623.21± 1.5251.35± 2.18 BCD Nets0.153± 0.00420.94± 0.5010.79± 2.1655.91± 2.07 DiBS0.180± 0.00440.74± 0.5027.78± 0.7465.63± 2.07 Table 3: Wall-clock time per train/evaluate on the Sachs dataset on an A100-40GB GPU. Time per run (s) SVI-DAG134.2± 30.9 ProDAG100.2± 4.1 BayesDAG266.5± 13.0 DDS65.1± 1.9 BCD Nets151.1± 2.5 DiBS5.6± 0.7 The flow cytometry dataset (Sachs et al., 2005) is a standard benchmark for causal discovery in biological systems, con- sisting of7466samples over11variables and18edges corre- sponding to phosphoproteins and phospholipids. Tables 1 and 2 provides the performances on Sachs dataset using DAG and CPDAG metrics. We can observe thatSVI-DAGoutperforms other algorithms in terms of F1 score for DAG metrics and AUROC for CPDAG metrics. Table 3 reports the computational time of the algorithms on this dataset. We can see comparable computational time of SVI-DAG with other algorithms. We note that all results reported for SVI-DAG in this paper correspond to a single choice of normaliz- ing flow—neural spline flows. The performance of SVI-DAG is largely governed by this choice, and our algorithm has been deliberately developed to be flow-agnostic, so that any normalizing flow can be substituted in. Hence, the comparisons presented here should not be taken too seriously. 9 Table 2: Performance on Sachs dataset. The average and standard errors are measured over 10-fold cross-validation splits of data. The best value of each CPDAG metric is indicated in bold. Brier score (↓)Exp. SHD (↓)Exp. F1 score (↑)AUROC (↑) SVI-DAG0.237± 0.01623.08± 0.8754.05± 2.0770.17± 2.28 ProDAG0.267± 0.00117.36± 0.0731.55± 0.6462.70± 0.90 BayesDAG0.247± 0.01026.38± 0.4250.10± 1.3662.88± 1.42 DDS0.240± 0.00523.02± 0.2943.46± 0.9168.91± 1.30 BCD Nets0.259± 0.00218.68± 0.1626.30± 1.6156.65± 1.62 DiBS0.215± 0.00335.49± 0.2254.46± 0.2966.80± 1.10 5 Conclusion We presentedSVI-DAG, a differentiable Bayesian approach to causal discovery using observational data and prior beliefs that uses normalizing flows to encode dependencies between edges for expres- sive and multimodal learning over DAGs. Empirically, we demonstrate thatSVI-DAGachieves better uncertainty quantifcation and competitive performance with SOTA Bayesian methods in terms of accuracy. 6 Acknowledgment The author gratefully acknowledges the computing resources provided by Purdue University, and thanks Prof. Ilias Bilionis for access to these resources and for the freedom to pursue this line of work during the author’s doctoral studies. References Ryan Prescott Adams and Richard S Zemel. Ranking via sinkhorn propagation. arXiv preprint arXiv:1106.1925, 2011. Bryan Andrews, Joseph Ramsey, Ruben Sanchez Romero, Jazmin Camchong, and Erich Kummerfeld. Fast scalable and accurate discovery of dags using the best order score search and grow shrink trees. Advances in neural information processing systems, 36: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. arXiv preprint arXiv:2106.07635, 2021. Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong. Bayesdag: Gradient-based posterior inference for causal discovery. Advances in Neural Information Processing Systems, 36:1738–1763, 2023. Taiyu Ban, Lyuzhou Chen, Xiangyu Wang, Xin Wang, Derui Lyu, and Huanhuan Chen. Differentiable structure learning with partial orders. Advances in Neural Information Processing Systems, 37: 117426–117455, 2024. Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization. Advances in Neural Information Processing Systems, 35:8226–8239, 2022. Yoshua Bengio, Michael Cohen, Damiano Fornasiere, Joumana Ghosn, Pietro Greiner, Matt Mac- Dermott, S ̈ oren Mindermann, Adam Oberman, Jesse Richardson, Oliver Richardson, et al. Su- perintelligent agents pose catastrophic risks: Can scientist ai offer a safer path? arXiv preprint arXiv:2502.15657, 2025. Yoshua Bengio, Stephen Clare, Carina Prunkl, Maksym Andriushchenko, Ben Bucknall, Malcolm Murray, Rishi Bommasani, Stephen Casper, Tom Davidson, Raymond Douglas, et al. International ai safety report 2026. arXiv preprint arXiv:2602.21012, 2026. 10 David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518):859–877, 2017. Bertrand Charpentier, Simon Kibler, and Stephan G ̈ unnemann. Differentiable DAG Sampling. International Conference on Learning Representations, 2022. Chris Cundy, Aditya Grover, and Stefano Ermon. Bcd nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095–7110, 2021. Tristan Deleu, Ant ́ onio G ́ ois, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. Bayesian structure learning with generative flow networks. In Uncertainty in Artificial Intelligence, pages 518–528. PMLR, 2022. 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. Advances in Neural Information Processing Systems, 36:31204–31231, 2023. Anish Dhir, Matthew Ashman, James Requeima, and Mark van der Wilk. A meta-learning approach to bayesian causal discovery. arXiv preprint arXiv:2412.16577, 2024. Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019. E Michael Foster. Causal inference and developmental psychology. Developmental psychology, 46 (6):1454, 2010. Nir Friedman and Daphne Koller. Being bayesian about network structure. a bayesian approach to structure discovery in bayesian networks. Machine learning, 50(1):95–125, 2003. Christina Heinze-Deml, Marloes H Maathuis, and Nicolai Meinshausen. Causal structure learning. Annual Review of Statistics and Its Application, 5:371–391, 2018. Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Sch ̈ olkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008. Guido W Imbens. Potential outcome and directed acyclic graph approaches to causality: Relevance for empirical practice in economics. Journal of Economic Literature, 58(4):1129–1179, 2020. Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. arXiv preprint arXiv:1611.01144, 2016. Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. S ́ ebastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural dag learning. arXiv preprint arXiv:1906.02226, 2019. Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. Advances in neural information processing systems, 29, 2016. Lars Lorch, Jonas Rothfuss, Bernhard Sch ̈ olkopf, and Andreas Krause. DiBS: Differentiable Bayesian Structure Learning. Advances in Neural Information Processing Systems, 2021. Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Sch ̈ olkopf. Amortized inference for causal structure learning. Advances in Neural Information Processing Systems, 35: 13104–13118, 2022. Sindy L ̈ owe, David Madras, Richard Zemel, and Max Welling. Amortized causal discovery: Learning to infer causal graphs from time-series data. In Conference on Causal Learning and Reasoning, pages 509–525. PMLR, 2022. 11 Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint arXiv:1611.00712, 2016. Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. A graph autoencoder approach to causal structure learning. arXiv preprint arXiv:1911.07420, 2019. Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags. Advances in Neural Information Processing Systems, 33:17943–17954, 2020. Ignavier Ng, Shengyu Zhu, Zhuangyan Fang, Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient-based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pages 424–432. SIAM, 2022. Mizu Nishikawa-Toomey, Tristan Deleu, Jithendaraa Subramanian, Yoshua Bengio, and Laurent Charlin. Bayesian learning of causal structure and mechanisms with gflownets and variational bayes. arXiv preprint arXiv:2211.02763, 2022. Judea Pearl. Causality. Cambridge university press, 2009. Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Sch ̈ olkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, 15(1):2009–2053, 2014. Jonas Peters, Dominik Janzing, and Bernhard Scholkopf. Elements of causal inference: foundations and learning algorithms. MIT press, 2017. Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International conference on machine learning, pages 1530–1538. PMLR, 2015. Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Mu ̃ noz-Mar ́ ı, et al. Inferring causation from time series in earth system sciences. Nature communications, 10(1):2553, 2019. Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 2005. Bernhard Sch ̈ olkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612–634, 2021. Muralikrishnna G Sethuraman, Romain Lopez, Rahul Mohan, Faramarz Fekri, Tommaso Biancalani, and Jan-Christian H ̈ utter. Nodags-flow: Nonlinear cyclic causal structure learning. In International Conference on Artificial Intelligence and Statistics, pages 6371–6387. PMLR, 2023. Chandler Squires and Caroline Uhler. Causal structure learning: A combinatorial perspective. Foundations of Computational Mathematics, 23(5):1781–1815, 2023. Peter WG Tennant, Eleanor J Murray, Kellyn F Arnold, Laurie Berrie, Matthew P Fox, Sarah C Gadd, Wendy J Harrison, Claire Keeble, Lynsie R Ranker, Johannes Textor, et al. Use of directed acyclic graphs (dags) to identify confounders in applied health research: review and recommendations. International journal of epidemiology, 50(2):620–632, 2021. Ryan Thompson, Edwin V Bonilla, and Robert Kohn. Prodag: Projected variational inference for directed acyclic graphs. arXiv preprint arXiv:2405.15167, 2024. Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden. D’ya like dags? a survey on structure learning and causal discovery. ACM Computing Surveys, 55(4):1–36, 2022. Benjie Wang, Matthew R Wicker, and Marta Kwiatkowska. Tractable uncertainty for structure learning. In International Conference on Machine Learning, pages 23131–23150. PMLR, 2022. Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: Dag structure learning with graph neural networks. In International conference on machine learning, pages 7154–7163. PMLR, 2019. 12 Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji. Dags with no curl: An efficient dag structure learning approach. In International Conference on Machine Learning, pages 12156–12166. Pmlr, 2021. Cheng Zhang, Judith B ̈ utepage, Hedvig Kjellstr ̈ om, and Stephan Mandt. Advances in variational inference. IEEE transactions on pattern analysis and machine intelligence, 41(8):2008–2026, 2018. Muhan Zhang, Shali Jiang, Zhicheng Cui, Roman Garnett, and Yixin Chen. D-vae: A variational autoencoder for directed acyclic graphs. Advances in neural information processing systems, 32, 2019. Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018. Jincheng Zhou, Mengbo Wang, Anqi He, Yumeng Zhou, Hessam Olya, Murat Kocaoglu, and Bruno Ribeiro. Differentiable constraint-based causal discovery. arXiv preprint arXiv:2510.22031, 2025. 13 Appendix This Appendix includes additional details for the paper, “SVI-DAG: A Structured Variational Infer- ence Approach to Bayesian Causal Discovery”. A Theory A.1 Derivation of p(B,r) We start from the construction that enforces acyclicity by mapping the free edge matrixB ∈ 0, 1 m×m and the potentialsr ∈ R m to the adjacencyA = B⊙ M (r)whereM (r)is the acyclic mask induced by the order defined byr. The joint prior overB, the latent edge probabilities Π =π ij i̸=j , and r is p(B, Π,r) ∝ Y i̸=j Bernoulli(B ij | π ij ) Beta(π ij | α ij ,β ij ) N r | 0,σ 2 r I m , withB i = 0for alliand with independence betweenrand the set of free edges in the prior. The factors are Bernoulli(B ij | π ij ) = π B ij ij (1− π ij ) 1−B ij ,Beta(π ij | α ij ,β ij ) = π α ij −1 ij (1− π ij ) β ij −1 B(α ij ,β ij ) , whereBdenotes the Beta function. Our goal is to integrate outπ ij to get the marginal prior on(B,r). For a fixed ordered pair (i,j) with i̸= j we compute p(B ij ) = Z 1 0 Bernoulli(B ij | π ij ) Beta(π ij | α ij ,β ij ) dπ ij . Substituting the explicit forms gives p(B ij ) = 1 B(α ij ,β ij ) Z 1 0 π B ij +α ij −1 ij (1− π ij ) (1−B ij )+β ij −1 dπ ij . The integral equals a Beta function with updated parameters, which gives the standard Beta Bernoulli marginal p(B ij ) = B(B ij + α ij , β ij + 1− B ij ) B(α ij ,β ij ) . Using independence across ordered pairs in the prior and independence ofrfromBin the prior, we get the full marginal as p(B,r) = Y i̸=j B(B ij + α ij , β ij + 1− B ij ) B(α ij ,β ij ) N r | 0,σ 2 r I m . The diagonal constraintB i = 0excludes self loops and is implemented by masking the diagonal of B. Note that instead of p(B ij = 1) = p ij and p(B ij = 0) = 1− p ij we have p(B ij = 1) = α ij α ij + β ij , p(B ij = 0) = β ij α ij + β ij , for a chosen finite ν ij . This completes the derivation. A.2 Derivation of ELBO without the normalizing flow GivenA = B⊙M (r), which is a deterministic function of(B,r)alone, we have the full joint model as p(B, Π,r,θ,D) = p(B, Π,r) p(θ | A) p(D | A,θ). The posterior is p(B, Π,r,θ |D) = p(B, Π,r,θ,D) R p(B, Π,r,θ,D)dB dΠdrdθ . 14 Since this posterior is intractable, we approximate it with a variational distribution q φ (B, Π,r,θ) = q φ (B, Π,r) q φ (θ | A), A = B⊙ M (r). Hereφdenotes the collection of variational parameters forB,Π,r, andθ. We begin by writing the evidence as logp(D) = log R p(B, Π,r,θ,D)dB dΠdrdθ. By definition, the ELBO is ELBO(φ) = E q φ (B,Π,r,θ) [logp(B, Π,r,θ,D)− logq φ (B, Π,r,θ)]. Using p(B, Π,r,θ,D) = p(D | A,θ) p(θ | A) p(B, Π,r), we have ELBO(φ) = E q φ (B,Π,r,θ) [logp(D | A,θ) + logp(θ | A) + logp(B, Π,r)− logq φ (B, Π,r,θ)]. If we assume independence between structural variables and parameter variables in the guide, we can factorize: q φ (B, Π,r,θ) = q φ (B, Π,r)q φ (θ | A), so that logq φ (B, Π,r,θ) = logq φ (B, Π,r) + logq φ (θ | A). By linearity of expectation and the definition of the KL divergence (D KL (q∥p) = E q [logq− logp]), the ELBO can then be rearranged as: ELBO(φ) = E q φ (B,Π,r) h E q φ (θ|A) [logp(D | A,θ) + logp(θ | A)] − E q φ (θ|A) [logq φ (θ | A)] + logp(B, Π,r)− logq φ (B, Π,r) i = E q φ (B,Π,r) h E q φ (θ|A) [logp(D | A,θ)] i − E q φ (B,Π,r) h D KL (q φ (θ | A)∥p(θ | A)) i − D KL (q φ (B, Π,r)∥p(B, Π,r)). If we analytically marginalize Π in the prior to get p(B,r) (see Derivation A.1) we have p(B,r) = Y i̸=j B (B ij + α ij ,β ij + 1− B ij ) B (α ij ,β ij ) N (r | 0,σ 2 r I m ), whereBrepresents the Beta function. SinceΠis an auxiliary latent variable, marginalizing it out does not change the evidencep(D). We therefore perform variational inference directly on the collapsed model p(D,B,r,θ) = p(D | A,θ)p(θ | A)p(B,r), A = B⊙ M (r), and the corresponding guide q φ (B,r)q φ (θ | A). The ELBO of this collapsed model is ELBO(φ) = E q φ (B,r) h E q φ (θ|A) [logp(D | A,θ)] |z expected log likelihood −D KL (q φ (θ | A)∥p(θ | A)) |z KL divergence for θ i − D KL (q φ (B,r)∥p(B,r)) | z KL divergence for B and r . (A.1) A.3 Derivation of ELBO with Normalizing Flow We derive the ELBO when the edge logitsγare generated by a conditional normalizing flowT φ and acyclicity is enforced by construction. The construction uses a free edge matrixBwith zero diagonal and a vector of node potentials r ∈ R m . The adjacency matrix used inside the model is given by A = B⊙ M (r), M (r) = P (r)LP (r) T , whereP (r),L, andM (r)are as defined in the main text. Note that a factorization whereq φ (A) assumes independence betweenrandγprevents the model from capturing multiple modes in the DAG posterior. For example, if two DAGs fit the data equally well but imply different topological orders, an independent guideq φ (r)q φ (γ)forces the model to average edge probabilities across incompatible orders, which can give a dense graph that does not coincide with any high posterior DAG. So to represent multiple modes, we model dependence betweenrandγ. Letz ∼ p 0 (z)and 15 γ = T φ (z;r), where for eachr,T φ (·;r)is aC 1 diffeomorphism with nonzero Jacobian determinant. Assume0 < p(D) <∞,q ψ (r)≪ p(r),q φ (γ | r)≪ p(γ),q φ (θ | A)≪ p(θ | A), and integrability underq φ,ψ of all displayed log-density terms, including the log likelihood. The induced conditional guide over γ is given by the change of variables q φ (γ | r) = p 0 T −1 φ (γ;r) det ∂T −1 φ (γ;r) ∂γ = p 0 (z) det ∂T φ (z;r) ∂z −1 . The hard DAG adjacency is a deterministic function of B and r. The joint model is p(D,B,r,θ,γ) = p(D | A,θ)p(θ | A)p(B | γ)p(r)p(γ),(A.2) wherep(B | γ)is the Bernoulli conditional andp(r)is the potential prior. Note that while the guideq φ conditionsγonr, the priorp(r)p(γ)typically factorizes unless specific domain knowledge suggests otherwise. We have the joint guide as q φ,ψ (γ,B,r,θ) = q φ (γ | r)p(B | γ)q ψ (r)q φ (θ | A),(A.3) to avoid the collapse to a single mode over DAGs. By definition, ELBO(φ,ψ) = E q φ,ψ (γ,B,r,θ) h logp(D,B,r,θ,γ)− logq φ,ψ (γ,B,r,θ) i . Using Eq. A.2 and Eq. A.3 we expand ELBO(φ,ψ) = E q ψ (r)q φ (γ|r)p(B|γ)q φ (θ|A) h logp(D | A,θ) + logp(θ | A) + logp(B | γ) + logp(r) + logp(γ)− logq φ (θ | A)− logp(B | γ) − logq ψ (r)− logq φ (γ | r) i .(A.4) We first integrate overθconditionally on the hard adjacencyA. The terms that depend onθare logp(D | A,θ) + logp(θ | A)− logq φ (θ | A). Therefore ELBO(φ,ψ) = E q ψ (r)q φ (γ|r)p(B|γ) h E q φ (θ|A) logp(D | A,θ) − D KL q φ (θ | A)∥p(θ | A) + logp(B | γ)− logp(B | γ) + logp(r)− logq ψ (r) i − E q ψ (r) D KL q φ (γ | r)∥p(γ) .(A.5) The model and guide share the same hard Bernoulli conditionalp(B | γ), so its two log-density terms cancel exactly. Separating the remaining terms for r gives ELBO(φ,ψ) = E q ψ (r)q φ (γ|r)p(B|γ) h E q φ (θ|A) logp(D | A,θ) − D KL q φ (θ | A)∥p(θ | A) i − D KL q ψ (r)∥p(r) − E q ψ (r) D KL q φ (γ | r)∥p(γ) .(A.6) The domain informed prior is now incorporated in p(γ) (see Derivation A.4) as p(γ) = Y i̸=j 1 B(α ij ,β ij ) σ(γ ij ) α ij 1− σ(γ ij ) β ij , σ(x) = 1 1 + e −x , Using Eq. A.3, we reparameterize the expectation overq φ (γ | r)byz ∼ p 0 (z)andγ = T φ (z;r). The conditional KL divergence for the edge logits becomes D KL q φ (γ | r)∥p(γ) = E z∼p 0 h logq φ T φ (z;r)| r − logp T φ (z;r) i = E z∼p 0 h logp 0 (z)− log det ∂T φ (z;r) ∂z − logp T φ (z;r) i . Substituting above equation into Eq. A.6 and changing variables in the inner expectations we get the reparameterized bound ELBO(φ,ψ) = E r∼q ψ (r) E z∼p 0 E B∼p(·|T φ (z;r)) h E q φ (θ|A) logp(D | A,θ) |z expected log likelihood −D KL q φ (θ | A)∥p(θ | A) |z KL divergence for θ + log det ∂T φ (z;r) ∂z − logp 0 (z) + logp T φ (z;r) i − D KL q ψ (r)∥p(r) | z KL divergence for potentials . (A.7) 16 This is our revised ELBO. We defineT φ as either a single flow transformation or a composition T (K) φ K ◦·◦ T (1) φ 1 . In compositional case, the log determinant term in above equation decomposes as log det ∂T φ (z;r) ∂z = K X k=1 log det ∂T φ k (z (k−1) ;r) ∂z (k−1) , z (k) = T φ k (z (k−1) ;r), z (0) = z. A.4 Logistic–Beta Prior on Edge Logits induced by the Domain Informed Prior For each ordered pair(i,j)withi̸= j, let the edge probabilityπ ij ∈ (0, 1)carry a domain informed Beta prior π ij ∼ Beta(α ij ,β ij ), p(π ij ) = 1 B(α ij ,β ij ) π α ij −1 ij (1− π ij ) β ij −1 , with hyperparameters constructed from our beliefs p ij ∈ (0, 1) and concentration ν ij > 0 via α ij = ν ij p ij + 1, β ij = ν ij (1− p ij ) + 1. In particularα ij > 0andβ ij > 0(indeedα ij ,β ij ≥ 1), so the Beta density is proper. Let’s define the edge logit γ ij = logit(π ij ) = log π ij 1− π ij ⇐⇒ π ij = σ(γ ij ) = 1 1 + e −γ ij . The mapγ ij 7→ σ(γ ij )is aC 1 diffeomorphism fromRonto(0, 1), with strictly positive derivative, so the change of variables formula for densities applies. The Jacobian factor is dπ ij dγ ij = σ(γ ij ) 1− σ(γ ij ) ,(A.8) which is positive, so the absolute value may be dropped. By the change of variables formula, the prior density of edge logit γ ij induced by the Beta prior on π ij is p(γ ij ) = p(π ij ) π ij =σ(γ ij ) dπ ij dγ ij = 1 B(α ij ,β ij ) σ(γ ij ) α ij −1 |z from p(π ij ) 1− σ(γ ij ) β ij −1 |z from p(π ij ) σ(γ ij ) 1− σ(γ ij ) |z Jacobian in equation A.8 = 1 B(α ij ,β ij ) σ(γ ij ) α ij 1− σ(γ ij ) β ij . This is the Logistic-Beta density on edge logits. It is a proper density sinceγ ij 7→ σ(γ ij )is aC 1 diffeomorphism of R onto (0, 1) with Jacobian equation A.8, the substitution π ij = σ(γ ij ) gives Z R p(γ ij )dγ ij = Z 1 0 p(π ij )dπ ij = 1. Assuming independence across ordered pairs in the prior, p(γ) = Y i̸=j p(γ ij ) = Y i̸=j 1 B(α ij ,β ij ) σ(γ ij ) α ij 1− σ(γ ij ) β ij . A.5 Closed Form KL for Gaussian p(θ | A) Let d θ be the dimension of θ. Suppose q φ (θ | A) = N μ q (A), Σ q (A) , p(θ | A) = N μ p (A), Σ p (A) , whereΣ q (A) = diag σ 2 φ (A) andΣ p (A) = diag σ 2 p (A) are diagonal with strictly positive entries σ φ,k (A) > 0andσ p,k (A) > 0for allk(so thatΣ q , Σ p are positive definite and invertible), and 17 μ q (A) = μ φ (A). These conditions hold in our model, sinceσ φ = softplus(ρ φ ) > 0andσ p is a fixed positive prior scale. We derive a closed form for D KL q φ (θ | A)∥p(θ | A) = E q φ (θ|A) logq φ (θ | A)− logp(θ | A) . For any mean μ and positive definite covariance Σ, logN (θ;μ, Σ) =− 1 2 d θ log(2π) + log|Σ| + (θ− μ) ⊤ Σ −1 (θ− μ) . Therefore, logq φ (θ | A) =− 1 2 d θ log(2π) + log|Σ q | + (θ− μ q ) ⊤ Σ −1 q (θ− μ q ) , logp(θ | A) =− 1 2 d θ log(2π) + log|Σ p | + (θ− μ p ) ⊤ Σ −1 p (θ− μ p ) , where we suppress the explicit dependence onAinsideμ q ,μ p , Σ q , Σ p for readability. Now taking the expectation under q φ (θ | A), the d θ log(2π) terms cancel and we have D KL (q∥p) = E q [logq− logp] = 1 2 log|Σ p |− log|Σ q | + 1 2 E q (θ− μ p ) ⊤ Σ −1 p (θ− μ p ) − 1 2 E q (θ− μ q ) ⊤ Σ −1 q (θ− μ q ) . (A.9) Let δ = μ q − μ p . Since (θ− μ p ) = (θ− μ q ) + δ and Σ −1 p is symmetric, we have (θ− μ p ) ⊤ Σ −1 p (θ− μ p ) = (θ− μ q ) ⊤ Σ −1 p (θ− μ q ) + 2δ ⊤ Σ −1 p (θ− μ q ) + δ ⊤ Σ −1 p δ. TakingE q and using the Gaussian momentsE q [θ− μ q ] = 0andE q (θ− μ q )(θ− μ q ) ⊤ = Σ q (so that E q [(θ− μ q ) ⊤ M (θ− μ q )] = tr(M Σ q ) for any fixed M ) gives the standard identity: E q (θ− μ p ) ⊤ Σ −1 p (θ− μ p ) = tr Σ −1 p Σ q + δ ⊤ Σ −1 p δ.(A.10) Similarly, E q (θ− μ q ) ⊤ Σ −1 q (θ− μ q ) = tr Σ −1 q Σ q = d θ .(A.11) Now plugging Eq. A.10 and Eq. A.11 into Eq. A.9 gives D KL (q∥p) = 1 2 log|Σ p |− log|Σ q |− d θ + tr Σ −1 p Σ q + (μ q − μ p ) ⊤ Σ −1 p (μ q − μ p ) . (A.12) SinceΣ q = diag(σ 2 φ )andΣ p = diag(σ 2 p )are diagonal with strictly positive entries, we have Σ −1 p = diag(1/σ 2 p ) and log|Σ p |− log|Σ q | = d θ X k=1 logσ 2 p,k − d θ X k=1 logσ 2 φ,k = 2 d θ X k=1 log σ p,k σ φ,k , where the last equality uses logσ 2 = 2 logσ (valid since σ p,k ,σ φ,k > 0), together with tr Σ −1 p Σ q = d θ X k=1 σ 2 φ,k σ 2 p,k ,(μ q − μ p ) ⊤ Σ −1 p (μ q − μ p ) = d θ X k=1 (μ φ,k − μ p,k ) 2 σ 2 p,k . Substituting the diagonal expressions into Eq. A.12 and restoring the explicitAdependence we get, D KL q φ (θ | A)∥p(θ | A) = 1 2 d θ X k=1 " σ 2 φ,k (A) σ 2 p,k (A) + μ φ,k (A)− μ p,k (A) 2 σ 2 p,k (A) − 1 + 2 log σ p,k (A) σ φ,k (A) # . This is the desired closed form. 18 A.6 Pseudo Code and Computational Complexity Algorithm 1 SVI-DAG Algorithm Inputs:DatasetD; priors(α ij ,β ij ); flowT φ ; hypernetworkH φ ; basep 0 (z); iterationsT itr ; Sinkhorn iterationsK S ; objective drawsS MC ; temperaturesT > 0,τ > 0; particles K part ; step sizes η φ ,η r Outputs: Variational parameters φ and posterior particlesr (k) K part k=1 1: Initialize parameters φ and particlesr (k) K part k=1 ∼ N (0,σ 2 r I m ) 2: for t = 1,...,T itr do 3:Compute particle wise objectives and gradients 4:for k = 1,...,K part do 5:Compute P τ,K S (r (k) )←S K S exp(S 0 (r (k) )/τ ) . 6:Compute M (r (k) )← P (r (k) )LP (r (k) ) T . 7:Compute M τ,K S (r (k) )← P τ,K S (r (k) )LP τ,K S (r (k) ) T . 8:for s = 1,...,S MC do 9:Sample z (k,s) ∼ p 0 and set γ (k,s) ← T φ (z (k,s) ;r (k) ). 10: 11:Sample B (k,s) ← 1γ (k,s) + R (k,s) > 0.▷ Eq. 1 12:Sample ̃ B (k,s) ← σ((γ (k,s) + R (k,s) )/T ).▷ Eq. 1 13:Set A (k,s) ← B (k,s) ⊙ M (r (k) )⊙ (1 m×m − I m ). 14:Set ̃ A (k,s) ← ̃ B (k,s) ⊙ M τ,K S (r (k) )⊙ (1 m×m − I m ). 15:A (k,s) ST ← ̃ A (k,s) + sg(A (k,s) − ̃ A (k,s) ).▷ Eq. 2 16:(μ (k,s) φ ,ρ (k,s) φ )← H φ (A (k,s) ST ) and set σ (k,s) φ ← softplus(ρ (k,s) φ ). 17:Sample ε (k,s) ∼N (0,I d θ ) and set θ (k,s) ← μ (k,s) φ + σ (k,s) φ ⊙ ε (k,s) . 18:Draw a minibatchD (k,s) b ⊂1,...,N and compute b L (k,s) ELL ← N |D (k,s) b | m X i=1 X n∈D (k,s) b logp i x (n) i | g i (x (n) ;θ (k,s) i ) . 19:ComputeL (k,s) KL (φ;A (k,s) ST )▷ Eq. 4 20:L (k,s) flow ← log det ∂γ (k,s) ∂z (k,s) − logp 0 (z (k,s) ) + logp(γ (k,s) ). 21: b L (k,s) ← b L (k,s) ELL −L (k,s) KL +L (k,s) flow . 22: b L (k) ← S −1 MC P S MC s=1 b L (k,s) .▷ Eq. 5 23:Target score: g (k) ←∇ ST r (k) b L (k) + logp(r (k) ) . 24:ψ (k) ← vec(M (k) τ,K S ). 25:Update particles via SVGD in relaxed acylicity space 26:for k = 1,...,K part do 27:Φ (k) ← K −1 part P K part j=1 k M (ψ (j) ,ψ (k) )g (j) +∇ r (j) k M (ψ (j) ,ψ (k) ) . 28:r (k) ← r (k) + η r Φ (k) .▷ Eq. 6 29:Update φ with particles held fixed 30: b J block ← K −1 part P k b L (k) . 31:φ← φ + η φ ∇ ST φ b J block . 32: return φ,r (k) K part k=1 . We include the pseudo code forSVI-DAGin Algorithm 1. We choose neural spline flows (Durkan et al., 2019) forT φ as they provide expressive transformations while allowing tractable forward evaluation and log Jacobian computation. The most computationally intensive components are the formation of the relaxed acyclicity mask and particle interactions in the SVGD style update. For each particle, applyingK S Sinkhorn normalizations to anm × mmatrix costsO(K S m 2 ), givingO(K part K S m 2 )work per training iteration. This is followed by the dense multiplication 19 P τ,K S LP T τ,K S that forms the generally cyclic relaxed acyclicity mask, imposing a worst case cost ofO(K part m 3 ). Forming the hard permutation mask can instead exploit permutation structure and does not create this cubic term. The SVGD style update requires pairwise kernel evaluations and gradients across theK part particles. Since the kernel uses vectorized relaxed masksuof lengthm 2 , this step addsO(K 2 part m 2 ). The particle computations before the interaction step can be batched on a GPU. If we are to include theS MC hard/soft edge and adjacency constructions, the overall algorithmic complexity isO(K part m 3 + K part K S m 2 + K part S MC m 2 + K 2 part m 2 ), which aligns with standard cubic complexity of modern approaches to causal discovery (Thompson et al., 2024; Lorch et al., 2021; Bello et al., 2022; Annadani et al., 2023). This excludes the architecture and minibatch dependent costs of the flow, hypernetwork, likelihood, and their reverse mode derivatives, so it is not a bound on total runtime. B Supporting Statements and Proofs Proposition B.1. Let’s fixp ∈ (0, 1)and setα ν = νp + 1andβ ν = ν(1− p) + 1forν > 0. Letπ ν ∼ Beta(α ν ,β ν )and conditionally onπ ν we drawB ν ∼ Bernoulli(π ν ). Thenπ ν → pin probability asν →∞andP(B ν = 1) = E[π]→ p. Hence the Beta Bernoulli distribution converges to Bernoulli(p) in the large confidence limit. Proof. We have E[π ν ] = α ν α ν + β ν = νp + 1 ν + 2 → p, as ν →∞. Also Var(π ν ) = α ν β ν (α ν + β ν ) 2 (α ν + β ν + 1) = ν 2 p(1− p) + ν + 1 (ν + 2) 2 (ν + 3) ∼ p(1− p) ν → 0. The bias-variance identity gives E (π ν − p) 2 = Var(π ν ) + E[π ν ]− p 2 = Var(π ν ) + 1− 2p ν + 2 2 → 0. By Markov’s inequality, for every ε > 0, P |π ν − p| > ε ≤ E[(π ν − p) 2 ] ε 2 → 0. Hence π ν → p in probability. Finally, by the law of iterated expectations, P(B ν = 1) = E[P(B ν = 1| π ν )] = E[π ν ] = νp + 1 ν + 2 → p. Thus the marginal law of B ν is Bernoulli(E[π ν ]) and converges to Bernoulli(p). Proposition B.2. LetXbe a measurable space and letT :X →Ybe a measurable map. Letqand p be probability measures onX with pushforwards q T and p T onY . Then D KL (q∥p) ≥ D KL (q T ∥p T ),(B.1) with equality if and only if a version of the conditional distributions satisfiesq(·| y) = p(·| y)forq T almost everyy. Specifically, for the deterministic mapping(B,r)7→ A = B⊙ M (r), the divergence in the (B,r) space provides an upper bound for the divergence in the adjacency matrix space A. Proof. LetX ∼ qandY = T (X). By disintegration, there exist conditional probability measures q(·| y) and p(·| y) such that q(dx) = q(dx| y)q T (dy), p(dx) = p(dx| y)p T (dy). We begin by expressing the total KL divergence on the spaceX as an integral: D KL (q∥p) = Z X log q(dx) p(dx) q(dx). 20 By substituting the disintegrated forms into the integrand, we have: D KL (q∥p) = Z Y Z X|y log q T (dy)q(dx| y) p T (dy)p(dx| y) q(dx| y)q T (dy). We then use the properties of the logarithm to separate the marginal and conditional components: D KL (q∥p) = Z Y log q T (dy) p T (dy) q T (dy) + Z Y " Z X|y log q(dx| y) p(dx| y) q(dx| y) # q T (dy). The first term on the right hand side is the KL divergence between the pushforward measures, D KL (q T ∥p T ). The second term represents the expected conditional KL divergence, which we denote asE q T [D KL (q(·| y)∥p(·| y))]. The Gibbs inequality states that the KL divergence between any two probability measures is non negative. Therefore, if we assume the measures are well defined, then the conditional term satisfies: Z Y D KL (q(·| y)∥p(·| y))q T (dy)≥ 0. This non negativity implies thatD KL (q∥p) ≥ D KL (q T ∥p T )establishing the inequality. We can observe that for the inequality to become an equality, the conditional term must vanish. This occurs if and only ifD KL (q(· | y)∥p(· | y)) = 0forq T almost everyy, which is equivalent to the condition thatq(·| y) = p(·| y). In the context of our algorithm, this demonstrates that minimizing the divergence in the latent parameter space(B,r)minimizes a conservative upper bound on the divergence between the actual graph distributions. Theorem B.3. Letr ∈ R m have pairwise distinct entries and letB ∈0, 1 m×m have zero diagonal. Let’s define M (r) = P (r)LP (r) T and A = B⊙ M (r). Then the directed graph with adjacency A is acyclic. Moreover for every binary DAG there exist potentialsrand a free edge matrixBsuch that A = B⊙ M (r) (Annadani et al., 2023). Proof. Let’s define the permuted adjacency matrix sorted by the potentials r as A sorted = P (r) T AP (r). SubstitutingA = B⊙M (r) = B⊙ P (r)LP (r) T and usingP T (C⊙D)P = (P T CP )⊙(P T DP ) for permutation matrices, together with P (r) T P (r) = I m , we have A sorted = P (r) T B P (r) ⊙ L. SinceLis strictly lower triangular,A sorted is also strictly lower triangular (i.e.,(A sorted ) ij = 0for alli≤ j). A directed graph whose adjacency matrix can be permuted to be strictly lower triangular contains no cycles, because its edges respect the strict total order on vertices induced by the potentials r. Thus, the graph defined byAis a DAG. Conversely, letGbe an arbitrary DAG with binary adjacency matrixA true . SinceGis acyclic, it admits a topological ordering. Letσbe a permutation of1,...,mconsistent with it, so that every edge ofGis mapped, underσ, into the strictly lower triangular support ofL. We choose the potential vectorr(with pairwise distinct entries, so thatP (r) is well defined) such that its descending sort induces this permutationσ. Under this choice the mask M (r)contains ones at all positions compatible with the topological order, hence the support ofA true is a subset of the support of M (r). Therefore, choosing B = A true , B⊙ M (r) = A true ⊙ M (r) = A true . Thus, every binary DAG can be represented by this construction. Lemma B.4. Letm≥ 2and letr = (r 1 ,...,r m ) T ∈ R m have pairwise distinct entries. We define the score matrixS 0 (r) = ro T , whereo = (m,..., 1) T . Forτ > 0andK S ≥ 1, letP τ,K S (r)denote the output obtained afterK S alternating row and column normalizations applied toexp S 0 (r)/τ , and let P ⋆ τ (r) = lim K S →∞ P τ,K S (r) be its Sinkhorn limit, which exists and is the unique doubly stochastic solution of the entropic assignment problem becauseexp(S 0 (r)/τ )has strictly positive entries (Sinkhorn–Knopp). Let 21 P (r)be the permutation matrix that sorts indices by descending potential, and setM ⋆ τ (r) = P ⋆ τ (r)LP ⋆ τ (r) T andM (r) = P (r)LP (r) T , whereLis the strictly lower triangular matrix of ones. Then P ⋆ τ (r)→ P (r) and M ⋆ τ (r)→ M (r) entrywise as τ → 0. Equivalently, for any iteration scheduleK S (τ )→∞with P τ,K S (τ ) (r)−P ⋆ τ (r) max → 0 asτ → 0, one hasP τ,K S (τ ) (r)→ P (r)andP τ,K S (τ ) (r)LP τ,K S (τ ) (r) T → M (r)entrywise (Annadani et al., 2023). Proof. We denote byΠ m the set ofm× mpermutation matrices. For anyP ∈ Π m associated with a permutation σ of1,...,m, we have the Frobenius inner product ⟨P,S 0 (r)⟩ = m X i=1 m X j=1 P ij r i o j = m X i=1 r i o σ(i) . Sinceois strictly decreasing and the entries ofrare pairwise distinct, the rearrangement inequality implies that this sum is uniquely maximized overΠ m by pairing the largest component ofrwith the largest component ofo, the second largest with the second largest, and so on. The unique maximizer isP (r). First letK S → ∞at fixedτ > 0. By the Sinkhorn-Knopp theorem,P τ,K S (r) → P ⋆ τ (r), the unique doubly stochastic matrix solving the entropic assignment problem P ⋆ τ (r) = arg max P∈P m n ⟨P,S 0 (r)⟩ + τH(P ) o ,(B.2) whereP m is the Birkhoff polytope andH(P ) = − P i,j P ij logP ij is the Shannon entropy. Now letτ → 0. SinceHis bounded on the compact polytopeP m , any sequence of maximizers of equation B.2 converges to a maximizer of the linear objectiveP 7→ ⟨P,S 0 (r)⟩overP m . By the Birkhoff-von Neumann theorem and the previous paragraph this maximizer is the unique vertexP (r). Hence P ⋆ τ (r)→ P (r) entrywise as τ → 0. Finally, for any indices a,b∈1,...,m, M ⋆ τ (r) ab = m X i=1 m X j=1 P ⋆ τ (r) ai L ij P ⋆ τ (r) bj , M (r) ab = m X i=1 m X j=1 P (r) ai L ij P (r) bj . SinceLis fixed and the sums are finite, entrywise convergenceP ⋆ τ (r) → P (r)givesM ⋆ τ (r) → M (r)entrywise. The schedule statement then follows from the triangle inequality together with P τ,K S (τ ) (r)− P ⋆ τ (r) max → 0. C Evaluation Metrics We assess the quality of the learned posterior distribution using four metrics, each capturing a distinct aspects of structured recovery. Letq(G)represent the approximate posterior over graphs and letG GT represent the ground truth causal graph. For each metric, we drawN e posterior samplesG (i) N e i=1 withG (i) ∼ q(G) and report Monte Carlo estimates. Expected SHD (E-SHD). TheSHD(G,G GT )counts the minimum number of edge additions, deletions, and reversals required to transform the estimated graphGintoG GT . Because our method returns a distribution rather than a point estimate, we report the expected SHD under the approximate posterior: E-SHD = E G∼q(G) SHD(G, G GT ) ≈ 1 N e N e X i=1 SHD G (i) , G GT . Lower values indicate that posterior samples are, on average, structurally closer to the ground truth. This metric penalizes both missing and spurious edges as well as orientation errors, providing a comprehensive measure of structural accuracy. 22 Expected F1 score. For each posterior sampleG (i) , we treat every ordered pair(i,j)as a binary classification problem: an edge is either present or absent. Precision is the fraction of predicted edges that appear inG GT , recall is the fraction of true edges recovered byG (i) , and the F1 score is their harmonic mean. The expected F1 score averages this quantity over the approximate posterior: E-F1 = E G∼q(G) F1(G, G GT ) ≈ 1 N e N e X i=1 F1 G (i) , G GT . Higher values reflect better balance between precision and recall across posterior samples. Brier score.The Brier score evaluates the calibration of the posterior edge probabilities by measur- ing the mean squared error between the marginal posterior edge probabilities and the ground truth binary edge indicators. Letˆp ij = E G∼q(G) [A ij ]denote the marginal posterior probability of edge (i,j)and letA GT ij ∈0, 1 denote the corresponding ground truth edge indicator. The Brier score is defined as Brier = 1 m(m− 1) X i̸=j ˆp ij − A GT ij 2 . A lower Brier score indicates superior uncertainty quantification. A perfectly calibrated posterior that assigns probability one to true edges and zero to absent edges achieves a Brier score of zero. AUROC. The AUROC measures the discriminative ability of the posterior edge probabilities ˆp ij i̸=j in distinguishing true edges from absent ones, treating the marginal probabilities as scores for a binary classifier. It equals the probability that a randomly chosen true edge receives a higher posterior probability than a randomly chosen absent edge. An AUROC of100%indicates perfect separation, while50%corresponds to chance level discrimination. Unlike the Brier score, the AUROC is threshold free and invariant to monotone transformations of the scores, making it a complementary measure of ranking quality. CPDAG evaluation. In settings like linear synthetic data, where the DAG is identifiable only up to its MEC, point wise comparison against a single ground truth DAG would unfairly penalize methods that correctly recover the equivalence class but select a different representative. We therefore follow standard practice and convert both the ground truth DAG and each posterior sample to their corresponding CPDAGs (Peters et al., 2017) before computing the metrics above. A CPDAG represents a MEC by directing only those edges whose orientation is shared by every member of the class and leaving the remaining edges undirected. All four metrics are then evaluated at the CPDAG level, ensuring that methods are not penalized for orientation ambiguities that are inherently unresolvable from observational data alone. D Code and license We use the following open source repositaries for comparison with the baselines: • ProDAG: https://github.com/ryan-thompson/ProDAG.jl (MIT license). • BayesDAG: https://github.com/microsoft/Project-BayesDAG (MIT license). • BCD: https://github.com/ermongroup/BCD-Nets (No license included). • DIBS: https://github.com/larslorch/dibs (MIT license). • DDS/VI-DP-DAG https://github.com/sharpenb/Differentiable-DAG-Sampling (No license included). 23 E Detailed Model Specifications and Hyperparameters. Table E.1: Hyperparameters held fixed across the benchmark experiments. Experiment dependent settings are reported below the table and are intentionally excluded from the common entries. HyperparameterValue SVI-DAG Mini-batch size64 Parameter optimizerAdam (Kingma and Ba, 2014) Global gradient-norm bound1.0 Number of ELBO MC samples1 Straight-through hard-adjacency estimatorEnabled Concrete edge temperature T0.3 Final Sinkhorn temperature τ0.1 Number of Sinkhorn iterations K S 100 Sinkhorn implementationLog-space normalization Bayesian node modelTwo hidden layers of width 32 Node-model activationtanh Node-model weight priorN(0, 1) Normalizing-flow familyRational-quadratic neural-spline coupling flow (Durkan et al., 2019) Number of flow blocks5 Number of spline bins8 Spline domain[−5, 5] Number of SVGD particles K part 20 SVGD kernelRBF kernel on relaxed DAG-mask features Order-potential priorr ∼N(0,I), i.e., σ r = 1 Particle-gradient clippingPer-particle ℓ 2 -norm bound of 10 Base SVGD repulsion weight1.0 BayesDAG Maximum epochs for the selected final model150 Maximum epochs for each sparsity-grid proxy fit25 Training batch size128 Number of SG-MCMC chains10 Permutation and SEM SG-MCMC step size3× 10 −4 Variational-helper optimizerAdam, learning rate 5× 10 −3 Permutation SG-MCMC noise scale0.1 SEM-weight SG-MCMC noise scale0.01 Model familyNonlinear ICGNN Number of hidden layers2 Hidden width128 ActivationReLU Layer normalizationEnabled Residual connectionsEnabled Input-permutation conditioningDisabled Variational normalizationEnabled Number of Sinkhorn iterations100 Permutation-relaxation temperature0.2 Sparsity proxy grid4 geometrically spaced values Sparsity-grid range[10, 10 3 ] Posterior samples per proxy fit64 Validation holdout⌊0.1n⌋ observations Selection of λSmallest discrepancy between posterior mean edge count and the true sparsity level Continued on the next page 24 Table E.1 – continued from the previous page HyperparameterValue BCD Nets Number of optimization steps15,000 ELBO MC batch size256 Optimizer for P and LAdaBelief-style Optaxscalebybelief Learning rate10 −3 Number of outer ELBO samples1 Permutation-network architectureTwo-hidden-layer MLP Permutation-network hidden width128 Permutation-network activationGeLU Permutation-logit soft bound10 Gumbel–Sinkhorn temperature schedule30→ 10→ 1, with changes at 5,000 and 10,000 steps Sinkhorn doubly-stochastic tolerance0.01 Posterior permutation samplingHard, with straight-through gradients Posterior over LDiagonal Gaussian over lower-triangular weights and log-noise scales Log-standard-deviation soft bound10 Likelihood-noise parameterizationOne noise scale per node Lower-triangular weight priorHorseshoe Expected in-degree1 Horseshoe global scaleτ H = q d p log(1/q d ), where q d = min2/(d− 1), 1 Log-noise prior standard deviation3.0 Permutation-network ℓ 2 penalty0 DiBS Number of SVGD steps2,000 Number of particles K part 20 OptimizerRMSProp SVGD step size0.005 Likelihood-gradient MC samples128 Acyclicity-gradient MC samples32 Acyclicity inverse-temperature slopeβ(t) = t Gumbel–softmax temperature1.0 Latent dimensionk = d Latent prior standard deviation1/ √ k Latent-graph kernel componentFrobenius squared-exponential kernel Latent-kernel bandwidth5.0 Graph priorErd ̋ os–R ́ enyi, with expected degree parameter 2 ProDAG Maximum number of training epochs1,000 Early-stopping patience5 epochs OptimizerAdam (Kingma and Ba, 2014) Learning rate0.1 Number of ELBO MC samples100 Variational weight priorN(0, 1) Likelihood noise variance1 Initial projected adjacency W (0) 0 Initial acyclicity penalty μ (1) 1 Penalty updateμ (t+1) = μ (t) /2 Number of penalty stages10 Maximum inner projection iterations10,000 Training projection tolerance10 −2 Continued on the next page 25 Table E.1 – continued from the previous page HyperparameterValue Posterior projection tolerance10 −4 Acyclicity threshold0.1 Validation holdout⌊0.1n⌋ observations Number of sparsity-grid values10 Minimum sparsity radius λ min 0 Maximum sparsity radius λ max Mean posterior ℓ 1 -norm obtained with λ =∞ Selection of λMinimum validation reconstruction error DDS/VI-DP-DAG Maximum number of epochs100 Early-stopping patience20 validation evaluations Validation frequencyEvery 2 epochs Training batch size64 OptimizerAdam (Kingma and Ba, 2014) Masked-autoencoder learning rate10 −3 DAG-parameter learning rate10 −2 Training/validation split80%/20% Order relaxationTop-k SoftSort Sampling modeHard, with straight-through gradients Order-relaxation temperature1.0 Noise factor1.0 Masked-autoencoder architectureLinear Hidden-layer widths(16, 16, 16) Training lossELBO Edge regularizerBernoulli KL divergence Edge-regularizer weight0.1 Prior edge probability0.01 Experiment-dependent settings. The common table excludes parameters whose values changed between experiments. For the synthetic Erd ̋ os–R ́ enyi (ER) benchmarks,(p,s)denotes the number of nodes and directed edges, respectively. For reproducibility, the experiment-dependent settings were as follows. SVI-DAG.The iteration budget was6,000for the two-node prior-sensitivity experiments,2,500for the Sachs protein-signaling experiment, and1,500for the linear and nonlinear ER benchmarks with(p,s) ∈ (25, 40), (50, 80). The parameter learning rate and SVGD step size were (10 −3 , 10 −3 )for the two-node experiments and(3× 10 −3 , 10 −1 )for the ER and Sachs experiments. The flow-conditioner widths were(5, 5)for the two-node experiments and (64, 64) for the ER and Sachs experiments. The two-node prior-sensitivity experiments used no straight-through warm-up, full-adjacency conditioning of the node hypernetworks, and no scale normalization of the order potentials. The ER and Sachs experiments used a full-run straight-through warm-up, row-wise adjacency conditioning, and scale-normalized order potentials. The two-node experiments additionally capped the SVGD repulsion norm at the attraction norm; this cap was disabled for the ER and Sachs experiments. The initial Sinkhorn temperature was20for the Sachs experiment and0.1for the two-node and synthetic ER experiments; the final Sinkhorn temperature was0.1in every experiment. The two-node and Sachs experiments learned graph-conditioned likelihood-noise scales using a log-noise prior centered atσ = 0.1. The linear and nonlinear ER benchmarks with (p,s)∈(25, 40), (50, 80)instead used fixed likelihood noiseσ = 0.5. The weight-space KL coefficient was1for the two-node and Sachs experiments and0.01for all synthetic ER experiments. BayesDAG.The same nonlinear ICGNN architecture was used for the linear ER, nonlinear ER, and Sachs experiments. Thus, BayesDAG did not select a separate architecture according to whether the data-generating structural equations were linear or nonlinear. 26 DiBS.MarginalDiBS with the BGe score, inverse-temperature slopeα(t) = t, and the score- function gradient estimator was used for the linear ER benchmarks with(p,s) ∈ (25, 40), (50, 80)and for the Sachs experiment. JointDiBS with a nonlinear Gaussian likelihood was used for the nonlinear ER benchmarks with(p,s) ∈ (25, 40), (50, 80). The nonlinear model used a single hidden layer of10ReLU units, inverse-temperature slopeα(t) = 0.05t, and the reparameterization gradient estimator. It additionally used parameter-kernel bandwidth h θ = 500 and likelihood-noise standard deviation 0.1. ProDAG.The linear ProDAG model with projection step size1/dwas used for the linear ER benchmarks with(p,s) ∈ (25, 40), (50, 80)and for the Sachs experiment.The nonlinear ProDAG model was used for the nonlinear ER benchmarks with(p,s) ∈ (25, 40), (50, 80). It used one hidden layer of10ReLU units and projection step size 0.25/d, where d denotes the number of variables. BCD Nets and DDS.These methods used the same linear structural model for the linear ER, nonlinear ER, and Sachs experiments. In particular, neither wrapper selected a separate nonlinear architecture for the nonlinear ER benchmarks. F Broader impact & limitations F.0.1 Broader impact statement This work is concerned with understanding cause and effect relationships from data, with potential applications across empirical sciences, economics, epidemiology, and climate science. By explicitly quantifying epistemic uncertainty over causal structures and providing a principled mechanism to incorporate domain knowledge, our approach can support more calibrated decision making, expose biases present in the data, and allow reliable answers to causal queries in settings where data are scarce or the underlying DAG is not identifiable. As such, we envision this line of work to not have any significant negative impact. F.0.2 Limitations Though our approach has several strengths, it naturally has limitations. As with most Bayesian approaches that rely on variational inference, our framework lacks theoretical guarantees for exact posterior approximation. In particular, while the combination of normalizing flows and SVGD in the relaxed acyclicity space is designed to promote mode coverage, our method does not guarantee recovery of all DAGs within MEC. It can only aim to capture as many of them as possible, and some equivalent structures may remain underrepresented in the learned posterior. 27