BCD Nets: Scalable Variational Approaches for Bayesian Causal Discovery
Authors: Chris Cundy
Aditya Grover
Stefano Ermon
Venue: NeurIPS 2021 Poster (2021)
Keywords: Bayesian
DAG
Variational Inference
Causal Discovery

Abstract: A structural equation model (SEM) is an effective framework to reason over causal relationships represented via a directed acyclic graph (DAG).
Recent advances have enabled effective maximum-likelihood point estimation of DAGs from observational data. 
However, a point estimate may not accurately capture the uncertainty in inferring the underlying graph in practical scenarios, wherein the true DAG is non-identifiable and/or the observed dataset is limited.
We propose Bayesian Causal Discovery Nets (BCD Nets), a variational inference framework for estimating a distribution over DAGs characterizing a linear-Gaussian SEM.
Developing a full Bayesian posterior over DAGs is challenging due to the the discrete and combinatorial nature of graphs.
We analyse key design choices for scalable VI over DAGs, such as 1) the parametrization of DAGs via an expressive variational family, 2) a continuous relaxation that enables low-variance stochastic optimization, and 3) suitable priors over the latent variables.
We provide a series of experiments on real and synthetic data showing that BCD Nets outperform maximum-likelihood methods on standard causal discovery metrics such as structural Hamming distance in low data regimes. 

BCD Nets: Scalable Variational Approaches for Bayesian Causal Discovery
Chris Cundy!
Aditya Grover Stefano Ermon
Department Computer Science, Stanford University Facebook Al Research
University of California, Los Angeles
{ (cundy, ermon @cs. stanford. edu adi tyagees wwtggesssucll..cm ucla.
Abstract
A structural equation model (SEM) is an effective framework reason over causal relationships represented via a directed acyclic graph (DAG). Recent advances have enabled effective maximum-likelihood point estimation of DAGs from ob- servational data. However, a point estimate may not accurately capture the un- certainty in inferring the underlying graph in practical scenarios, wherein the true DAG is non-identifiable and/or the observed dataset is limited. We propose Bayesian Causal Discovery Nets (BCD Nets), variational inference framework for estimating a distribution over DAGs characterizing a Developing a full Bayesian posterior over DAGS challenging due to the the discrete and combinatorial nature of graphs. We analyse key design choices for scalable VI over DAGs, such 1) the parametrization of DAGs via an expressive variational family, 2) a continuous relaxation that enables low-veriance stochas- tic optimization, and 3) suitable priors over the latent variables. We provide a series of experiments on real and synthetic data showing that BCD Nets outper- form maximum-likelinn likelihood methods on standard causal discovery metrics such as structural Hamming distance in low data regimes.
SEM.
1 Introduction
One of the key uses of statistical methods is learning causal relationships from observed data a.k.a. causal discovery [39]. Causal models allow forecast the effects of interventions and counter- factuals several real-world domains, such as economic policy [57] and medicine [47]. Although early approaches to statistical inference emphasised that 'correlation is not causation' [15], it has since been shown that for certain families of data-generating processes, it is indeed possible to infer causal relationships from purely observational data [9].
One such widely studied data-generating process is the linear-Gaussia structural equation model (SEM) [37]. where the causal relationships between the random variables in the model can be rep- resented via a weighted directed acyclic graph (DAG). The value any variable the DAG of Sinear-Gaussin SEM is given by a lincar combination of the values of its parent nodes and ad- ditive noise. For causal discovery, naive Monte-Carlo sampling or enumeration of the possible DAGS quickly becomes intractable, since the number of possible DAGs over model grows super- exponentially with the number of variables [18]. A variety of methods have been developed over the years efficiently sample or optimize over DAGs [11, 53, 51, 45J. For example, a recent line of work scales to high dimensions by maximizing the likelihood of the model (MLE) over set of continuous relaxations of adjacency matrices using gradient-baser methods and specialized DAG regularization terms. [62, 61, 35, 59].
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

D3
31
0
0
0
0
31
22
1.3
0
D22
J33
0.3
2.5
0
22
24
F4
0
5.7
0
X44
Figure I: A Structural Equation Model (SEM). Left: DAG with 4 nodes. Right: Linear-Gaussiar
SEM (with (0,
However, the majority of the aforementioned works for scalable causal discovery in SEMs focus recovering point estimates for the underlying DAG via MLE. In many practical scenarios, however, a point estimate fails to reflect the uncertainty in inferring the underlying DAG. This includes scenarios where the true DAG is non-identifiable given (infinite) observational data as well as practical limitations due to an imperfect optimization algorithm, model mismatch, or simply finite dataset. In any of the above scenarios, it is desirable to obtain explicit posterior distribution over the unobserved DAG instead of single point estimate [19]. Causal inference increasingly being applied in situations with important real-world consequences, where such Bayesian estimation procedure could be useful to sample and reason over alternative generative mechanisms for observed data.
To achieve this goal, we propose Bayesian Causal Discovery Nets (BCD Nets), an algorithmic framework for Bayesian causal discovery in linear-Gaussian SEMs based on modern variational inference [6]. Classical approaches to Bayesian causal discovery struggle in higher dimensions [23, 16], with limited improvements via Monte Carlo approximations [22, 58]. In order to scale our variational method to high dimensions, we address several design challenges. First, we describe an expressive variational family of factorized posterior distributions weights and noise variance) using deep neural networks. The factorization exploits the decompo- sition of DAGs into triangular matrices and permutations for specifying the distribution weights and node orderings respectively. Further, for low-variance stochastic optimization of the variational objective, we exploit recent advances modeling and reparametrizing distributions over permutations via continuous relaxations [31, 28]. Finally, we employ a horseshoe prior [7] on the edge weights, which promotes sparsity.'
the SEM parameters (edge
edge
We evaluate BCD Nets for causal discovery on a range of synthetic and real-world data benchmarks. Our experiments demonstrate that with finite datasets, there is considerable uncertainty in the in- ferred posterior over DAGs. Using BCD Nets, we are able to effectively quantify the uncertainty and significantly outperform competing estimators [35, 58] the standard Structured Hamming Distance metric, especially in low data regimes.
2 Preliminaries
2.1 Linear-Gausiia Structural Equation Models
A structural equation model (SEM) is collection of random variables ..,....22 zd associated with a directed acyclic graph (DAG) G with d nodes [37]. The SEM consists of a series of equations r, = (a()))))))) where Pa(x;) gives the values of the parents of the jth node in G and e; is noise variable. For a linear-Gaussian SEM, the equations fj are linear and the noise ej is additive Gaussian. Considering X [ xa] as a vector IRd and W as the weighted adjacency matrix of G. we can see that X must satisfy X = WT X + e, where ~ NO, Σ) is the additive noise vector and Σ = diag {oî .o.,...,o7 ..) is a diagonal noise covariance matrix. We will denote the setting ad
when all noise variances are equal i.e.,
"equal variance" setting and "non-equal variance" otherwise. Figure/ shows an illustration. Linear SEMs have been used to model microarray data [38] and protein pathways [13], among many other systems.
'Code is available github, com/ wwwongroup/BDD--Meee.

---

