DiBS: Differentiable Bayesian Structure Learning
Authors: Lars Lorch
Jonas Rothfuss
Bernhard Schölkopf
Andreas Krause
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: causality
structure learning
Bayesian networks
Bayesian inference

Abstract: Bayesian structure learning allows inferring Bayesian network structure from data while reasoning about the epistemic uncertainty---a key element towards enabling active causal discovery and designing interventions in real world systems. In this work, we propose a general, fully differentiable framework for Bayesian structure learning (DiBS) that operates in the continuous space of a latent probabilistic graph representation. Contrary to existing work, DiBS is agnostic to the form of the local conditional distributions and allows for joint posterior inference of both the graph structure and the conditional distribution parameters. This makes our formulation directly applicable to posterior inference of nonstandard Bayesian network models, e.g., with nonlinear dependencies encoded by neural networks. Using DiBS, we devise an efficient, general purpose variational inference method for approximating distributions over structural models. In evaluations on simulated and real-world data, our method significantly outperforms related approaches to joint posterior inference.



DiBS: Differentiable Bayesian Structure Learning
Lars Lorch ETH Zurich
Jonas Rothfuss
ETH Zurich Zurich, Switzerland j rothfusseiff wwtpassinnffeehhh..cmmmm
Zurich, Switzerland
lars. wor.chrrnniff. .wt...c.t.cccccc
Bernhard Schölkopf
Andreas Krause ETH Zurich Zurich, Switzerland krauseaßethz.ch ch
MPI for Intelligent Systems
Tübingen, Germany
bs@tuebingen.c mpg.mp...... de
Abstract
Bayesian structure learning allows inferring Bayesian network structure from data while reasoning about the epistemic uncertainty-4 key element towards enabling active causal discovery and designing interventions in real world systems. In this work, we propose a general, fully differentiable framework for Bayesian structure learning (DiBS) that operates in the continuous space of latent probabilistic graph representation. Contrary to existing work, DiBS is agnostic to the form of the local conditional distributions and allows for joint posterior inference of both the graph structure and the conditional distribution parameters. This makes our formulation directly applicable posterior inference of complex Bayesian network models, e.g., with nonlinear dependencies encoded by neural networks. Using DiBS, we devise an efficient, general purpose variational inference method for approximating distributions over structural models. In evaluations on simulated and real-world data, our method significantly outperforms related approaches to joint posterior inference.'
Introduction
Discovering the statistical and causal dependencies that underlie the variables of Mlaa-ggenerttin system is of central scientific interest. Bayesian networks (BNs) [1] and structural equation models are commonly used for this purpose [2-5]. Structure learning, the task of learning BN from observations of its variables, well-studied, but computationally very challenging due to the combinatorially large number of candidate graphs and the constraint of graph acyelicity.
While structure learning methods arrive at single plausible graph or its Markov equivalence class (MEC), e.g., [6-9], Buyesian structure learning aims to infer a full posterior distribution over BNs given the observations. distribution over structures allows quantifying the epistemic uncertainty and the degree of confidence in any given BN model, e.g., when the amount of data is small. Most importantly, downstream tasks such as experimental design and active causal discovery rely on posterior distribution over BNs to quantify the information gain from specific interventions and uncover the causal structure a small number experiments [10-15].
A key challenge in Bayesian structure learning is working with posterior over BNs a distribution over the joint space of (discrete) directed acyclic graphs and (continuous) conditional distribution parameters. Most of the practically viable approaches to Bayesian structure learning revolve around
'Our Python JAX implementation of DiBS is available https: / /tt///ggttttttt http::::///gttuu.mmm/////r////// www.:/assloorrc//......
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Markov chain Monte Carlo (MCMC) sampling in combinatorial spaces and bootstrapping of classical score and constraint- based structure learning methods, e.g., in causal discovery [10-15]. However, these methods marginalize out the parameters and thus require closed form for the marginal likelihood of the observations given the graph to remain tractable. This limits inference to simple and by now well-studied linear Gaussian and categorical BN models [16-18] and makes it difficult to infer more expressive BNs that, c.g., model nonlinear relationships among the variables. Due to the discrete nature of these approaches, recent advances in approximate inference and gradient-based optimization could not yet be translated into similar performance improvements in Bayesian structure learning.
In this work, we propose novel, fully differentiable framework for Bayesian structure learning (DiBS) that operates in the continuous space of a latent probabilistic graph representation. Contrary to existing work, our formulation is agnostic to the distributional form of the BN and allows for inference of the joint posterior over both the conditional distribution parameters and the graph structure. This makes our approach directly applicable to more flexible BN models where neither the marginal likelihood nor the maximum likelihood parameter estimate have closed form. We instantiate DiBS with the particle variational inference method of Liu and Wang [19] and present general purpose method for approximate Bayesian structure learning. In our experiments on synthetic and real-world data, DiBS outperforms all alternative approaches to joint posterior inference of graphs and parameters and when modeling nonlinear interactions among the variables, often by a significant margin. This allows to narrow down plausible causal graphs with greater precision and make better predictions under interventions an important stepping stone towards active causal discovery.
2 Background
Bayesian networks A Bayesian network (G, b) models the joint density p(x) of a set of d vari- ables x = Fi:d using (1) a directed acyclic graph (DAG) G encoding the conditional independencies of and (2) parameters @ defining the local conditional distributions of each variable given its parents in the DAG. When modeling p(x) using a BN, each variable is assumed to be independent of its non-descendants given its parents, thus allowing for compact factorization of the joint p(x G, O) into a product of local conditional distributions for each variable and its parents in G.
Bayesian inference of BNs Given independent observations D = -xx(), ......; x()),, we consider the task of inferring full posterior density over Bayesian networks that model the observations. Following Friedman and Koller [20], given a prior distribution over DAGs p(G) and a prior over BN parameters p( G), Bayes' Theorem yields the joint and marginal posterior distributions
p(G, θ (D) x p (G) G))(D)) O)
p D) x PG)))p(D
(2)
where p(D G) = ( G, d)e is the marginal likelihood. Thus, p(G D)) in (2) is only tractable in special conjugate cases where the integral over θ can be computed in closed form. The Bayesian formalism allows us to compute expectations of the form
E E D) [ O)
or E D) [[G]]]
(3)
for any function f of interest. For instance, to perform Bayesian model averaging, we would use J[G, θ) = p(x G, θ) or (G) = p(x G), respectively [21, 22). In active learning of causal BN structures, a commonly used is the expected decrease in entropy of G after intervention [10-12, 14, 15). Inferring either posterior is computationally challenging because there are (dd2(a) possible DAGS with d nodes [23]. Thus, computing the normalization constant p(D) is generally intractable.
Continuous characterization of acyelic graphs Orthogonal to the work on Bayesian inference, Zheng et al. [9] have recently proposed a differentiable characterization of acyclic graphs for structure learning. In this work, we adopt the formulation of Yu et al. [24], who show that a graph with adjacency matrix G E (0, 1] 1 does not have any cycles if and only if ((G) O, where
GGG) : tr (1 GG) dd
(4)
If ((G) > O, the function can be interpreted as quantifying the cyclicity or non-DAG-nest of G. Follow-up work has leveraged this insight to model nonlinear relationships [24-27], time-series data
[28), in the context of generative modeling [29, 30), for causal inference [31-33], and contributed
to its theoretical understanding [34, 35]. So far, connection to Bayesian structure learning has been missing.
2

---

3
Related Work
The existing literature on Bayesian Structure Learning predominantly focuses on inferring the marginal graph posterior G D). Since this requires p(2) G) to be tractable, inference is limited to BNs with linear Gaussian or Categorical conditional distributions [16-18]. By contrast, the formulation we introduce overcomes this fundamental restriction by allowing for joint inference of the graph and the parameters, thereby facilitating the active (causal) discovery of more expressive BNs.
MCMC Sampling from the posterior
graphs is the most general approach to approximate Eatropolli-Hassting in the
Bayesian structure learning, Structure MCMC (MC*) [36, 37] performs space of DAGs by changing one edge time. Several works try to remedy its poor mixing behavior [38-40]. Alternatively, order MCMC draws samples in the smaller but still exponential space of node orders, which typically requires a hard limit on the maximum parent set size [20]. Attempts correct for its unintended structural bias are themselves NP-hard to compute and/or limit the parent size [38, 41]. By performing variational inference circumvents such mixing issues and parent size limitations.
continuous latent space, the method we propose
Bootstrapping The nonparametrie DAG bootstrap [42] performs model averaging by bootstrapping
D, where each resampled data set is used to learn a single graph, e.g. using the GES or PC algorithms I6. 7]. The obtained set of DAGs approximates the posterior by weighting each unique graph by its unnormalized posterior probability. In simple cases, a closed-form maximum likelihood parameter estimate may be used to approximate the joint posterior [14], but only if p(D G) is tractable in the first place.
Exact methods A few notable exceptions use dynamic programming to achieve exact marginal inference in time (d2d), which is only feasible for d S 20 nodes [43, 44]. In special cases, e.g., for tree structures or known node orderings, exact inference be performed more efficiently [45, 46].
4
A Fully Differentiable Framework for Bayesian Structure Learning
4.1
General Approach
With the goal of moving beyond the restrictive conjugate setups required by most discrete sampling methods for Bayesian structure learning, we propose to transfer the posterior inference task into the latent space of a probabilistic graph representation. Our result- ing framework consistent with the original Bayesian structure learning task in (3), enforces the acyclicity of G via the latent space prior, and provides the score the continuous latent pos- terior, thus making general purpose inference methods applicable off-the-shelf..
Z
X
N
Generative model of BNs with latent variable Z. This formulation generalizes the standard Bayesian setup in (1) where only G. θ, and x are modeled explicitly.
Figure 1:
Without loss of generality, we assume that there exists latent variable Z that models the generative process of G. Specifically, the default generative model described in Section 2 is generalized into the following factorization:
(5) Figure 1 displays the corresponding graphical model. The following insight provides us with an equivalence between the expectation we ultimately want to approximate and an expectation over the posterior of the continuous variable Proposition 1 (Latent posterior expectation). Under the generative model in (5), it holds that
pZZ G. o, D) p( Z)p(O) ())(D|G, +)
Ep(G 12 [/(G)p(D) GG]
(a) B D) [[[]]] (y(Z)
and
Bp 2) [p(D] G)
Ep(ai [f(G, O)p G)p(D (AADD|GG))) G.O)]
(b) EnG, D) [[(,,) Epzz,e D)
E)(()) [p(0] G) G, e]
proof is provided in Appendix A.1. Rather than approximating P (D) or p(G. θ DD) our goal will be to infer p(Z (D) pZZ, θ (D) instead, which by Proposition 1 allows to compute
3

---

expectations of the form in (3). In the following, we first discuss how to define the two factors p(G Z) and p(Z) in way that models only directed acyclic graphs. Then, we provide the details necessary to perform black box inference of the posteriors of the continuous latent variable Z.
4.2 Representing DAGs in a Continuous Latent Space
Generative model of directed graphs We define the latent variable as consisting of two embedding matrices U, V E [ i... Z = [U, V.. Building on Kipf and Welling [47], we propose
to use bilinear generative model for the adjacency matrix G E {0, jddd of directed graphs using the inner product between the latent variables in Z:
d d
Po Z)
G
Po (gg (us,v)) vj) with (a 1 up vy) (u, vj)
(6)
where aa (x) 1/(1 exp(-oar) denotes the sigmoid function with inverse temperature a. We denote the corresponding matrix of edge probabilities in G given Z by GZ E [0, 1jdxd with
(Z) == Pa (gig
u1 Vz)
Since model acyclic graphs, which do not contain self-loops, we Pa (a 1 Z) = 0. The latent dimensionality k trades off the complexity of the variable interactions with tractability during inference. For A: 2 d. the matrix of edge probabilities (Z) is not constrained in rank, and the generative model in (6) can represent any adjacency matrix without self-loops. That said, the size of the latent representation Z only grows (dd-k) k) and, in principle, kt can be chosen independently of d. Note that the formulation in (6) models directed graphs since aa V2) # Ja( v:). We can even interpret u, and v: as node embeddings that may encode more information than mere edge probabilities, e.g., for graph neural networks. The fact that P. (G Z) is invariant to orthogonal transformations of U and V is not an issue in practice.
Acyclicity via the latent prior distribution A major constraint in learning BNs is the acyclicity of G. With a latent graph model p(G Z) in place, we design the prior p( to act as a soft constraint enforcing that only DAGs are modeled. Specifically, we define the prior of Z as the product of independent Gaussians with a Gibbs distribution that penalizes the expected cyclicity of G given Z:
Pa Z x exp (-3 EEGGI [MG]]) I[5( 0,
(8)
Here, h is the DAG constraint function given in (4) and $ is inverse temperature parameter controlling how strongly the acyclicity is enforced. As 3 + co, the support of ps (Z) reduces to all Z that only model valid DAGs. The Gaussian component ensures that the norm of Z is wellebehaved. Traditional graph priors of the form GG) x f() that induce, e.g., sparsity of G, can be flexibly incorporated into (Z) by means of an additional factor involving (G.. (Z)) Ep (G)) Unless highlighting specific point, we omit writing the hyperparameters a and /3 to simplify notation.
4.3 Estimators for Gradient-Basse -aased Bayesian Inference
The extended generative model in (5) not only allows us to incorporate the notoriously difficult acyclicity constraint into the Bayesian framework. By rephrasing the posteriors p(G (D) and p(G, θ D) in terms of p(Z D) and pZZ, + D), respectively, also makes Bayesian structure learning amenable to variational inference techniques that operate in continuous space and rely on the gradient of the unnormalized log posterior, also known as the score. The following result provides us with the score of both latent posteriors. Detailed derivations can be found in Appendix A.2.
Proposition 2 (Latent posterior score). Under the generative graph model defined in (5), the gradient of the log posterior density of Z is given by
Vz EpIG [PD G
(a)
Vz loggopZZ D) Vz log
(9)
ExG|z) p(D [G]
which is relevant when the marginal likelihood p(D G) can be computed efficiently. In the general case, when inferring the joint posterior of Z and θ the gradients of the log posterior are given by
Vz E z) ( e D G)]
(b) Vzlogp(Z. θ D) = Vz lZ) +
(10)