2.2 Causal Discovery
Given some dataset X = X2......... drawn i.id. from a linear- Gaussian SEM, we are interested in inferring the adjacency matrix W and the diagonal entries of the covariance matrix Σ.
Maximum Likelihood Estimation (MLE) This approach obtains point estimate for W and E by maximizing the likelihood of the dataset X Using the fact that W is constrained to be a DAG, can be shown that 7- WT is always invertible.?. Hence, we can rearrange to obtain X = (7- - -
obtaining X as the product of a matrix and Gaussian vector. It follows that X itself has a Gaussian distribution, X ~ (0, 9-1), with precision matrix θ = (I - W-E W) T_ This gives a joint log-likelihood over the dataset Xi as
log logp(X?;E, W)
n 2
(logdet O dlog(2T))
1 2
Ox,,
(I)
where X, runs over the points in the dataset and θ is given above. The space of DAG adjacency matrices W is characterised by the acyclicity constraint. When viewed as an subset of the space of all adjacency matrices (Rdxd), this is a piecewise linear manifold with a number of facets that grows super-exponentially in the dimension (roughly d!2d//2), rendering any approach based on enumeration over all DAGs intractable. Many approaches have been proposed to scale causal discovery via MLE to high-dimensiooal data. This includes recent approaches which relax the DAG constraint to larger, continuous set of adjacency matrices, such as Rdddd in [62] or the Birkhoff' polytope of doubly-stochastic matrices in [5], which enables the use of tradient-based optimization.
However, there are a number of challenges with using point estimators for causal discovery. Fun- damentally, the true SEM parameters identifiable from observational data only under specific conditions. In fact, the map between the the data distribution parameters e and the SEM parameters (W, } may not be bijective even in the limit of infinite data and an oracle optimizer. Theoreti- cal results on identifiability in Amear-Gaussian SEMs currently hold only under restricted settings; notably, this only includes the equal variance setting [38]. Further, even if the noise variances equal in the ground-truth SEM, we are limited in practice by the finite size of the dataset and MLE could still converge to an incorrect solution. Finally, such point estimators are typically not robust to potential misspecifications in the likelihood model, e.g., a non-Gaussian noise model, non-linear cause-effect relationships, etc. For the above scenarios, it is useful characterize the uncertainty in estimating the SEM parameters to guide downstream analysis. For example, in medical applications giving a distribution of possible causal pathways, with quantified uncertainty, could be much more useful than a point estimate of the most likely pathway.
Bayesian Estimation. In contrast to point estimators, Bayesian methods explicitly characterize the uncertainty in the estimated parameters [19]. That is, we treat the unknown SEM parameters
{ }} as random variables associated with a prior distribution p(W, E). The likelihood model P(XT|E, W) follows the same expression the RHS in equation (1), Given the prior and the likeli- hood, we obtain a posterior distribution p(W, Σ [ over { [W. Σ} via Bayes rule, which quantifies the uncertainty in estimating the SEM parameters. As long as the likelihood model is well-specified and the prior includes the ground -ruth SEM parameters in its support, as the dataset size increases the posterior will concentrate around set of SEM parameters. In the equal variance case these will be the true SEM parameters. In the non-equal variance case the posterior will concentrate on the set of DAG parameters quasi-equivalent to the ground truth defined in [35]. These the set of DAG parameters which generate data with the same covariance as the ground truth, and are hence indistinguishable based on data alone. We discuss this further in the appendix, Section F.
The key challenge in Bayesian estimation tractable computation of the posterior distribution in high-dimensional spaces. With the exception of specific prior and likelihood families e.g., conjugate distributions, computing the posterior is typically intractable. In the next section. we present a variational framework for scalable approximation of the posterior for Bayesian causal discovery.
Wriitng W PLP from Section 3. I T (1 PL PL = P(I L -_)PP. Now (I L upper- trian-triangula with unit diagonal, det(1 )T))
L = 1. Thus (1 (1-
( W has unit determinant, is invertible.
3

---

Causal Discovery via Bayesian Causal Discovery Nets
As discussed previously, we are interested in learning the posterior distribution p(W, X over the unknown SEM parameters {W, Σ} given an observed dataset X Unlike point estimators, such a posterior distribution will allow us to quantify the uncertainty estimation. Our framework, Bayesian Causal Discovery Nets (BCD Nets) allows us to tractably estimate this posterior.
As first step, our approach involves parametrizing the adjacency matrix W the product of per- mutation matrix P and a strictly lower-triangular matrix L, so that W = -PLPT. In graphical terms, L is weight matrix for canonical DAG with fixed ordering, while pre- and post-multiplication by P and PT modifies the ordering of nodes. L is parameterised by vector of weights / IRd(d-1)/2 and the constraint that P is a permutation ensures that W is the adjacency matrix of a DAG.
Our goal is to obtain the posterior distribution p(P, L, Σ|XT). Due to the intractable partition func- tion, we cannot directly compute the posterior. We turn to variational inference to deliver a tractable approximation to the posterior [24]. The key idea here is to cast inference as optimization prob- lem, wherein we approximate the true posterior with a tractable family of distributions qo (o(P, P, L, Σ) parameterized by ø and optimize these parameters to minimize the KL divergence between the approximate and true posterior distributions:
DRL. (Po L, L, xf))
(o L,,P,L, 2) p L, E)
(.B(LL, (ogp(XT(P,,,E)) L, log
logg(X?)
ELBO(Ø)
Hence, minimizing the KL divergence above corresponds to maximizing the evidence lower bound (ELBO) w.r.t. variational parameters ø. With a sufficiently expressive variational family from which to choose q, maximizing the ELBO recovers the true posterior 9a (P, L, Σ) = p(P, L, E|X}).
In practice, we face important modeling choices which have a substantial impact on the quality of the posterior obtained, as well as the difficulty of optimizing the ELBO. These include the variable ordering to use when factorizing qo (P. I., Σ) using the chain rule, choice of variational family for the individual (conditional) factors, as well the prior distribution p(P. E).. We discuss these algorithmic design choices next.
3.1 Factorization of Approximate Posterior
Approaches to variational inference with multiple sets of latent variables often use mean-field factorization [24], in our case corresponding to go P. L, =
This mean-field approach can often simplify the optimization of the ELBO, but severely limits the expressiveness of the approximate posterior. For example, consider two-dimensional linear-Gaussial SEM with non- equal variances (and therefore has non-uniquely identifiable parameters). Under infinite data, the posterior density concentrates on the two (observationally undistinguishable) maximum-likede likelihood solutions: edge I I2 with some weight ajand edge posterior concentrates to a bimodal distribution, with density around the region (41 P1) and around (12, Pz). A meeanfieel factored posterior qs (P) cannot represent this correlated density. Em- pirically we observe that such a factored posterior leads to a worse ELBO, illustrated ablation experiments in Section 5.6.
I2 with another weight l2. The
For BCD Nets, we use factorization qs (P, E. = qo qo E) [,, E), sampling I. and Σ jointly first, then conditionally sampling P based on these values, using a neural network to learn the parameters of the conditional distribution qo (PI,, E). This leads to an ELBO
(, log (P(LL, Σ)
log (L, p(L, 2)
Ep Ep-gs(-1L,S))
(3)
D(PL,
3.2 Variational Families
An important design choice in variational methods is the variational family used. A distribution 4o over latents z used in ELBO optimization must support two operations: drawing sample = ~ fa- and computing log (2) 1(z). Depending on the prior distribution it may be additionally

---

possible to compute the term Ex [log (o(2) log = DXL 19s. p] in closed form, possibly reducing variance compared Monte Carlo estimates [43]. It also desirable that sampling qo can be written as 9o(7), 2o. i.e. that is obtained by sampling from fixed distribution 4o and transformed through a parameterized, differentiable sampling path Pg. This lets us use pathwise gra- dient estimators [34], which typically have lower variance than the score-function alternatives [60].
3.2.1 Distribution over Weights & Noise Variances
We consider two different variational families for the distribution over weights and noise variances, Ls,L, Σ), depending on the modeling assumptions over the noise. Under the equal variance model- ing assumption, we parameterize our variational family as a (diagonal covariance) normal distribu- tion, with & directly encoding the mean and variance of the ddd 1)/2 random variables for L and the single random variable for We use this simple distribution since expect the posterior over to be relatively unimodal in the identifiable equal variance case.
In the non-equal variance case, and in the experiments using real-world data, we use normaliz- ing flow [41] for (L, Σ). We expect that in this case the distribution over L could be much more complicated, and so more expressive density model is desirable. We use continuously indexed nor- malizing flows [10], a recently-developer family of flows offering good performance on multimodal densities. As desired, both the normal distribution and flows have pathwise gradient estimators.
3.2.2 Distribution over permutations
Since the set of dddimensional permutation matrices Pa is discrete and its size scales combinatorially with d, it challenging to specify a variational family of distributions over Pa that permits both density estimation and sampling for low-variance stochastic optimization of the ELBO objective in equation (3). Since permutations are discrete, pathwise gradient estimators do not exist. Hence, we consider relaxations to distributions over permutations. Our base distribution is the Boltzmann distribution over Pa. parametrised by T E Rddd with probability (P)(P) x exp(T, , P) for P E Pa-
Density estimation. Computing the partition function for the Boltzmann distribution, Epep, (P) involves an expensive enumeration and is therefore intractable evaluate in high dimensions. In order to approximate the partition function, we follow [28] in noting that the partition function is equal to the matrix permanent perm (exp T). This can in turn be approximated tractably via the Bethe permanent, denoted as perm g (exp T). The Bethe permanent is known to satisfy log perm T d log 2 < log perm n T < log perm T, so that the density will be over-estimated, by more than factor of 1 log 2 [3]. We refer the reader to appendix C in [32] for efficient implementation of the Bethe permanent estimator based on message passing.
Pathwise Gradient Estimation. Exact sampling from the Boltzmann distribution is challenging for similar reasons as exact density estimation. Moreover, even tractable low-rank approximations to the Boltzmann distribution based on Gumbel-Matching distributions [55] are not useful they involve non-differentiable operations and so cannot be used to derive a pathwise gradient estimator. Instead, use relaxation to the Gumbel-Matching distribution, the Gumbel-Sinkhom distribution [31].
To draw a sample from the Gumbel-Sinkhorn distribution with parameters T, we calculate S((T + 7)/r), with S the Sinkhorn operator [50], a matrix of i.idd standard Gumbel noise and 7 a tem- perature hyperparameter. S(T) returns the fixed point obtained from repeated row and column normalization, starting from the elementwise exp of T. In the limit of an infinite number of itera- tions, this returns doubly stochastic matrix. As T approaches zero, the samples approach samples from Pa. with a distribution given by the Gumbel-Matching distribution. A proof of this fact is given in the appendix of [31]. As the Sinkhorn algorithm is a differentiable function of standard Gumbel noise, Himbel-Sinkhorn distribution. Additional implementation details for our Sinkhorn approach are in the appendix. Section B.
use a pathwise gradient estimator of gradients involving samples from the
3.3 Prior Distributions
A key aspect of any probabilistic model is the choice of prior distribution for the unknown param- eters. The prior incorporates domain knowledge into the problem. Moreover, specific choices of prior can be computationally friendly.
5

---

We show in the appendix (Section A) that if we choose the prior over edge weights be an isotropic Gaussian, we can analytically marginalize out the weights. only requiring the distribution over P to characterise the full posterior. Once we have P. it is straightforward to obtain T, since it regression problem which can be solved tractably [54]. Although
Gaussian Prior:
very convenient to avoid modeling a distribution over L, in practice we find that the Gaussian generative assumption on the weights is not particularly useful for datasets which we would like to analyse, for which the underlying DAGs typically sparse. A sample from the distribution over DAGs with Gaussian edge weights will likely have many large- or moderate- edges.
Previous work finding the maximum-likelie likelihood solution of equation (1) has added a term A|| WIl1 penalizing the L1 norm of the adjacency matrix [35, 62], which can be interpreted as imposing an isotropic Laplace prior. The Laplace prior is known to induce sparsity in the posterior [54], but there is generally no way to choose the regularization coefficient À without cross-validation.
Laplace Prior.
Horseshoe Prior. Given the limitations of the above choices of priors, we instead propose to use horseshoe prior on L. The horseshoe prior has a sharp peak at zero and relatively flat tails which tend to induce sparsity in the posterior while not significantly penalizing larger coefficients [7]. Mathematically, variable 3; has horseshoe distribution it the result of first drawing a random variable X C+ (0. 1) from a half-Cauchy distribution, then sampling 3g ~ (o, 1²) The parameter 17 can be adjusted to encode the prior belief on the degree of the DAG generating the data. Based roughly on [40], we suggest a rule of thumb of setting 1) = p/dv/n, with the prior belief of the average degree of the DAG. This results in a sparsity prior that penalises more stringently with more data, similarly to the BIC score [46].
3.4 Overall Approach
Incorporating all the above design choices, we obtain our overall algorithm for BCD Nets. The pseudocode is shown in Algorithm 1. Here, qo [L, Σ) is parameterized as either normal distribution or normalizing flow depending on the modeling assumption of equal non-equal variances. The distribution es (PIL, Σ) is Gunbel-Sinkhorn relaxation of the Gumbel-Matching distribution over permutations, parameterized by function ha conditioned on L. For the function hg. we use a simple two-layer multi-layer perceptron. For stochastic optimization w.r.t. this distribution, we additionally find useful to use the straight-through gradient estimator [4]. This means that on the forward pass of the backpropogation algorithm, we obtain the - O limiting value of the Sinkhorn relaxation using the Hungarian algorithm [27], giving a hard permutation P. On the backward pass the gradients are taken with respect to the finite-r doubly-stochastic matrix P. The prior is given by p(P, L, E)= = ((P)p(L)p((5 where p(L) is horseshoe prior, p(P) is a uniform prior over permutations and p(E) is relatively uninformative Gaussian prior on log
Σ.
Algorithm 1: Bayesian Causal Discovery Nets (BCD Nets)
Input : data XT, Gradient-based optimizer step, temperature hyperparameter r Initialize parameterized distribution qg. neural network h& (L, while not converged do
Draw L, >.
go 2)
Compute logits T = h (L,
Draw 7 E Raddd i.i.d from standard Gumbel
Compute soft P- S((T TT), hard P Hungarian(())
Compute g = V [ELBO()] from equation (3) with sampled P, L, using P in the forward pass and P in the backward pass Update e via step using gradient
end
4 Related Work
Structure Learning for Bayesian Networks: The field of Bayesian structure learning investigates how to infer the structure of Bayesian networks from data. Learning the structure of graphical models from data is known to be NP-hard [8]. Nevertheless, two main families of approaches have

---

been developed to tackle this problem. The first are constraint beatraiint---edddddd approaches [52], aiming to infer conditional independences from the data and learn a structure consistent with these independences.
The second family of approaches are score-based methods [30, 26], which assign a score to a pro- posed DAG based on the goodness-of-fit to the data, and search over the space of DAGS to find a structure which maximizes the score. The combinatorial size and discrete nature of the set of possible DAGs in dimensions present the main difficulties to this approach. Several competing scores have been proposed [46, IJ. which generally have an either implicit or explicit term penaliz- ing complexity to avoid overfitting to the data. Recent work has focussed on developing approximate algorithms that can return high-scoring graph with probabilistic guarantees [36] and on exact meth- ods for graphs under structural assumptions. such as an upper bound on the number of parents of a node [12. 45], similar to the classic Chow-Liu algorithm [9]. Using these techniques it is possible to search over all graphs where nodes have at most two parents for dimension
1000 [48].
Another set of works improve sampling from distributions over DAGS, e.g. approaches sampling over orderings instead of over DAGs [53]. Given a partieular network structure, there are several approaches to learning DAG edge weights, even in the presence of latent variables [33].
Continuous Relaxations for Structure Learning: The problem of learning the structure of prob- abilistic graphical models on undirected graphs can be formulated as a convex optimization prob- lem [17], which allows fast inference even for large graphs. However, the acyclicity constraint for learning directed graphs means continuous optimization cannot be applied directly. In recent years, eaase approaches structure learning have emerged, based the idea of learning the adjacency matrix corresponding to a DAG via optimization over a matrix W E Rdxd, instead of in the set of adjacency matrices corresponding to DAGs. This approach was introduced in Zheng et al. [62], with the least-squares loss || - W 11 used to measure the fit of the data to W. To encourage W approximately form an adjacency matrix, penalty h(W) was added the loss, where hWW) = 0 only if W is a DAG. Originally the penalty h(W) Tr Trewew wew was used, and subsequently different penalties have been proposed [61].
Particularly of interest to us is GOLEM [35], which specifically investigates relaxed optimization approaches for the finear-Gausian SEM that we consider. They point out the similarity of the squared loss to the log-likelihood in equation (1), and show that the additional term log (T W) also serves as a penalty for non-DAG W. This allows optimization over W maximize the likelihood without the augmented Lagrangian optimization procedure used in Zheng et al. [62].
Learned Orderings:
There are strong connections between our SEM setting and the problem of
learning an ordering for autoregressive models. An SEM is similar to an autoregressive model [2] in the sense that the value at the ith index depends on (possibly) all of the values before i. It has been observed that certain orderings for autoregressive models can lead to higher likelihood [20], and subsequent work [25, 28] has developed methods to find the orderings which lead to the best performance from the autoregressive model. The closest to our work is Li et al. [28], which uses latent permutations to learn ordering for an autoregressive model. However, this approach doesn't 't allow reparameterised gradients due to the hard assignment for the autoregressive order, and the REINFORCE [60] gradient estimator is used instead of lower-variance reparameterised gradients.
5 Experiments
In this section study the empirical performance of our method. On synthetic data, we show that our distributional approach outperforms baselines, including the maximum likelihood approach, in the low-data regime in terms of edges identified correctly. On a real-word protein dataset [44], we also identify more edges correctly compared to the maximum-likelibon causal inference problem, we outperform competing methods. Finally, we carry out an ablation with various parts of our algorithm, showing the degradation in performance if we remove key parts. In all cases we use the structural Hamming distance (SHD) [56] to quantify how close certain DAG is to another. The SHD measures how many insertions, deletions or flips are required to turn one graph into another. We also use the CPDAG SHD (SHD-C) [56] which computes SHD up to the equivalence class of completed partially directed acyelic graphs (CPDAG)s [42]. We give additional experimental details in the appendix B. In every case when we give the SHD produced by our model, this represents the SHD marginalized over our posterior, obtained by taking 100 samples of W from
method. On a toy
7

---

ER1, Gaussian Noise
ER2, Gaussian Noise
ERI. Gumbel Noise
ER2, Gumbel Noise
250
250
250
250
Bila Nets joung
GAEGET
Structure Hemminn Distance
200
but Rusckin
200
200
200
150
150
100
100
50-
50
50-
50
0
0-
20
40
60
20 40
60
20
40
60
20
40
60
Figure 2: Distance between true graph and that estimated from BCD Nets, our variational approach, compared to several baselines described in the main text. Lower is better.
our converged posterior, computing the SHD and taking the average of these SHDs. Error bars are given by 5 random seeds for the data generation choice of data subset.
5.1 Synthetic Data
We study our model's performance in the low-data regime, with = 100 data points. We expect that when there little data, the maximum-lii -ikeliiooo DAG may not be particularly representative of the DAGs that could have plausibly made the data, and so a Bayesian approach which can assign posterior density to multiple candidates may correctly identify more edges. In this experiment, we draw Erdös-Renyi graphs [14] with average degree equal to or 2, denoted ER1 and ER2. The weights are generated uniformly in (-2. -0.5) U {0.5. 2} following [35]. The data is generated in the equal-variance case so that the ground-truth is identifiable and the SHD is meaningful, We evaluate with Gaussian noise as well as with Gumbel noise to test how our method behaves under a misspecified likelihood. The distribution over the weights also induces some misspecification, since the weights are not drawn from the modelled prior of a horseshoe.
We compare against GOLEM [35]. which attempts to find the maximum-likedn likelihhodd solution of equation (1) using the continuous relaxation introduced in [62]. We found that the hyperparameters given by the authors did not lead to optimal performance with low data, so we found the best ones via a valss-vvliiation scheme given in appendix B. We also include results from GADGET [58], a recent highly-optimized state-of-the-arr Monte Carlo method for Bayesian causal discovery for DAGs, using 300,000 iterations with 48 parallel chains. To compare against a baseline that is widely used in practice, we also evaluate the DirectLiNGAM method [49], designed for non- Gan-Gansiinn likelihoods. Also included are two methods designed specifically for Gaussian DAGs, the greedy directed acyclic graph search with equal error variance (GDSEEV) method [38]. and the main method ( ("Ghoshal') from [21]. For the latter, we were unable to find an implementation by the authors and used our own implementation. We use the equal- vgrianva formulations of GOLEM and BCD Nets.
Results are shown in figure 2. The null baseline is a method that simply predicts no edges. We observe that we reliably obtain significantly lower SHD than MLE approaches, especially at higher d and degree where we tend to recover a significant fraction of the edges. For example, in the Gaussian noise case with average degree 2. for dimension 64 we get an SHD of 59 + 14 compared to GOLEM's 205 + 30. This corresponds to double the true positive rate while having a false discovery rate ohe-third of GOLEM's. We give full plots of the true positive rate, false positive rate and false discovery rates in the appendix, Section D. The Bayesian GADGET method also obtains better results than MLE methods, with BCD Nets' match to this exact sampling method suggesting that our approach finds a variational posterior close to the real posterior. Compared to GADGET, our method is generally faster (see table 2) and gives an explicit parametric form for the posterior, instead of simply giving samples. Furthermore, the ELBO is a natural metric to evaluate convergence, while for GADGET it is unclear how many iterations are required to mix?. DirectLiNGAM is competitive on the Gumbel likelihood, fitting its linear non- -anussssi framework.
An earlier version of this paper reported worse SHDs for GADGET after using 50,000 steps. Increasing the number of iterations by 6 times resulted in much better mixing
8

---

5.2 Protein Dataset
We also evaluate on a benchmark protein signalling dataset [44]. The d 11-dimensional dataset consists of n = 853 observations, with an expert-provided ground-truth graph. The true structure has 17 edges. We train on random draws of 100 observations. The results are shown in table For GOLEM, both the equal-variance and non-equal variance maximun-likelimee worse than predicting no edges at all. Our equal variance method does not perform much better. Our non-equal variance method performs much better than maximum likelihood, obtaining an SHD of 15 from only 100 samples, close to the SHD of 14 obtained by GOLEM-NV when using the entire 853 data points (reported in Ng et al. [35]). NOTEARS, from [62] performs better than GOLEM. Meanwhile, the score-based GES method also achieves an SHD-C of 14, showing that our method approaches the performance of methods which enumerate all graphs (but scale badly). The non- probabilistic DLiNGAM method achieves similar SHD, although has the lowest SHD-C, indicating gets some edges correct but the direction wrong. Finally, GADGET [58] achieves slightly lower SHD of 13.9, within error of our approach and GES.
methods perform
5.3 Causal Inference
To illustrate how our distributional approach could be used for causal inference, we test the ability of our method to make interventional predictions. In this experiment we generate a synthetic ER graph of degree 1 as in Section 5.1. We then choose a random edge in the graph between nodes i. with I; + Zz.- Using the ground-truth parameters W*, we choose random value a and sample zy do(x; alW*, ,E), by directly intervening in the data- generating process. For estimated W. E we can also sample I ~ do(x; alW, Σ). For our distributional approach, we marginalize over the final posterior distribution of parameters qo (W, Σ), drawing zy from the distribution with probability Ew,z~q, [P(do(z; a W. Σ))] We then compare the sampled empirical distribution of zj to the ground truth interventional distribu- tion. Our marginalization over the posterior allows to get significantly closer the interven- tional distribution, measured by the Wasserstein distance (e.g. 0.25 for the Thear-Gaussian graph at d 64, compared to 2.8 for GOLEM). Full results as function of d are given in appendix E.
Table I: Causal discovery approaches on the pro- tein dataset with reduced data (n. 100)
# Edges
SHID
SHD-C
GOLEM-EV GOLEM-NV NOTEARS
1,5 1.3
18.5 + 1.3
18.5 + 1.3
1.5 + 1.3
18.5 + 1.3
18.5 + 1.3
18.5
16.5 + 0.9
GES
12
14.6 2.0
PC
0.5
1.6
GADGET
0.5
1.29
0.66
DLINGAM
4.6
14.8
0.5
BCD Nets-EV
11.3 1.2
0.3
19.4
BCD Ncts-NV
9.0
14.7
14.0 + 1.0
5.4 Model Running Time
We expect the time taken to train our model will asymptotically scale as (dd), similarly to maxinun-likelihood methods [35, 62]. In table we give the time taken train to convergence as we vary the dimension d. The time taken for GOLEM includes the time required to choose the sparsity parameter À via cross-validation. BCD Nets didn't need any cross-validation to choose sparsity parameters. All the methods were run on the same hardware, a single Nvidia 2080Ti GPU with 16 CPUs. Our method takes more time to converge than GOLEM. This not surprising, since we are training a neural network with many Sinkhorn iterations per optimization step. The Ghoshal algorithm [21] is very fast in comparison, consisting only of a precision estimation step then d stages of matrix multiplication. At high d, GADGET is slower than our method. We found that speeding up GADGET by reducing the number of iterations resulted in dramatically reduced performance.
5.5 Increasing Dataset Size
We expect that our method performs best in the low-data regime, where Bayesian approach with priors has an advantage over -an-asaei methods. To show this, we perform experiments with increased amounts of data, results of which are shown in figure 3. We see that the advantage of BCD Nets is reduced as the quantity of data used increased. We suspect that our variational method suffers from a more challenging optimization problem as the posterior becomes more peaked around the correct solution. Methods such as amortized inference or sequential Monte- Carlo may help to combat this problem and improve our model at larger n.

---

ER1, Caussion Noise
GOLEN
BCD Nets IDung
GADGET
Nail Buseline
GOSSEV
60
Chochel
28
200
403
600 82o
1000
Figure 3: Structural Hamming Distance to the ground truth as a function of n, on 64-dimensional, degree Erdös-Renyi graphs. With large amounts of data, baselines catch up BCD Nets.
5.6 Ablation
Finally we perform an ablation of our algo- rithm to illustrate the importance of architec- tural choices. We use the d = 32 case, average degree I, with Gaussian likelihood. We report the change in performance with a factored pos- terior (Mean-Field) and with Laplace prior on the weights (Laplace) instead of Horse- shoe prior. We also report performance with a fixed 100 Sinkhorn steps (Sinkhorn-100), as opposed the adaptive number that between the sample covariance and the covariance induced by the sampled parameters (E, T., P). which incorporates how well the posterior approximates the distribution of I., unlike SHD. We give the results in table and 4, in the appendix. Changing any of the parts of the algorithm result in SHD increasing from 11 to around 30. Interestingly, under the Laplace prior the SHD high but the KL-divergence is only modestly higher than with the Horseshoe. This indicates a learned DAG which can generate the data, but with extra edges due to an ineffective sparsity prior.
Table 2: The computing time required, in minutes,
converge a solution for several approaches.
16
32
64
128
GOLEM
25
30
40
65
90
GADGET GHOSHAL
40
150
385
635
1200
<1
<1
<1
3
15
BCD Nets (Ours)
50
160
300
350
900
in our approach. We also report the KL--divergence
Conclusion
We introduce Bayesian Causal Discovery Nets, a framework for performing Bayesian causal discov- ery for linear-Gaussian structural equation models. BCD Nets exploit recent advances in variational inference to flexibly model the posterior distribution over SEM parameters given data, outperform- ing methods that only return point estimates. An explicit posterior is also useful in the high-stakes, low-data regimes where causal inference is increasingly used [57, 47].
On the other hand, while indications are that our method may be robust to a small amount of mis- specification, the validity of the linear modelling assumption should be carefully considered in plications. Assuming a linear relationship is present where highly non-linear one exists could lead to harmfully incorrect inferences, particularly for minority groups with complex and under- studied variable interactions. Furthermore, our approach assumes that no unobserved confounders are present: variables which influence the observed variables but are not themselves observed. Since it is difficult to obtain useful inferential results under the effects of possibly arbitrary unknown con- founders, the assumption of no unobserved confounders is standard in many areas of causal infer- ence [38]. However, this assumption is unlikely hold exactly real applications, and so the possible influence of unobserved confounders must be seriously considered and their influence re- duced as much as possible before using our method. Future work could explore improvements to inference speed by replacing the Sinkhorn with more efficient algorithms from optimal transport.
10

---

Acknowledgements
This research supported by NSF(#1651565,, #1522054, #1733686), ONR (N000141912145) AFOSR (FA95501910024), ARO (1911N9-21-1-0122)) and the Sloan Fellowship. We thank the anonymous reviewers for their helpful comments.
Funding Transparency Statement
1. Funding: This research was supported by NSF(#1651565, #1522054, #1 #1733686), ONR (N000141912145), AFOSR (FA95501910024), ARO (W911NF-21-1-0122) and the Sloan Fellow- ship
2. Competing Interests: There competing financial interests that could be perceived to influence the contents of this work.

---

References
[1] Hirotogu Akaike. Information theory and an extension of the maximum likelihood principle. In Selected Papers of Hirotugu Akaike, pages 199-213. Springer, 1998. 7
[2] Hirotugu Akaike. Fitting autoregressive models for prediction. Annals of the institute of Sta- tistical Mathematics, 21(1):243-22 1969. 7
[3] Nima Anari and Alireza Rezaei. A tight analysis bethe approximation for permanent. In FOCS, 2019. 5
[4] Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating propagating gradients through stochastic neurons for conditional computation. arXiv, 2013. 6, 20
[5] Tolga Birdal and Umut Simsekli. Probabilistic permutation synchronization using the rieman- nian structure of the birkhoff polytope. In CVPR, pages 11105-11116, 2019.
[6] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. statisticians, JASA, 2017. 2
Variational inference: A review for
[7] Carlos M Carvalho, Nicholas Polson, and James G Scott. Handling sparsity via the horse- shoe. In AISTATS, 2009. 2, 6
[8] David Maxwell Chickering. Learning bayesian networks is NP-complete. In Learning from Data. Springer, 1996.
[9] Chow and Cong Liu. Approximating discrete probability distributions with dependence trees. IEEE transactions Information Theory, 14(3):462-467, 1968.
[10] Rob Cornish, Anthony Caterini, George Deligiannidis, and Arnaud Doucet. Relaxing bijectiv- ity constraints with continuously indexed normalising flows. ICML, 2020.
[11] James Cussens. Bayesian network learning with cutting planes.
arXiv preprint
arXiv: /202.3713 2012.
[12] James Cussens, Matti Järvisalo, Janne H. Korhonen, and Mark Bartlett. Bayesian Network Structure Learning with Integer Programming: Polytopes, Facets and Complexity. JAIR, Jan- uary 2017. 7
[13] Marco F Eigenmann, Prectam Nandy, and Marloes H Maathuis. Structure learning of linear Gaussian structural equation models with weak edges. arXiv preprint 1rXiv:1707.075600 2017.
2
[14] Paul Erdös and Alfréd Rényi. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5716,17-60, 1960.
[15] Ronald Fisher. Lung cancer and cigarettes? Nature, 182(4628);108-108, 1958.
[16] Tiago Fragoso, Wesley Bertoli, and Francisco Louzada. Bayesian model averaging: A systematic review and conceptual classification. International Statistical Review, 2018. 2
[17] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432-44) 2008.
[18] Nir Friedman and Daphne Koller. Being Bayesian about network structure. A Bayesian ap- proach to structure discovery in Bayesian networks. Machine learning, 2003. I
[19] Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. Bayesian Data Analysis, Third Edition. CRC Press, 2013. 2,
[20] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle. Made: Masked autoen- coder for distribution estimation. In ICML, 2015.7 7
[21] Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In International Conference on Artificial Intelligence and Statis- fics, pages 1466-1475. PMLR, March 2018. 8,
12

---

[22] Shengtong Han, Hongmei Zhang, Ramin Homayouni, and Wilfried Karmaus. An efficient bayesian approach for gaussian bayesian network structure learning. Communications in Statisties-Simulatiion and Computation, 2017. 2
[23] David Heckerman, Christopher Meek, and Gregory Cooper. A bayesian approach to causal discovery. Computation, causation, and discovery, 19:141-166, 1999. 2
[24] Tommi S Jaakkola and Michael Jordan. Bayesian parameter estimation via variational meth- ods. Statistics and Computing. 10(1):25-31 -7, 2000.
[25] Ajay Jain, Pieter Abbeel, and Deepak Pathak. Locally Masked Convolution for Autoregressive Models. In UAI, 2020.
[26] Daphne Koller and Nir Friedman. Probabilistic Graphical Models: Principles and Techniques. MIT press, 2009.
[27] Harold W Kuhn. The Hungarian method for the assignment problem. Naval research logistics quarterly, 2(1-2):83-97, 1955. 6, 20
[28] Xuanlin Li, Brandon Trabucco, Dong Huk Park, Yang Gao, Michael Luo, Sheng Shen, and Trevor Darrell. Discovering autoregressive orderings with variational inference. In ICLR, 2021. 2, 5, 2.
(29] Roderick J Little and Donald B Rubin. Causal effects in clinical and epidemiological studies via potential outcomes: Concepts and analytical approaches. Annual review of public health, 21(1):121-145, 2000.
[30] Brandon Malone, Kustaa Kangas, Matti Järvisalo, Mikko Koivisto, and Petri Myllymäki. As- asl: Algorithm selection with auto-sklearn. In Open Algorithm Selection Challenge 2017. PMLR, 2017.
[31] Gonzalo Mena, David Belanger, Scott Linderman, and Jasper Snock. Learning Latent Permu- tations with Gumbel-Sinkhom Networks. In ICLR, 2018. 5,5, 20
[32] Gonzalo Mena, Erdem Varol, Amin Nejatbakhsh, Eviatar Yemini, and Liam Paninski. Sinkhorn permutation variational marginal inference. AABI, 2020. 5
[33] Edgar C Merkle and Yves Rosseel. Blavaan: Bayesian structural equation models via param- eter expansion. JOURNAL OF STATISTICAL SOFTWARE, 85(4):1-30, 2018. 7
[34] Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih. Monte Carlo Gradient Estimation in Machine Learning. arXiv:/906e /00,:/2,00/0055,, June 2019. 5
[35] Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the Role of Sparsity and DAG Con- straints for Learning Linear DAGs. NeurIPS, 2020. 1, 2, 3, 6, 7, 8, 9, 19. 23
[36] Siqi Nie, Cassio de Campos, and Qiang Ji. Learning Bayesian networks with bounded tree- width via guided search. In AAAI, 2016. 7
[37] Judea Pearl. Causality, Cambridge university press, 2009. 1,2
[38] Jonas Peters and Peter Bühlmann. Identifiabilityy of Gaussian structural equation models with equal error variances. Biometrika, 101(1):219-22) 2014. 2, 3, 8, 10, 24
[39] Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. JMLR, 2014.
[40] Juho Piironen and Aki Vehtari. On the hyperprior choice for the global shrinkage parameter in the horseshoe prior. In AISTATS, 2017. 6
[41] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In ICML, 2015.
[42] Thomas Richardson, Peter Spirtes, et al. Ancestral graph markov models. The Annals of Statistics, 2002. 7
13

---

[43] Geoffrey Roeder, Yuhuai Wu, and David K Duvenaud. Sticking the landing: Simple, lower- variance gradient estimators for variational inference. In NIPS, 2017.
[44] Karen Sachs, Omar Perez, Dana Pe'er, er, Douglas A. Lauffenburger, and Garry P. Nolan, Causal protein-signaling networks derived from multiparameter single-cell data. Science (New York, N.Y.), April 2005. 7.9
[45] Mauro Scanagatta, Cassio Polpo de Campos, Giorgio Corani, and Marco Zaffalon. Learning bayesian networks with thousands of variables. In NIPS, 2015. 7,7
[46] Gideon Schwarz et al. Estimating the dimension of model. Annals of statistics, 6(2):461-46) 1978. 6, 7, 24
[47] Uri Shalit, Fredrik D. Johansson, and David Sontag. Estimating individual treatment effect: Generalization bounds and algorithms. arXiv: 03966 [cs, stat], June 2016. 1, 10
[48] Nuala A Sheehan, Mark Bartlett, and James Cussens. Improved maximum likelihood recon- struction of complex multi-generational pedigrees. Theoretical population biology, 97:11-19, 2014.7
[49] Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O Hoyer, and Kenneth Bollen. Directlingam: A direct method for learning linear non- con-gaussian structural equation model. The Journal of Machine Learning Research, 12:1 1248, 2011. 8
[50] Richard Sinkhorn. A relationship between arbitrary positive matrices and doubly stochastic matrices. The annals of mathematical statistics, 35(2):876-879, 1964. 5
[51] Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs Social science computer review, 1991.
[52] Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. Springer-Verlag, 1993. 7
[53] Mare Teyssier and Daphne Koller. Ordering-based search: A simple and effective algorithm for learning Bayesian networks. In UAI, 2005. 1,
[54] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodologicul), 58(1):267-288, 1996. 6
[55] Jakub M Tomezak. On some properties of the low-dimensional gumbel perturbations in the perturb-and-man model. Statistics & Probability Letters, 115:8-15, 15, 2016. 5
[56] loannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006. 7
[57] Hal R Varian. Causal inference in economics and marketing. PNaS, 113, 2016. 1, 10
[58] Jussi Viinikka, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. Towards scalable bayesian learning of causal dags. NeurlPS, 33, 2020. 2, 8,8.9 9
[59] Dennis Wei, Tian Gao, and Yue Yu. DAGs with no fears: A closer look at continuous opti- mization for learning bayesian networks. NewrlPS, 2020. 1
[60] Ronald Williams. Simple statistical gradient-following algorithms for connectionist rein- forcement learning. Machine learning, 8(3-4):229 256, 1992. 5,7
[61] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG Structure Learning with Graph Neural Networks. arXiv: 1904. 10098 fes, stat), April 2019. 1,7
[62] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. Dags with no tears: Con- tinuous optimization for structure learning. In NeurlPS, 2018 1,3, 6, 7, 8,
[63] Juntang Zhuang, Tommy Tang, Sekhar Tatikonda, Nicha Dvornek, Yifan Ding, Xenophon Papademetris, and James S Duncan. Adabelief optimizer: Adapting stepsizes by the belief in observed gradients. In NeurlPS, 2020. 19
14

---

Checklist
15

---

1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] Yes, the conclusion we discuss how we are limited to linear relationships. We also emphasise throughout the paper that in many settings, the true DAG parameters are nonidentifiable from observational data
(c) Did you discuss any potential negative societal impacts of your work? [Yes] In the conclusion we discuss how wrongly assuming linear relationships could lead to incor- rect inferences, and how those could be amplified in minority groups
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] We conform to the ethics review guidelines
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] We give the full preconditions for theorem to hold in the supplementary material
(b) Did you include complete proofs of all theoretical results? [Yes] We give full proof of theorem in the supplementary material
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed reproduce the main exper- imental results (either in the supplemental material or as a URL)? [No] Given time constraints, we were not able to get the code shape to release
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] In the additional material discuss all hyperparameters, in- cluding for the valss-validation for GOLEM
(c) Did you report error bars (e.g., with respect to the random seed after running exper- iments multiple times)? [Yes] Yes, all results have error bars from running with several random seeds
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] This is described in the appendix, section I
4. If you are using existing assets (e.g., code, data, models) or curating/relcasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] Yes, we used the data producing code from Ng al and they are cited
(b) Did you mention the license of the assets? [Yes] We described the license
(c) Did you include any new assets either in the supplemental material or as a URL? [No]
(d) Did you discuss whether and how consent obtained from people whose data you're using/curating? [N/A] Under the license there is no need to get consent from the github maintainers long as we properly cite
(e) Did you discuss whether the data you are using/curating contains personally identifi- able information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
16


Reviews:
[{"review_id": "vm91MJ_1lZE", "reviewer": "Reviewer_uT5D", "summary": "### 1. Summary\n\nThe paper proposes a novel framework in the field of (causal) structure learning, using a variational inference method.\n\nThis framework is proposed to be able to provide estimates of distribtutions in the graph regime of the DAG to be learned.\n", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### 2. Rationale for the score\n\nThe problem chosen is a well researched problem. Showing that indeed state of the art methods can be outperformed is remarkable.\n\nThe use of synthetic and real world data for experiments corresponds with the usual approach in this field of research.\n\nThe paper is well structured, graphical contributions in the main part as well as in the appendix serve to enable a good understanding of the approach done.\n\nThe novelty of providing distributions instead of point estimates is a good contribution in this field of research.\n\n### 3. Positive aspects\n\n- Using the Structural Hamming Distance is a good choice here to be able to compare true and learned DAGs\n\n- The ablations performed are well chosen and provide a good benefit to the experiments done\n\n- The CO2 impact mentioned in the appendix is a nice detail\n\n- The \"causal inference\" section is a nicely chosen addition to the structure learning parts.\n\n\n### 4. negative aspects\n\n- I would have expected a codebase provided in the supplementing material for the implementation of the algorithm\n\n- More discussion of the results and outlooks regarding special corner cases (e.g. hidden confounders (=unobserved variables) in the true graphs) would be a nice addition", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I would suggest an additional paragraph about potential negative societal impact", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "JBxz4ZviJQU", "reviewer": "Reviewer_nWpC", "summary": "The paper proposes a variational inference approach for approximate Bayesian causal discovery. In that, it focuses on modelling distribution over linear-Gaussian (LG) SEMs. The LG-SEMs are represented as tuples of an observation noise matrix  $\\Sigma$ and a weighted adjacency matrix $W = PLP^\\top$  and, which are modelled as combination of strictly lower-triangular real-valued matrix $L$ followed by a variable permutation represented by permutation matrix $P$. The proposed algorithm maximizes the ELBO w.r.t. to the parameters $\\phi$ of the variational posterior $q_\\phi(P|L, \\Sigma) = q_\\phi(P|L, \\Sigma) q_\\phi(L, \\Sigma)$. Here, $q_\\phi(P|L, \\Sigma)$ is modelled as conditional Bolzmann distribution where re-parametrized sampling of permutation matrices $P$ is facilitated via the Gumble-Sinkhorn relaxation - this way pathwise gradient estimators of the ELBO can be obtained. In the experiment section, the proposed method is compared to three related approaches on simulated Erdos-Renyi graphs and on the well-known protein-signalling task.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Writing:** The paper is well-written, and very easy to follow. Though, some parts of the paper are a bit verbose, particularly Section 2, 3.3 and 5. Making them more concise would free up space that could be used to e.g. strengthen the experiment section.\n\n**Relevance:** The paper addresses a relevant problem of facilitating Bayesian inference over SCMs. The key motivation for a Bayesian treatment of causal discovery is to be able to quantify the epistemic uncertainty over SCMs - in principle the proposed method is able to do so.\n\nUnfortunately, the paper hardly evaluates the uncertainty quantification aspect of the proposed method. The TP and FP rate are displayed in the appendix. Why not report the AUROC as done by most papers in the field of Bayesian Causal Discovery? Then the information is compressed into a single number and you can present the results in the main paper. \n\nAlso Fig. 5 in the appendix is a nice result - but why doesn’t it include GADGET and DLiNGAM? I recommend adding the corresponding results for the other two methods and moving the result to the main paper. Alternatively, you could report the interventional log-likelihood as done in [3] - this nicely evaluates the methods towards the purpose of reasoning about interventions under uncertainty - the main purpose of Bayesian causal discovery.\n\n**Originality & Significance:** None of the technical pieces in the paper is particular novel / original. However, the paper skilfully assembles recent method such as the  re-parameterized sampling of permutation matrices via the Gumbel-Sinkhorn distribution, the approximation of the matrix permanent, the horseshoe prior to enforce sparsity etc., to develop a tractable variational inference algorithm over SCMs. Thus, the overall algorithm can certainly be considered a contribution to the research field. A particular advantage of the method is that it does not need to explicitly enforce a acyclicity constraint like the NOTEARS method [1].\n\n**Quality & soundness:** The proposed method is sound and well-motivated. Unfortunately, the experiment section is a bit weak. Two of the three baseline methods do not even model a distribution over DAGs - both output a point estimate. The only really relevant baseline is GADGET. It would be good to include simple baselines that can capture epistemic uncertainty such as bootstrapped versions of standard algorithms such as LiNGAM, GES or NOTEARS (e.g. see [2]). Moreover, the ER-1 graphs used in the experiment section are very sparse. It would be good to include an experiment with harder problems such as ER-3 or ER-4 where predicting the null graph is not already a strong baseline.\n\nA major concern of mine is the lack of code which makes the method & experiment results very hard to reproduce. Since this paper is a resubmission and there is more than a week after the deadline to prepare and polish the supplementary material, failing to provide usable and well-documented code is definitely a red flag. \n\n**Overall Assessment:** All in all, the method proposed in the paper is relevant, sound, well-presented and certainly of interest to the causal inference community - not a ‘game changer’ but a solid contribution.  Unfortunately, the paper falls short of a convincing empirical evaluation and fails to provide code. In the current form I do not see the paper above the acceptance threshold for NeurIPs. However, if the authors sufficiently address the mentioned shortcomings I am happy to change my assessment.\n\n[1] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. 2018. “DAGs with NO TEARS: continuous optimization for structure learning”, NeurIPS, 2018.\n[2] Agrawal, Raj, et al. \"Abcd-strategy: Budgeted experimental design for targeted causal structure discovery.\", AISTATS, 2019.\n[3] Lorch, Lars, et al. \"DiBS: Differentiable Bayesian Structure Learning.\" arXiv preprint arXiv:2105.11839 (2021).\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper discusses the limitations of the linear Gaussian SEM assumption. It would be good to also discuss how restrictive the chosen variational family $\\{q_\\phi(P, L, \\Sigma): \\phi \\in \\Phi\\}$ is.  \n\nSince the paper concerns algorithmic aspects of approximate Bayesian causal discovery, it is generally applicable in many fields - thus, I do no expect any immediate societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4-5 hours"}, "extra_text": {}}, {"review_id": "HfzoOAsaEYb", "reviewer": "Reviewer_MkjZ", "summary": "This submission presents an approach to quantify the uncertainty in estimating linear DAG structure using variational inference. In particular, the authors present a set of modeling design choices that (i) lead to tractable gradients of the ELBO, and (ii) are robust to certain classes of model misspecification. They go on to demonstrate that their approach outperforms maximum likelihood-based approaches on synthetic and semi-synthetic data.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall this submission addresses an important and often overlooked task of quantifying uncertainty in structure learning. As the authors note, quantifying uncertainty in this way is particularly important in settings where the number of observations are limited, or when structure can not be uniquely identified from observational data asymptotically. In addition to the problem being relevant, the design decisions are reasonable and well motivated. \n\nHowever, I believe the submission could be significantly improved by (i) expanding the baselines and metrics in the empirical evaluation and (ii) clarifying the necessary assumptions relating the proposed approach to causality.\n\nEmpirical evaluation:\n\nAs the authors note, probabilistic approaches to structure learning have many benefits over maximum likelihood approaches. However, the authors do not evaluate their claim that variational approaches provide scalability benefits over e.g. sampling-based approaches to probabilistic structure learning. As one of several example baselines to compare against, Mohammadi and Wit, 2015 propose a birth-death MCMC algorithm for a very similar (or identical) problem formulation of linear Gaussian DAGs. In addition to this (potentially strong) baseline, it would be good to compare against a more naive sampling-based approach such as a random walk MCMC over graph permutations and weights.\n\nIn addition to new baselines, the experiments could be improved by evaluating the performance at various observational sample sizes for identifiable graphs. My expectation is that we should see the maximum likelihood approaches and the variational approaches converge to the same performance as n gets very large. This would emphasize the domain under which the proposed method is an appropriate choice, which I believe is large d and small n.\n\nDespite these critiques, I am very happy to see that the authors performed an ablation study of their proposed method. This honest empirical assessment is a very good trend for the machine learning research community.\n\nClarity:\n\nAs the authors claim that their method can \"infer causal relationships\" it is extremely important to specify the necessary assumptions for this interpretation upfront. As an example, are the authors assuming causal sufficiency (i.e. no latent confounders)? \n\nMinor critique:\n\nThere are many spelling and grammatical errors throughout the draft that negatively impact its quality. I highly recommend taking a thorough editorial pass through future versions of the submission. This has not affected my recommendation.\n\nA. Mohammadi, E. C. Wit. Bayesian Structure Learning in Sparse Gaussian Graphical Models. Bayesian Anal. 2015\n\nUpdate: I have increased my score after seeing the authors' response.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No, these authors have not considered limitations and potential negative societal impact of their work. That being said, the methodological focus of this submission means that this is no more a concern than any other work.", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "6AecmwjqqnK", "reviewer": "Reviewer_PQX5", "summary": "This paper introduces a causal discovery method for linear Gaussian SEM. The authors extend the variational inference to learn the causal structure for linear SEM data.", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper considers a fundamental problem of learning causal relationships in the linear model with gaussian distributions.  \nThe paper builds on prior work in the fields of variational inference and bayesian estimation. \\\nThe authors put some effort into learning the posterior distribution $p(W,\\Sigma|X)$ by variational inference.\\\nThe paper provides the details of the critical distribution choices for the estimation of the parameters. \\\nThe experiments are quite thorough, and the paper is well-organized.\n \nHowever, there are some concerns or suggestions:\n\n -It is unfair to compare with the DLiNGAM. This method was designed for the non-Gaussian dataset. The authors may need to compare with methods that are designed for linear Gaussian data, such as [1,2]\n \n1. J. Peters and P. Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1): 219--228, 2014.\n \n2. Ghoshal A, Honorio J. Learning linear structural equation models in polynomial time and sample complexity. International Conference on Artificial Intelligence and Statistics. PMLR, 2018: 1466-1475.\nBesides,\n \n-The algorithm is evaluated on synthetic cases with only 100 datapoints. It would have been more convincing with more data points. \n \n-Why the PC and GES algorithms are not presented in the synthetic data? What is the reason?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This article introduces an approach to causal graph inference using variational Bayesian inference.  The paper focuses on linear-Gaussian structural equation models, and employs a carefully designed variational approximation to the posterior on directed acyclic graphs (DAGs) to handle the challenging combinatorial nature of the graph space; the proposed approach is called LiGa-VI.  The posterior is parametrized in terms of $(P,L,\Sigma)$ where $P$ is a permutation matrix for the vertex ordering, $L$ is a lower triangular matrix such that $P L P^\texttt{T}$ is the graph adjacency matrix, and $\Sigma$ is a diagonal covariance matrix for the Gaussian model.  For the variational approximation $q_\phi(P,L,\Sigma)$, the authors use a normal distribution or normalizing flow for $q_\phi(L,\Sigma)$ and a Gumbel-Sinkhorn relaxation of the Gumbel-Matching distribution over permutations for $q_\phi(P|L,\Sigma)$. Gradient-based optimization is used to optimize the variational parameter $\phi$ to fit the posterior by minimizing KL divergence.  Experiments are performed to assess performance versus several competing methods on synthetic and real data.

The variational approximation used over DAG space is interesting and innovative.  The experimental performance appears to be surprisingly good in terms of concentration near the true graph.  Overall, the reviewers were quite positive about the paper, but I have some significant concerns.

1) No improvement in performance as the sample size $n$ increases? I'm puzzled by the performance as a function of the sample size $n$ in the additional figure supplied by the authors in their reply to Reviewer PQX5 (https://www.dropbox.com/s/clzi1tzp0jgdwr0/new_fig_n.png?dl=0).  Why does the proposed method (LiGa-VI) not benefit from having more data?  Similarly, why does the Bayesian method GADGET not benefit from having more data after $n$ increases beyond $300$?  Other methods such as Ghosal's are clearly benefiting from having more data.  This makes me concerned that the good performance of LiGa-VI here is solely due to something trivial like the choice of prior.  This needs to be explained.  Also, I highly recommend adding this plot to the paper.

2) I'm skeptical about Figure 2.  To have such a big difference in performance, it seems like something about this simulation is cherry picked to yield good performance.  Any Bayesian method with the same prior should yield roughly the same accuracy here, up to the accuracy of the posterior approximation.  Is there a big difference in the priors used by GADGET and LiGa-VI?  Or is GADGET not sampling well from the posterior?  It is fine to present simulations that showcase the proposed method, but they need to be put in context to show the relative strengths *and* weaknesses of the method.

3) The computation time of LiGa-VI is not that attractive (Appendix H, Table 4), taking around 2-10x the time required for GOLEM even though GOLEM is using cross-validation while LiGa-VI is not.  I found this disappointing since I would have expected a variational approach to be more computationally attractive.  I would highly recommend including the computation time for the other competing approaches, especially GADGET since it employs MCMC for Bayesian inference.

4) Mathematical notation and unclear exposition. Some of the notation and descriptions were not at all clear, such as "$q_\phi(L,\sigma) = \mathcal{N}(\mu_\phi,\sigma_\phi)$" for the variational approximation to the posterior on $(L,\Sigma)$ --- the left-hand side is a density while the right-hand side is a distribution, where does $L$ enter into the right-hand side, and how does this yield a distribution over $\sigma$?  Please make sure all mathematical expressions are properly written for clear exposition. I was also confused by this phrase in Appendix H: "we are training a neural network with many Sinkhorn iterations per optimization step", since I didn't see where LiGa-VI is using a neural network in the description of the method.  This needs to be clarified.


Decision: Accept (Poster)