---

E B Vep(O,D G)
(c)
Ve logp D)
(11)
By p D G)]
where p p D G) - p (D
G, h) which is efficient to compute by construction.
The expectations have tractable Monte Carlo approximations because sampling from G (Z) in (6) simple and parallelizable. The gradient terms of the form VzEng z)][]] I], which also appear inside Vz logp(Z); (Z), be estimated using two different techniques, depending the BN model inferred.
Differentiable (marginal) likelihood
Using the Gumbel-softmax trick [48, 49], we can separate the randomness from Z when sampling from p(G Z) and obtain the following estimator:
VzEmG 2) [P(D] G)] Ep)) V (GGPDDG))==....22 Vz G, (L, Z)
(12)
where L ~ Logistic(O, 1)dx i...d The matrix-valued function (.) is defined elementwise as
G (L, 2) = { (4g - ou? vz)
f i
(13)
The estimator applies equally to p D G) place of p(D G). For the estimator to be well- defined, p(D) G) or O D G), respectively, needs to be differentiable with respect to G. More specifically, the gradient Gp G) or Vcp(D, θ G) needs to be defined when lies on the interior of [0, 1j4xd and
its discrete endpoints. This depends on the parameterization of the BN model we want to infer. In case p(D G) or p( O D G) is only defined for discrete G, it is possible evaluate GG G) or V GGp(O, O D G) using hard Gumbel-max samples of G (i.e., with T oo) and a straight- through gradient estimator. Since the DAG constraint h is differentiable, the Gumbel-softma trick can always be applied inside Vz log P(Z). In practice, we always use T
Non-differentiable (marginal) likelihood In general, Vap(D G) or Vap(D, θ G) depending on the inference task might be available ill-defined. In this setting, the score function estimator provides us with a way to estimate the gradient we need [50]:
V E z) [P(D G)] PPp (M(DIG) b) Vz log Z))
(14)
The estimator likewise applies for p G) in place of p(D G). Here, b constant with respect to G that can be used for variance reduction [51], and Vz log p(G Z) is trivial to compute. The derivations of both (12) and (14), alongside more detailed discussion, can be found in Appendix B.
5 Particle Variational Inference for Structure Learning
In the previous section, we have proposed differentiable formulation for Bayesian structure learning (DiBS) that is agnostic to the form of the local BN conditionals and, more importantly, translates learning discrete graph structures into an inference problem over the continuous variable Z. In the following, we overcome the remaining challenge of inferring the intractable DiBS posteriors p(Z) D) and p(Z. ZZ θ D) by employing Stein variational gradient descent (SVGD) [19]. gradient- based and general purpose variational inference method. The resulting algorithm infers a particle approximation of the marginal or joint posterior density over BNs given observational data.
SVGD for posterior inference Since Proposition 2 provides us with the gradient of the latent poste- rior score functions, we can apply SVGD off-the-shell SVGD minimizes the KL divergence to a target distribution by iteratively transporting a set of particles using a sequence of kernel-based transforma- tion steps. We provide more detailed overview of SVGD in Appendix C. Following this paradigm for DiBS, we iteratively update a fixed set (m MM or {Z(m) θ (m) mm=1 to approximate p(Zi [Z] D) or pZZ θ D), respectively. If the BN model we aim to infer has properly-definet likelihood gradient with respect to G, we use the Gumbel- estimator in (12) to approximate the posterior score. Otherwise, we resort to the score function estimator in (14). We use a simple kernel for SVGD:
M
1 Z
( (Z,
cxp
+ exp
||e 0'
(15)
7=
with bandwidths 72, 70. For inference of p(Z D), we leave out the second term involving O O'. While Z is invariant to orthogonal transformations, more elaborate kernels that are, e.g., invariant to such transforms empirically perform worse in experiments.
5

---

Algorithm 1 DiBS with SVGD [19] for inference of p(G G
Input: Initial of latent particles (Z(")) J.M
Jm=1. kernel k, schedules for cu,, Bt, and stepsizes y/: Output: Set of discrete graph particles (G!m)) Jmm1 approximating p(G D))
1: Incorporate prior belief of p into (Z)
D. See Section 4.2
2: for iteration i = 0 to T 1 do
3: Estimate score Vz log p ()) given (9) for each z!-))
See (12) and (14)
particle m
M
5:
Z
Z (m + 1/ Φ (m)
SVGD step
where xx(+)
1
(8) k(Z
(D) + V V k(Z V Vz("0 logp(Z,")
M
6: return oGoo (Z (m)
)
b See (16) and (17)
Annealing a and | The latent variable Z not only probabilistically models the graph G, but can also be viewed as a continuous relaxation of G. with a trading off smoothness with accuracy. As + oC, the sigmoid (+) converges to the unit step function. Hence, as a - oo in the graph model Pa (G ZZ) in (6), the expectations in Proposition simplify to:
BpG D)
f(G)
+ pptz D) [())(())]
(16)
ExG.e f(G, + Epzz,1))) (( 6)]
where Gx denotes the single limiting graph implied by Z = UU Vj and defined as
if u > 0 and otherwise G. {
(17)
See Appendix A.3. In this case, ( Z) converges to representing only a single graph. To be able to invoke this simplification, we anneal a - oo over the iterations of SVGD and, upon termination, convert the latent variables Z to the single discrete G (Z). Furthermore, we similarly let B - co in the latent prior pa (Z) over the iterations to enforce that the latent representation of G only models DAGs. By Equation (16), the resulting DAGs form consistent particle approximation of p( (D) or pl G. D), respectively. Algorithm summarizes DiBS instantiated with SVGD for inference of p G D). The general case of inferring P(G, θ ( is given in Algorithm 2 of Appendix D.
D).
Single-particle approximation SVGD reduces to regular gradient ascent for the maximum a posteriori estimate when transporting only single particle [19]. In this special case, DiBS with SVGD recovers some the existing continuous structure learning methods: gradient ascent on a linear Gaussian likelihood solves optimization problem similar to NOTEARS [9]. The cyclicity penalizer acts analogously. However, not only does DiBS automatically turn into a full Bayesian approach when using more particles, it is also not limited to settings such as linear Gaussian conditionals, where the adjacency matrix and the parameters O can be modeled together by a weighted adjacency matrix.
Weighted particle mixture In high dimensional settings. it may be beneficial to move beyond uniform weighting of the inferred particles of BN models to approximate p( (D) or p(G. 9 DD). We consider a particle mixture that weights each particle by its unnormalized posterior probability P(G. D) or p( θ D), respectively, under the BN model. While we do not have a strong theoretical justification for the weighting, will unique as a consequence of the super-exponte exponenxonnt large space of DAGs, which may result in crude approximation of their posterior probability mass function. In our experiments, DiBS and its instantiation with SVGD are used interchangeably. and DiBS+ denotes the weighted particle mixture.
motivation is that most of the particles in the empirical distribution
6
Evaluation on Synthetic Data
6.1 Experimental Setup
Synthetic data We compare DiBS to set of related methods in marginal and joint posterior
inference of synthetic linear and nonlinear Gaussian BNs. Our setup follows [9, 24, 27, 52], who consider inferring BNs with Erdös- Rénós-Réénii.e and scale-free random structures [53, 54]. For each graph,

---

E-SHD
AUROC
neg. MLL
neg. I-MLL
60
Erdös R-nyinnyi
0.96
900
50-
1200
0.88
800
40-
0.80
1050
700
0.72
900
600
20
0.64
750
504
0.96
960
780
Scale-free
40-
0.88
880
30-
0.80
800
20-
0.72
0.64
720
10-
MCa
BPC
BGES
DiBs
DIBS4
Maa
BPC BGES
DiBs
DiBS
Mca
BPC
EGES
DIBS
DIBS
Mca BPC BGES
DiBS
DIBS4
Figure 2: Marginal posterior inference of 20-node linear Gaussian BNs using the BGe marginal likelihood. Higher scores on AUROC and lower scores on E-SHD, neg. MLL, neg. I-MLL are
preferred. DiBS+ performs competitively across all metrics, in particular E-SHD and neg. (I-)MLL
here with d E (20, 50) nodes and 2d edges expectation, we sample a set of ground truth parameters and then generate training, held-out, and interventional data sets. In all settings, use N 100 observations for inference, emulating the use case of Bayesian structure learning where the uncertainty about the graph structure is significant.
Graph priors For Erdös-Rényí Rényi graphs, all methods use the prior p(G) x qllG)i (1 -i)-G
capturing that each edge exists independently w.p. [53]. For scale-free graphs, we define the prior GGG) x I (1+ G; ||1)- analogous to their power law degree distribution p(deg) ~ deg [54]. Here, G; the i-th column of the adjacency matrix. DiBS implements either prior by using the corresponding term above as additional factor in p(Z) with G Ga (Z) (see Section 4.2).
Metrics Since neither density nor samples of the ground truth posteriors are available for BNs of d E (20,50) variables, we follow the evaluation metrics used by previous work. We define the expected structural Hamming distance to the ground truth graph G under the inferred posterior as
E-SHD(p. G') = Σ D(G SHD(G,
(18)
where SHD G*) counts the edge changes that separate the essential graphs representing the MECs of G and G [8, 55]. In addition, we follow Friedman and Koller [20] and Ellis and Wong [41] and compute the area under the receiver operating characteristic curve (AUROC) for pairwise edge predictions when varying the confidence threshold under the inferred marginal p(gr = 1 (D). Finally, following Murphy [10], we also evaluate the ability to predict future observations by computing the average negative (marginal) log likelihood on 100 held-out observations Diest,
neg. LL(p, Dlest)
Σ P(G. θ | log (DDies G,O
G, O)
(19)
When inferring P D). the corresponding neg. MLL metric instead uses p(Diett G). Analogously, we also compute the interventional log likelihoods I-LL and I-MLL, relevant metric in causal infer- ence [10, 12]. Here, an interventional data set (Dint I) is instead used to compute p(Dint G, O, and p( Dint G, in (19), respectively. Scores are the average of 10 interventional data sets. All reported metrics this section are aggregated for inference of 30 random synthetic BNs.
In the remainder, DiBS is always run for 3,000 iterations and with k = d for inference of d-variable BNs, which leaves the matrix of edge probabilities unconstrained in rank. We discard DiBS particle in the rare case that a returned graph is cyclic. Complete details on Gaussian BNs, the evaluation metries, hyperparameters, and all baselines can be found in Appendix E.
6.2 Linear Gaussian Bayesian Networks
Marginal posterior inference For linear Gaussian BNs, we first evaluate the classical setting of inferring the marginal posterior p(G (D since p(D|G G) can be computed in closed form. To this end, employ the commonly used Bayesian Gaussian Equivalent (BGe) marginal likelihood, which
7

---

E-SHD
AUROC
neg. L.
neg. I-LL
1.0
Brdös- Brdös-Rényj Rényi
75
2500
60
0.9
2000
2400
45
0.8
1500
1800
30-
0.7
1000
1200
15-
500
600
60
2000
1800
0.96
Scale-free
45
0.88
1600
1500
1200
1200
30-
0.80
900
15-
0.72
800
600
M-Mc.
G-MC.
BPC BGES
DiBs
DIDS4
M-MG. G-MO9 BFC EGES.
DiBs
DIBS
400
M-MC G-MCR BPC BGES' DIBS
DiBS
300
M-MC G-MOs DPC* BGES DiBs DIBS
Figure 3: Joint posterior inference of graphs and parameters of linear Gaussian networks with d = 20 nodes. DiBS+ performs best across all of the metrics. BGES*, the next-best alternative, yields substantially worse performance in E-SHD, i.e., in recovering the overall graph structure and MEC. Recall that higher AUROC and lower E-SHD, neg. LL, and neg. I-LL scores are preferable.
scores Markov equivalent structures equally [16, 17]. The form of the BGe score requires DiBS to use the score function estimator in (14).
We compare DiBS with the nonparametric DAG bootstrap [42] using the constraint-baser PC [7] and the score-based GES [6] algorithms (BPC and BGES). For MCMC, we only consider structure MCMC (MC*) [37] comparison. Order MCMC or hybrid DP approaches bound the number of parents and thus often exclude the ground truth graph a priori, especially for scale-free BN structures. Burn-in and thinning for MC are chosen to make the wall time comparable with DiBS run on CPUs. In the remainder of the paper, each method uses 30 samples to approximate the posterior over BNs.
Figure 2 summarizes the results for 30 randomly generated BNs with d 20 nodes. We find that DiBS+ performs well compared to the other methods, all of which were specifically developed for the marginal inference scenario evaluated here. DiBS+ appears to be preferable to DiBS.
Joint posterior inference When inferring the joint posterior p(G, θ I we can employ a more explicit representation of linear Gaussian BNs, where the conditional distribution parameters are standard Gaussian. Here, DiBS can leverage the Gumbel-softmax estimator in (12) because pi G, θ 'D) is well-defined when G lies on the interior of [0. 1]dxd (see Appendix E. To provide a comparison with DiBS in the absence of an applicable MCMC method, we propose two variants of MC3 (M-MC*) jointly samples parameters and structures,
MC as baselines. Metropolis-Hastinge and Metropolis- wietlii--iitti---ee Gibbs MC (G-MC3) alternates in proposing structure and parameters [56]. Moreover, extend the bootstrap methods by taking the closed-fomm maximum likelihood estimate [57] as the posterior parameter sample for a given graph inferred using the BGe score (BPC* and BGES*), approach taken in, e.g., causal BN learning [14].
Figure 3 shows the results for = 20 nodes, where E-SHD and AUROC are computed by empirically marginalizing out the parameters. We find that DiBS+ is the only considered method that performs well across all of the metries, often outperforming the baselines by a significant margin. As for marginal posterior inference of linear Gaussian BNs, DiBS+ performs slightly better than DiBS.
6.3 Nonlinear Gaussian Bayesian Networks
We also consider joint inference of nonlinear Gaussian BNs where the mean of each local conditional Gaussian parameterized by 2-layer dense neural network with five hidden nodes and ReLU activation functions (see Appendix 1. Since the marginal likelihood does not have a closed form, we are unable use BPC* and BGES* as means of comparison. Figure displays the results for d = 20 variables, where a given BN model has ( 2,220 weights and biases. Analogous to joint inference of linear Gaussian BNs, DiBS and DiBS+ outperform the MCMC baselines across the considered metrics. To the best of our knowledge, this is the first time that such nonlinear Gaussian BN models have been inferred under the Bayesian paradigm, which opens up exciting avenues the active learning of more sophisticated causal structures.
8

---

E-SHD
AUROC
neg. LL
neg. I-LL
Erdös-Rede Rényi
100
0.90
1.2x109
1.2x10$
0.75
8x10
8x104
0.60
4x10
4x104
40
20
0.45
0-
122110
1.2x10
Scale-free
0.90
9x104
9x104
0.75
6x10
6x104
30
0.60
3x10
3x10
20
0.45
M-MC
G-MC
DIBS
DiBS
M-MC
G-MOu
D
DiBS+
M-MC. G-MCa DiBS
DiBS
M-MC G-MO3 DiBS DiBS
Figure 4: Joint posterior inference of nonlinear Gaussian BNs with 20 nodes, Here, the mean of
the local conditional distribution of each node is parameterized by a 2-layer neural network with five hidden nodes. DiBS and DiBS+ perform favorably across the board, partieularly in the graph metrics.
Appendix F complements Sections 6.2 and 6.3 with details on computing time and efficient imple- mentation of DiBS with SVGD, showing that GPU wall times of the above inference experiments lie on the order of seconds or only few minutes. In addition, Appendix G provides results for = 50 variables, where DiBS+ likewise performs favorably when jointly inferring p(G, θ D). For marginal posterior inference of p(G. D) under the BGe marginal likelihood, DiBS appears to require more Monte Carlo samples compensate for the high variance of the score function estimator in this nigh-dimensional setting.
6.4 DiBS with SVGD: Additional Analyses and Ablation Studies
Having compared DiBS and its instantiation with SVGD to existing approaches, we finally devote Appendix H to empirically analyzing some properties of the algorithm. One of our key results is that, all other things held equal, substituting our inner product model in (6) with (gts Z) = aa Pa(2)) where single scalars encode the edge probabilities, results in significantly worse evaluation metrics.
We additionally study the uncertainty quantification in (non adentiieteee edge structures and show the effects of reducing the latent dimensionality i or the number of iterations T. Our findings suggest that reducing either hyperparameter still allows for competitive posterior approximations and enables trading off posterior inference quality with computational efficiency, e.g., in large-scale applications.
7 Application: Inferring Protein Signaling Networks From Cell Data
A widely used benchmark in structure learning is the proteomics data set by Sachs et al. [3]. The data contain = 7, 766 continuous measurements of d = 11 proteins involved in human immune system cells as well as an established causal network of their signaling interactions.
Table 1: Inference of protein signaling pathways with Gaussian BNs. Metrics are the mean + SD of 30 random restarts.
E-SHI
AUROC
MC
34.0 0.7
0.616 = 0.027
We infer both linear and nonlinear Gaussian BNs with Erdös-Rényi graph priors exactly as in Section 6. The AUROC results in Table 1 indicate that the posterior by DiBS under the BGe model provides the most calibrated edge confidence scores. Not penalizing model complexity as much, the marginal BGe posterior of DiBS (DiBS+) av-
BPC
25.5
0.566 + 0.020
BGES
33.7
0.641 0.030
DIBS
0.5
0.647 0.047
DiBS+
4..7
0.629 + 0.045
M-MC"
37.3
0.551 0.078
G-MCa
3.2
0.527 = 0.067
DIBS
23.4
0.598 + 0.052
erages a high expected number of 39.0 (35.0) edges, com- pared to 12.7 (14.2) and 12.6 (14.2) for its joint posteriors over linear and nonlinear BNs, respectively. Appendix provides further details and analyses on this matter. The E-SHD scores show that, among all the methods, DiBS is closest in structure to the consensus network when performing joint inference with nonlinear Gaussian BNs. This further highlights the need for nonlinear conditionals and joint inference of the graph and parameters in complex real-worke settings.
DiBS+
2.7
0.557 + 0.052
M-MC.
2.02
0.526 0.084
G-MC²
35.1 3.2
0.540 + 0.080
9 DiBS
22,6 - 0.5
0.577 0.039
DIBS+
22.8 1.9
0.535
0.041
Lineur Gaussiun BN; graph only
marginal lik.
5 Lincar Gaussian BN: graph and parameters jointly Nonlinear Gaussian BN; graph and parameters jointly

---

8 Conclusion
We have presented general, fully differentiable approach inference of posterior distributions BNs. Our framework is based on a continuous latent representation of DAGG, whose posterior can be equivalently -nftrred-without loss of generality and using existing black box inference methods. While we have used SVGD [19] for this purpose, our general approach could also be instantiated with, e.g., gradient-based sampling methods that rely on the score of the target density [58, 59]. This may improve upon the asymptotic runtime of DiBS with SVGD, which scales quadratically in the number of sampled particles. We expect that our end-to-end approach can be extended to handle missing and interventional data as well as to amortized contexts, where rich unstructured data is available.
Broader impact Our work is relevant to any scientific discipline that aims at inferring the (causal) structure of a system or reasoning about the effects of interventions. algorithms and decisions are grounded in the structural understanding of a data- gata-gennrating system and take into account the epistemic uncertainty, we expect them to be more robust and have fewer unforeseen side-effects. However, the assumptions allowing for causal interpretation of DAGs, e.g., the absence of unmeasured confounders, are often untestable and to be taken with care [60], particularly in aufett--ritical and societally-sensitivee sensitive applications. Hence, while potential misuse can never be ruled out, our presented method predominantly promises positive societal and scientific impact.
Acknowledgments and Disclosure of Funding
This project received funding from the Swiss National Science Foundation under NCCR Automation under grant agreement 51NF40 180545, the European Research Council (ERC) under the European Union's 's Horizon 2020 research and innovation program grant agreement no. 815943, and was supported with compute resources by Oracle Cloud Services. This work was also supported by the German Federal Ministry of Education and Research (BMBF): Tübingen AI Center, FKZ: 18031, and by the Machine Learning Cluster of Excellence, EXC number 2064/1 Project number 390727645. We thank Nicolo Ruggeri and Guillaume Wang for their valuable feedback.
References
[1] Judea Pearl. Probubilistic reasoning in intelligent systems: networks of plausible inference. 1988.
[2] Dana Pe' er, Aviv Regev, Gal Elidan, and Friedman. Inferring subnetworks from perturbed expression profiles. Bioinformatics, 17:S215-S224, 2001.
[3] Karen Sachs, Omar Perez, Dana Pe'er, Douglas Lauffenburger, and Garry P Nolan. Causal protein- signaling networks derived from multiparameter single-cell data. Science,
2005.
[4] Chikako Van Koten and AR Gray. An application bayesian network for predicting object-oriented software maintainability. Information and Software Technology, 48(1):59-67, 2006.
[5] Philippe Weber, Gabricla Mcdina-Olite Oliva, Christophe Simon, and Benoit Tung. Overview on bayesian networks applications for dependability, risk analysis and maintenance arcas. Engineering Applications of Artificial Intelligence, 25(4):671 1--622, 2012.
[6] David Maxwell Chickering. Optimal structure identification with greedy search. J. Mach. Learn. Res., 3:
507-554, March 2003.
[7] Peter Spirtes, Clark N. Glymour, and Richard Scheines. Causation, prediction, and search. Adaptive computation and machine learning. MIT Press, Cambridge, Mass, 2nd ed edition, 2000.
[8] Toannis Tsamardinos, I Laura E Brown, and Constantin Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
[9] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[10] Kevin P Murphy, Active learning of causal bayes net structure. 2001.
[11] Simon Tong and Daphne Koller. Active learning for structure in bayesian networks. In International joint conference artificial intelligence, volume 17, pages 863 869, 2001.
[12] Hyunghoon Cho, Bonnie Berger, and Jian Peng. Reconstructing causal biological networks through active learning. PloS one, 1061), 2016.
10

---

[13] Robert Osazuwa Ness, Karen Sachs, Parag Mallick, and Olga Vitek. A bayesian active learning experimen- tal design for inferring signaling networks. In International Conference on Rexearch in Computational Molecular Biology. pages 134-156. Springer, 2017.
[14] Raj Agrawal, Chandler Squires, Karren Yang. Karthikeyan Shanmugam, and Caroline Uhler. ABCD- Strategy: Budgeted experimental design for targeted causal structure discovery. In Proceedings of Machine Learning Research, volume 89, pages 3400-3409. PMI.R, 16-18 Apr 2019.
[15] Julius von Kügelgen, Paul K Rubenstein, Bernhard Schölkopf, and Adrian Weller. Optimal experimental design via bayesian optimization: active causal structure learning for gaussian process networks. arXiv preprint arXiv:1910.0396 2019.
[16] Dan Geiger and David Heckerman. Learning gaussian networks. In Proceedings of the Tenth International Conference on Uncertainty in Artificial Intelligence, UAI'94, page -35-243, San Francisco, CA. USA, 1994.
[17] Dan Geiger and David Heckerman. Parameter priors for directed acyelic graphical models and the charac- terization of several probability distributions. Ann. Statist., 30(5):1412-1444 October 2002. Publisher: The Institute of Mathematical Statistics.
[18] David Heckerman, Dan Geiger, and David M. Chickering. Learning Bayesian Networks: The Combination of Knowledge and Statistical Data. Machine Learning, 20(3):197-244 243, September 1995.
[19] Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS'16, 16, page 2378-2386, 2016.
[20] 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, January 2003.
[21] David Madigan and Adrian E Raftery. Model selection and accounting for model uncertainty in graphical models using occam's 's window. Journal of the American Statistical Association, 89(428) 1994.
[22] David Madigan, Jonathan Gavrin, and Adrian Raftery, Eliciting prior information to enhance the predictive performance of bayesian graphical models. Communications in Statisticx-Theur and Methods, 24(9):2271-2292. 1995.
[23] R.W. Robinson. Counting labeled acyclic digraphs. New Directions in the Theory of Graphs, 1973.
[24] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN; DAG structure learning with graph neural networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 7154-7163, PMLR, 09-15 Jun 2019.
[25] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural dag learning. In International Conference on Learning Representations, 2020.
[26] Ignavier Ng. Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. A graph autoencoder approach to causal structure learning. In NeurlPS 2019 Workshop "Do the right thing": machine learning and causal inference for improved decision muking, 2019.
(27] Xun Zheng, Chen Dan, Bryon Aragam, Pradcep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, pages 3414 3425. PMLR, 2020.
[28] Roxana Pamfil, Nisara Sriwattnnaworcchai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 1595-160S. PMLR, 2020.
[29] Mengyue Yang, Furui Liu, Zhitang Chen, Xinwei Shen, Jianye Hao, and Jun Wang. Causalvae: Disen- tangled representation learning via neural structural causal models. arXiv preprint arXiv: 2004.08697, 2020.
[30] Muhan Zhang, Shali Jiang, Zhicheng Cui, Roman Garnett, and Yixin Chen, D-vae: A variational autoencoder for directed acyclic graphs. In Advances in Neural Information Processing Systems, volume 32, 2019.
[31] Nan Rosemary Kc, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Chris Pal, and Yoshua Bengio. Learning neural causal models from unknown interventions. arXiv preprint arXiv: 1910.01075, 2019.
[32] Ignavier Ng. Zhuangyan Fang, Shengyu Zhu, Zhitang Chen, and Jun Wang. Masked gradient-based causal structure learning. arXiv preprint arXiv: 1910008577 2019.
[33] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. arXiv preprint arXiv:2007.0173 2020.
[34] Dennis Wei, Tian Gao, and Yue Yu. Dags with no fears: A closer look at continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 2020.
11

---

[35] Ignavier Ng. Sébastien Lachapelle, Nan Rosemary Ke, and Simon Lacoste-Julien. On the convergence of continuous constrained optimization for structure learning. arXiv preprint arXiv:2011.111 1.1101011111,,,11 2020.
[36] David Madigan, Jeremy York, and Denis Allard. Bayesian graphical models for discrete data. International Stutistical Review/Revue Internationale de Statistique, pages 215-232, 1995.
[37] Paolo Giudici and Robert Castelo. Improving markov chain monte carlo model search for data mining. Machine learning, 50(1-2):127-158, 2003.
[38] Daniel Eaton and Kevin Murphy, Bayesian strueture learning using dynamic programming and MCMC. In Proceedings of the Twenty-Thiie Conference on Uncertainty in Artificial Intelligence (UAT), 2007.
[39] Marco Grzegorzzyk and Dirk Husmeier. Improving the structure MCMC sampler for Bayesian networks by introducing new edge reversal move. Machine Learning, 71(2)265, April 2008
[40] Jack Kuipers and Giusi Moffa, Partition meme for inference on acyelic digraphs. Journal of the American Statistical Association, 1
2017.
[41] Byron Ellis and Wing Hung Wong. Learning causal bayesian network structures from experimental data. Journal of the American Statistical Association, 103(482):778-788 2008
[42] N. Friedman, M. Goldszmidt, and A. J. Wyner. Data analysis with Baycsian networks: A bootstrap approach. In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, 1999.
[43] Mikko Koivisto and Kismat Sood. Exact bayesian structure discovery in bayesian networks. Journal of Machine Learning Research, 5(May):549-55)) -73, 2004.
[44] Mikko Koivisto. Advances in exact bayesian structure discovery in bayesian networks. In Proceedings of the Twenty-Second Conference on Uncertainty in Artificaal Intelligence, page 241-248, 2006.
[45] Marina Meila and Tommi S Jaakkola. Tractable bayesian learning of tree belief networks. Statistics and Computing, (16):77-22, 2006.
[46] Denver Dash and Gregory F Cooper. Model averaging for prediction with discrete bayesian networks. Journal of Machine Learning Research, 1(Sep))1177-1100) 2004.
[47] Thomas N Kipf and Max Welling. Variational graph auto-encoders. arXiv preprint arXiv: 1,00000 2016.
[48] Chris J. Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: continuous relaxation of discrete random variables. In Proceedings International Conference on Learning Representations, 2017.
[49] Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparametrization with gumhel-softmax. In Proceed- ings International Conference on Learning Representations 2017, April 2017.
[50] Ronald J Williams. Simple statistical gradient-foloownng algorithms for connectionist reinforcement learning. Machine learning, 8(3-4):229-256 1992.
[51] Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih. Monte carlo gradient estimation in machine learning. Journal of Machine Learning Research, (112)):1162, 2020.
2020.
[52] 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. 2020.
[53] P. Erdōs and A Rényi On random graphs. Publicationes Mathematicae, 6:290-297, 1959.
[54] Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439): 509-512, 1999.
[55] Steen A Andersson, David Madigan, Michael D Perlman, et al. A characterization of markov equivalence classes for acyclic digraphs. Annals of statistics, 25(2):505-541, 1997.
[56] W Keith Hastings. Monte carlo sampling methods using markov chains and their applications. 1970.
[57] Alain Hauser and Peter Bühlmann. Jointly interventional and observational data: estimation of interven- tional markov equivalence classes directed acyelic graphs. Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 291-318 2015.
[58]
Radford M Neal et al. Meme using hamiltonian dynamics. Handbook of Markov chain Monte Carlo, 2 (11):,, 101
[59] Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML'11, page 68 681-688, 2011.
[60] A. Philip Dawid, Beware of the DAG! In Proceedings of Workshop on Causaliry: Objectives and Assessment at NIPS 2008, pages 59-86, 2010.
[61] Jonas Peters and Peter Bühlmann. Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1):219225 1282 2014.
12

---

[62] Jack Kuipers, Giusi Moffa. and David Heckerman. Addendum on the scoring of Gaussian directed acyclic graphical models. Ann. Statist., 1661,,1689-116,, August 2014. Publisher: The Institute of Mathematical Statistics.
[63] Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez- Pooee--Pz.,aa and Michèle Sebag. Structural agnostic modeling: Adversarial learning of causal graphs. arXiv preprint arXiv: 04999, 2018.
[64]
Chris Sherlock, Paul Fearnhead, Gareth O Roberts, et al. The random walk metropolis: linking theory and practice through a case study. Statistical Science, 25(2):172-190) 2010.
[65] Dorit Dor and Michael Tarsi, A simple algorithm to construct consistent extension of partially oriented graph. 1992.
[66] Yangbo He, Jinzhu Jia, and Bin Yu. Counting and exploring sizes of markov equivalence classes of directed acyelic graphs. The Journal of Machine Learning Research, 16(1):2589-260% 2015.
[67] Diviyan Kalainathan and Olivier Goudet. Causal discovery toolbox: Uncover causal relationships in python. arXiv preprint arXiv: / 1903, 02278, 2019. URL https://githut con./FenTeechhooutiions httpe//eooooooaaeooooeaeeaaaaeee
(cuusalDiscoverytlolooo
[68] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: compos- able transformations of Python+NumPyy programs, 2018. URL http: //tht::///gthbbb com/ gowgle/oge/aaa
[69] Yoshua Bengio. Tristan Deleu, Nasim Rahaman, Nan Rosemary Ke. Sebastien Lachapelle, Olexa Bila- niuk, Anirudh Goyal, and Christopher Pal. A meta-transfer objective for learning to disentangle causal mechanisms. In International Conference Learning Representations, 2020.
13


Reviews:
[{"review_id": "nLfG89kcPeC", "reviewer": "Reviewer_GxXy", "summary": "This work proposes a new approach to posterior approximation in the context of structure learning based on the Stein Variational Gradient Descent (SVGD) procedure and the differentiable acyclicity constraint inspired from Zheng et al 2018. A continuous latent variable Z (parent of G only) is introduced in the model to allow the computation of the gradient of its posterior (with respect to Z) which in turn enables one to use SVGD to sample from P(G|D). This gradient can be estimated using the Gumbel-softmax trick or score function estimator. The framework is validated on synthetic data as well as the protein signaling dataset from Sachs et al..", "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": "Originality:\n\nThe method proposed is original, however, I believe two papers strongly related to this work have been omitted, namely, “Masked Gradient-Based Causal Structure Learning” by Ng et al. and “Differentiable Causal Discovery from Interventional Data” by Brouillard et al.. Both works leverage the Gumbel-softmax trick to learn an adjacency matrix (like this work), although without framing their approach as Bayesian (thus without learning a posterior distribution over graphs). The latter work frames its continuous constrained optimization problem as a stochastic relaxation of the original discrete combinatorial problem. Given the strong resemblance (when using a single particle), I believe it is important that the paper contrasts these approaches with theirs. An empirical comparison (with the single particle version) would also be interesting.\n\nQuality:\n\nThe method proposed is well motivated and technically sound. I have not reviewed the proofs however. This is an interesting direction and I believe this a good paper. However, there are some things I would like the authors to clarify/adjust, mainly regarding experiments. I am open to revise my score once the following concerns are addressed:\n\nL219: Could you provide some motivation for this “weighted particle mixture” modification? Does SVGD allow for this kind of modification or is it a contribution from the authors? This is especially important since it appears to be performing better in some cases.\n\nI am under the impression that the metrics reported will favor methods which underestimate the uncertainty in the posteriors. For instance, if the samples from the approximate posterior (wrongfully) concentrate around the MAP estimator, E-SHD will be potentially lower than it would have been if the samples were coming from the actual posterior. I am not well-versed in Bayesian model selection so I don’t know what are the typical ways to evaluate the quality of the posterior, but I suspect it needs to somehow take into account proper uncertainty quantification. This aspect is studied in Appendix E.4 for a very special case of a graph with 4 nodes, allowing to compute the exact posterior and compare it with the one approximated by DiBS. I noticed that the number of SVGD iterations in this study is much smaller than in the experiment presented in the main paper (500 against 3000). I would like to see how the numbers from Table 4 change if one makes 3000 SVGD iterations. Does the uncertainty estimation remain the same? Given the importance that this work puts on being able to estimate p(G|D), I believe it would be important to make this kind of analysis on more graphs and more methods to see how DiBS compare. These experiments should also appear in the main text given their importance.\n\nL268: Is 30 samples from the posterior enough to estimate the various performance metrics? It would be interesting to see how much the metrics vary when one repeats the sampling. Do the performance metrics change a lot when using more samples from the posterior?\n\nClarity:\n\nThe paper is well written and I enjoyed reading it. Here are a few things that might still be clarified:\n\nIn Section 4.3, titled “Black box bayesian inference”, I was not sure what the actual algorithm suggested was. I believe the title is a bit misleading since no specific inference procedure is suggested in this section. The first paragraph of the following section also implicitly mentions that some inference procedure was suggested above: “... translates learning discrete graph structures into an inference problem over the continuous variable Z.” As an outsider to the Bayesian approach, the only inference procedure I thought of after reading only 4.3 is *maximum a posteriori estimation* via gradient ascent. Later reading the paragraph at line 212 confirmed my guess. Please clarify.\n\nEquation (15) felt a bit hard to interpret given that no explanation of the stein variational gradient descent was provided. Maybe adding some general explanation of this method in the background section would help, or maybe just refer to Algorithm 1 at the very beginning of Section 5.\n\nMinor:\n\nEquation (8): Why using E_{G | Z}[h(G)] instead of h(G_\\alpha(Z)) ? After all, h was designed to take as input weighted adjacency matrices. \n\nEquation (17) if <u_i, v_j> = 0, the limit is ½ right? Probably not relevant since this concerns a set of measure zero. \n\nFor each method, which thresholds were explored to compute the AUROC, and how were they selected? Also regarding the AUROC, how is p(g_ij = 1 | D) estimated? I assume by counting the proportion of graphs with g_ij = 1 among the 30 samples from the posterior? \n\nConfused by ijk in equation (8), I thought Z was a k x 2d matrix (line 129) \n\nHow is E-SHD estimated? I am guessing via Monte Carlo estimation using the samples from the posterior. To clarify.\n\nSignificance:\n\nI believe proper quantification of uncertainty in structure learning is important and this work makes an interesting contribution in that direction by using recent developments in continuous optimization for structure learning. This work has potential for impact.\n\n\n\n", "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"}, "extra_text": {}}, {"review_id": "k-Bz6xwVOu", "reviewer": "Reviewer_ZLxN", "summary": "This paper proposes a method called DiBS (differentiable Bayesian structure learning) for learning the structure of a Bayesian network (BN) from data based on a fully differentiable reparameterization of the relevant likelihood expressions, allowing for the use of continuous optimization-based Bayesian inference algorithms. This method allows for posterior inference of graph structure alone, or joint inference over structure and conditional distribution parameters. First, they augment the standard generative model of data from a Bayesian network with a continuous-valued latent variable of node embeddings such that the probability of each edge in the graph adjacency matrix can be modeled by the inner product similarity of node embeddings with a temperature parameter. The authors use a previously-introduced function of the adjacency matrix that calculates graph cyclicity to include a prior that encourages acyclic structure. They derive expressions that reformulate the necessary score functions in terms of gradients with respect to the new continuous latent graph representation, using the Gumbel-softmax trick to reparameterize the gradient where applicable and the score function gradient estimator otherwise. The method uses Stein variational gradient descent as the inference algorithm, optimizing a set of particles for the latent graph representation while annealing the temperature parameters to result in a final set of particles that correspond to adjacency matrices of valid directed acyclic graphs. A variant of the base method called DiBS+ that weights each particle by its unnormalized posterior probability is also explored. On synthetic data from linear Gaussian BNs, DiBS and DiBS+ show comparable performance to greedy search-based (BPC, BGES) and sampling-based ($\\text{MC}^3$) algorithms across a variety of metrics when inferring structure alone, and comparable or better performance when jointly inferring structure and parameters. On synthetic data with nonlinear local conditionals parameterized by a neural network, DiBS/DiBS+ perform comparably to or better than sampling-based methods for joint inference across metrics. These results are also reflected on real data from a protein signaling network, where DiBS/DiBS+ achieve the best recovery of the true network and the most well-calibrated edge confidence scores in the joint inference setting.", "questions": "", "limitations": "", "rating": 8, "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": "This paper presents a very interesting approach for joint Bayesian inference of structure and parameters of graphical models using continuous optimization-based methods, with the overall presentation of the method being clear with sufficient detail. While the performance of DiBS is comparable to other methods for marginal posterior inference in the linear Gaussian setting, it still provides a competitive alternative to existing approaches. DiBS seems to particularly excel when applied to joint inference of structure and parameters, especially with nonlinear conditionals, which is relevant for a lot of real-world applications and presents this method as a potentially valuable tool for practitioners interested in this kind of inference for their data. Given that DiBS is also agnostic to the particular inference algorithm used (as the authors mention in their Conclusion), any advancements in Bayesian inference algorithms that require a differentiable likelihood could be applied to improve upon this framework. Overall, I believe this is a strong paper and recommend it for acceptance.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors briefly discuss the potential negative impacts, specifically the risk in taking a causal interpretation of the results of this method applied to real data without considering the necessary assumptions.", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns come to mind.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "L6bgJc0eJ6W", "reviewer": "Reviewer_voY1", "summary": "This paper proposes a fully differentiable (i.e., continuous-space) method for Bayesian learning of Bayesian network structures as well as parameters. The proposed method (DiBS) is compared against several existing state-of-the-art methods on simulated and real data sets, and the results demonstrate very competitive performance.", "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": "This paper proposes a fully differentiable (i.e., continuous-space) method for Bayesian learning of Bayesian network structures as well as parameters. Authors propose a seemingly simple generative model (Fig. 1): a generative graph model is motivated by variational graph autoencoders and utilizes recent results on continuous characterization of acyclic graphs that are incorporated as priors in the generative model. Model inference is implemented with the Stein variational gradient descent, with the necessary technical terms being derived in the manuscript/appendix. The proposed method (DiBS) is compared against several existing state-of-the-art methods on simulated and real data sets, and the results demonstrate very competitive performance.\n\nThe manuscript is very clear and easy to read. To the best of my knowledge, this manuscript proposes a novel technique for structure/parameter learning for BNs that goes beyond the previous methods; it could even be fair to say that perhaps this manuscript opens a new direction for Bayesian learning for BN structure/parameters and for causal inference in general. Having worked on the BN learning in the past, I must say that I read the manuscript with great interest.  I have no real concerns about this work: all technical details seem correct, and simulations and results support the claims made in the manuscript.  I list some comments/suggestions and minor comments below that the authors can implement when revising their manuscript. \n\n1. Overall the literature review appears to be very comprehensive. While it is true that most of the MCMC techniques for BN structure inference assume a closed form for the marginal likelihood (which makes it easy to implement e.g. Metropolis-Hasting algorithm over structures), previous methods also include reversible jump MCMC methods for inferring both the structure and parameters (proposed at least in the context of dynamics BNs). \n\n2. Although Fig. 2 does not contain the most impressive results, I found it very convincing that the proposed method performs comparably with (or in terms of neg. (I-)MLL measures even slightly better than) the MC^3 for linear Gaussian BNs. \n\n3. Page 7, row 268: why did you use so few samples (S=30)?\n\n4. Appendix D.5, results for d=50: can you discuss about possible reasons why the method performs worse for the case where closed form for the MLL is available? \n\n5. Sec. 7, phosphoproteomics data analysis: I would recommend authors to consider using neg. log.likel. measures on hold-out data in addition to comparing against the \"ground truth network structure\", as the ground truth may not be entirely correct (e.g. missing links). \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "6d23Jsl2lf", "reviewer": "Reviewer_aXYX", "summary": "This paper addresses the structure learning problem, in which the goal is to find a directed graphical model (structure and parameters) that accurately models an observed dataset. The authors work in the Bayesian structure learning paradigm, defining a prior distribution over graphical models and aiming to infer a posterior. Rather than defining a discrete prior over the unknown graph structure $G$, the authors define a prior over a continuous latent variable $Z$, and a conditional distribution on $G$ given $Z$, such that the data likelihood, $P(D \\mid Z) = E_{G \\mid Z}(P(D \\mid G))$, is differentiable with respect to $Z$. This enables gradient-based inference methods to be applied to the problem of inferring $Z$, marginalizing $G$—if the necessary gradients can be efficiently estimated. In principle, any gradient-based inference method could be applied, but the authors spell out details of a Stein Variational Gradient Descent procedure. This proposed technique, and a variant in which an additional reweighting step is applied, are evaluated on synthetic data, against MCMC-based and non-parametric bootstrapping baselines.", "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": "This paper presents a novel approach to the important but difficult problem of Bayesian structure learning. Bayesian approaches to the problem are conceptually appealing (define a prior over graphs, and the obvious data likelihood), but hard to implement in practice, because they involve combinatorial search over a high-dimensional discrete space. This paper applies recent advances in differentiable characterizations of graph \"acyclicness\" to formulate a structure learning model to which (some) gradient-based inference methods can be applied.\n\nAlthough various components of this approach have been previously proposed (e.g., structure learning as continuous optimization, and \"inner product\" priors for edge probabilities in graphs), this work presents a valuable synthesis of these ideas into a coherent model and inference algorithm for Bayesian structure learning. The empirical validation is thorough, and I appreciated the variety of different metrics used to evaluate performance.\n\nOne aspect I found unsatisfying was the move from DiBS to DiBS+. The unweighted particle collection is already meant to represent the posterior; wouldn't weighting it by the unnormalized target \"double-count,\" further concentrating the inferred posterior? That this performs better than DiBS using the metrics investigated here suggests either that something is off about the metrics (i.e., they do not reward proper uncertainty?), or that the inference algorithm is under-confident / producing higher-uncertainty estimates than it should. Or am I misunderstanding something? One question for the authors is whether you have tried applying a similar re-weighting to samples from the MC^3 baselines you consider (weighting according to the posterior densities they target).\n\nOther minor questions:\n\n* In the BGe marginal likelihood, is the likelihood 0 for graphs G with cycles?\n* Maybe I missed it, but could you report -log p(D^{test} | G*, Theta*) for the true data-generating graph? This may help calibrate the negative log likelihood numbers you report.\n* In Section 7, where does \"the consensus network\" come from?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Broader impact is explicitly discussed (with a focus on the potential positives), but I don't see an explicit discussion of limitations. Some questions that might be worth considering:\n\n* Is it really the case that other black-box inference methods would be easy to apply to your examples (L332)? For example, I don't believe Hamiltonian Monte Carlo can easily make use of the gradient estimators you give in Section 4.3.\n\n* The score function estimator (Eq 14) is often described as having unacceptably high variance. You appear to have used it successfully in Section 6.2, but do you expect it will scale to e.g. models with parameters $\\Theta$ that can't be analytically marginalized, or to higher-dimensional data? ", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4.5 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers agree that the proposed use of continuous embeddings for Bayesian learning of directed graphical models is innovative and potentially high-impact.  Reviews generally have very positive comments about the manuscript's clear presentation and promising experimental results.  In future revisions, please do better discuss the motivations for the "DiBS+" variant, and clarify the points about technical details and evaluation metrics that were discussed with reviewers.
Decision: Accept (Spotlight)