Efficient Neural Causal Discovery without Acyclicity Constraints
Authors: Phillip Lippe
Taco Cohen
Efstratios Gavves
Venue: ICLR 2022 Poster (2022)
Keywords: Causal discovery
structure learning

Abstract: Learning the structure of a causal graphical model using both observational and interventional data is a fundamental problem in many scientific fields. A promising direction is continuous optimization for score-based methods, which, however, require constrained optimization to enforce acyclicity or lack convergence guarantees. In this paper, we present ENCO, an efficient structure learning method for directed, acyclic causal graphs leveraging observational and interventional data. ENCO formulates the graph search as an optimization of independent edge likelihoods, with the edge orientation being modeled as a separate parameter. Consequently, we provide for ENCO convergence guarantees under mild conditions, without having to constrain the score function with respect to acyclicity. In experiments, we show that ENCO can efficiently recover graphs with hundreds of nodes, an order of magnitude larger than what was previously possible, while handling deterministic variables and discovering latent confounders.

Published as a conference paper at ICLR 2022
EFFICIENT NEURAL CAUSAL DISCOVERY WITHOUT ACYCLICITY CONSTRAINTS
Phillip Lippe University of Amsterdam
Taco Cohen
Efstratios Gavves
Qualcomm AI Research" tacos@qti qua comm. com
University of Amsterdam QUVA lab
QUVA lab p. p.lippe@uva.mn nl
e w.gaveeeuva...
ABSTRACT
Learning the structure of a causal graphical model using both observational and interventional data is a fundamental problem in many scientific fields. A promising direction is continuous optimization for score-based methods, which, however, re- quire constrained optimization to enforce acyclicity or lack convergence guarantees. In this paper, we present ENCO, an efficient structure learning method for directed, acyclic causal graphs leveraging observational and interventional data. ENCO formulates the graph search as optimization of independent edge likelihoods, with the edge orientation being modeled as a separate parameter. Consequently, we provide for ENCO convergence guarantees when interventions all variables are available, without having to constrain the score function with respect to acyclicity. In experiments, show that ENCO can efficiently recover graphs with hundreds of nodes, an order of magnitude larger than what was previously possible, while handling deterministic variables and discovering latent confounders.
INTRODUCTION
Uncovering and understanding causal mechanisms is important problem not only in machine learning (Schölkopf et al., 2021; Pearl, 2009) but also in various scientific disciplines such as computational biology (Friedman et al., 2000; Sachs et al., 2005), epidemiology (Robins et al., 2000; Vandenbroucke et al., 2016), and economics (Pearl, 2009; Hicks et al., 1980). A common task of interest causal structure learning (Pearl, 2009; Peters et al., 2017), which aims at learning directed acyelic graph (DAG) in which edges represent causal relations between variables. While observational data alone is in general not sufficient identify the DAG (Yang et al., 2018; 8; Hauser & Bühlmann, 2012), interventional data can improve identifiability up to finding the exact graph (Eberhardt et al., 2005; Eberhardt, 2008). Unfortunately, the solution space of DAGs grows super-exponentially with the variable count, requiring efficient methods for large graphs. Current methods are typically applied few dozens of variables and cannot scale so well, which is imperative for modern applications like learning causal relations with gene editing interventions (Dixit et al., 2016; Macosko et al., 2015).
promising new direction for scaling up DAG discovery methods are continuous-opiimization methods (Zheng et al., 2018; 2020; Zhu et al., 2020; Ke et al., 2019; Brouillard et al., 2020; Yu et al., 2019). In contrast to score-based and constrained-baser (Peters et al., 2017; Guo al., 2020) methods,
methods reinterpret the search over discrete graph topologies as a continuous problem with neural networks as function approximators, for which efficient solvers amenable. To restrict the search space to acyclic graphs, Zheng et al. (2018) first proposed to view the search as a constrained optimization problem using an augmented Lagrangian procedure to solve it. While several improvements have been explored (Zheng et al., 2020; Brouillard et al., 2020; Yu et al., 2019; Lachapelle et al., 2020), constrained optimization methods remain slow and hard to train. Alternatively, it is possible apply a regularizer (Zhu et al., 2020; Ke al., 2019) to penalize cyclic graphs. While simpler to optimize, methods relying on acyclicity regularizers commonly lack guarantees for finding the correct causal graph, often converging to suboptimal solutions. Despite the advances, beyond linear, continuous settings (Ng et al., 2020; Varando, 2020) continuous optimization methods still cannot scale to more than 100 variables, often due difficulties in enforcing acyclicity.
"Qualcomm Al Research initiative of Qualcomm Technologies, Inc.

---

Published as a conference paper at ICLR 2022
In this work, we address both problems. By modeling the orientation of an edge separate parameter, we can define the score function without any acyclicity constraints or regularizers. This allows for unbiased low-variance gradient estimators that scale learning to much larger graphs. Yet, we are able to intervene on all variables, the proposed optimization is guaranteed to converge to the correct, acyelic graph. Importantly, since such interventions might not always be available, we show that our algorithm performs robustly even when intervening on fewer variables and having small sample sizes. We call our method ENCO for Efficient Neural Causal Discovery.
We make the following four contributions. Firstly, we propose ENCO, a causal structure learning method for observational and interventional data using continuous optimization. Different from recent methods. ENCO models the edge orientation as a separate parameter. Secondly, we derive unbiased, low-variance gradient estimators, which is crucial for scaling up the model to large numbers of variables. Thirdly, we show that ENCO is guaranteed to converge to the correct causal graph if interventions on all variables are available, despite not having any acyclicity constraints. Yet, we show in practice that the algorithm works on partial intervention sets well. Fourthly, extend ENCO to detecting latent confounders. In various experimental settings, ENCO recovers graphs accurately, making less than one on graphs with 1,000 variables in less than nine hours of computation.
2 BACKGROUND AND RELATED WORK
2.1 CAUSAL GRAPHICAL MODELS
A causal graphical model (CGM) is defined by a distribution P over a set of random variables X [X,,, [Xi,.., NN} and a directed acyclic graph (DAG) G (V. Each node i E V corresponds to the random variable X; and an edge (i, j) E E represents a direct causal relation from variable X, to Xy: X; + Xy. A joint distribution is faithful to a graph G all and only the conditional independencies present in IP are entailed by G (Pearl, 1988). Vice versa, distribution PP is Markov to graph G if the joint distribution can be factorized as p( II P; (x.,)pa(XX)) where pa(X;) is the set of parents of the node i in G. An important concept which distinguishes CGMs from standard Bayesian Networks are interventions. An intervention on variable X, describes the local replacement of its conditional distribution (X.: (//////aXXiXa (pa(X.)) by a new distribution X(XX,)pp pp(XX))),, )) X, is thereby referred to as the intervention target. An intervention is "perfect" when the new distribution is independent of the original parents, iee. p( X; |aa Xx)) = p(X,).
2.2 CAUSAL STRUCTURE LEARNING
Discovering the graph G from samples of a joint distribution IP is called causal structure learning or causal discovery, fundamental problem in causality (Pearl, 2009; Peters et al., 2017). While often performed from observational data, i.e. samples from (see Glymour et al. (2019) for an overview), we focus in this paper on algorithms that recover graphs from joint observational and interventional data. Commonly, such methods are grouped into constraint-ban bassdd and score- based approaches.
Constraint-based methods use conditional independence tests to identify causal relations (Monti et al., 2019; Spirtes et al., 2000; Kocaoglu et al., 2019; Jaber et al., 2020; Sun et al., 2007: Hyttinen al., 2014). For instance, the Invariant Causal Prediction (ICP) algorithm (Peters et al., 2016; Christina et al., 2018) exploits that causal mechanisms remain unchanged under an intervention except the one intervened on (Pearl, 2009; Schölkopf al., 2012). We rely on a similar idea by testing for mechanisms that generalize from the observational to the interventional setting. Another line of work extend methods working on observations only interventions by incorporating those as additional constraints in the structure learning process (Mooij et al., 2020; Jaber et al., 2020).
Score-based methods, on the other hand, search through the space of all possible causal structures with the goal of optimizing a specified metric (Tssmardiioos et al., 2006; Ke et al., 2019; Goudet et al., 2017; Zhu et al., 2020). This metric, also referred to as score function, is usually a combination of how well the structure fits the data, for instance in terms of logllikeli as well as regularizers for encouraging sparsity. Since the search space of DAGs is super-exponential in the number of nodes, many methods rely greedy search, yet returning graphs in the true equivalence class (Meek, 1997; Hauser & Bühlmann, 2012; Wang et al., 2017: Yang et al., 2018). For instance, GIES (Hauser & Bühlmann, 2012) repeatedly adds, removes, and flips the directions of edges in proposal graph
2

---

Published as a conference paper at ICLR 2022
Distribution fitting
Graph fitting
X1
NN fa.
-)
Altemate between both steps
X
X2
Ysz
X2
NN fe.
p(Ksl.))
OTS
PIS
X3
Sax
P(G)
Figure 1: Visualization of the two training stages of ENCO, distribution fitting and graph fitting, on an example graph with 3 variables X2, The graph on the right further shows how the parameters y and 0 correspond to edge probabilities. We learn those parameters by comparing multiple graph samples on how well they generalize from observational to interventional data.
until no higher-scoring graph can be found. The Interventional Greedy SP (IGSP) algorithm (Wang et al., 2017) is a hybrid method using conditional independence tests in its score function.
methods score-based methods that avoid the combinatorial greedy search over DAGs by using gradient-basee based methods (Zheng et al., 2018; 18; Ke et al., 2019; Lachapelle al., 2020; Yu al., 2019; Zheng et al., 2020; Zhu al., 2020; Brouillard et al., 2020). Thereby, the adjacency matrix is parameterized by weights that represent linear factors or probabilities of having an edge between a pair of nodes. The main challenge of such methods is how to limit the search space to acyelic graphs. One common approach is to view the search as constrained optimization problem and deploy an augmented Lagrangian procedure to solve it (Zheng et al., 2018; 2020; Yu et al., 2019; Brouillard et al., 2020), including NOTEARS (Zheng et al., 2018) and DCDI (Brouillard et al., 2020). Alternatively, Ke et al. (2019); Ng et al. (2020) propose to use regularization term penalizing cyclic graphs while allowing unconstrained optimization. However, the regularizer must be designed and weighted such that the correct, acyelic causal graph is the global optimum of the score function.
Continuous-optimizz
EFFICIENT NEURAL CAUSAL DISCOVERY
3.1 SCOPE AND ASSUMPTIONS
We consider the task of finding a directed acyelic graph G= (V, F)) with N variables of an unknown CGM given observational and interventional samples. Firstly, we assume that: (1) The CGM is causally sufficient, i.e., all common causes of variables are included and observable; (2) We have N interventional datasets, each sparsely intervening on different variable; (3) The interventions are "perfect" and "stochastic", meaning the intervention does not set the variable necessarily to a single value. Thereby, do not strictly require faithfulness, thus also recovering some graphs violating faithfulness. We emphasize that we place no constraints on the domains of the variables (they can be discrete, continuous, or mixed) or the distributions of the interventions. We discuss later how to extend the algorithm to infer causal mechanisms in graphs with latent confounding causal variables. Further, we discuss how extend the algorithm to support interventions to subsets of variables only.
3.2 OVERVIEW
ENCO learns a causal graph from observational and interventional data by modelling a probability for every possible directed edge between pairs of variables. The goal is that the probabilities corresponding to the edges of the ground truth graph converge to one, while the probabilities of all other edges converge to zero. For this to happen, we exploit the idea of independent causal mechanisms (Pearl, 2009; Peters et al., 2016), according to which the conditional distributions for all variables in the ground-truth CGM stay invariant under an intervention, except for the intervened ones. By contrast, for graphs modelling the same joint distribution but with flipped or additional edge, this does not hold (Peters et al., 2016). In short, we search for the graph which generalizes best from observational to interventional data. To implement the optimization, learning stages, that is distribution fitting and graph fitting, visually summarized in Figure 1. Ke et al. (2019) proposed similar two- who-phase framework, but importantly differ in the graph parameterization.
alternate between two
Distribution fitting trains a neural network fo. per variable X; parameterized by de to model its observational, conditional data distribution p(X.l...). The input to the network are all other variables,
3

---

Published as a conference paper at ICLR 2022
X-1- For simplicity, we want this neural network to model the conditional of the variable X, with respect to any possible set of parent variables. We, therefore, apply a dropout-like scheme to the input to simulate different sets of parents, similar as (Ke et al., 2019; Ivanov et al., 2019; Li et al., 2020; Brouillard et al., 2020). In that case, during training, we randomly set an input variable X, to zero based the probability of its corresponding edge Xj X;, and minimize
min ExEM |- log fo, (X; M O X-)),,
(I)
where M₁ ~ Ber(p(X; + X,)). For categorical random variables X;, we apply a softmax output activation for s and for continuous ones, we use Normalizing Flows (Rezende & Mohamed, 2015).
Graph fitting uses the learned networks to score and compare different graphs on interventional data. For parameterizing the edge probabilities, we use two sets of parameters: E RNN represents the existence of edges in graph, and R-N the orientation of the edges. The likelihood of an edge is determined by p(X;
X,) a (7i5) (00r), with (....) being the sigmoid function and 8i = -0 The probability of the two orientations always sum to one. The benefit of separating the edge probabilities into two independent parameters y and @ is that it gives more control over the gradient updates. The existence of an (undirected) edge can usually be already learned from observational or arbitrary interventional data alone, excluding deterministic variables (Pearl, 2009). In contrast, the orientation can only be reliably detected from data for which an intervention is performed on its adjacent nodes, i.e., X; or Xj for learning Ozg. While other interventions eventually provide information on the edge direction, e.g., intervening on a node X₁ which is a child of X; and a parent of Xj. we do not know the relation of X2 to X, and X, at this stage, as we are in the process of learning the structure. Despite having just one variable for the orientation, Yas and Tiss are learned as two separate parameters. One reason is that on interventional data, an edge can improve the log-liken li-iikeeihood estimate in one direction, but not necessarily the other, leading to conflicting gradients.
We optimize the graph parameters 7 and by minimizing
N + Asparse ΣΣ
غ = (Eipp(())Bpp E E (C)
Σ Lo(X))
(m)) - o(1)))
where pT( is the distribution over which variable to intervene on (usually uniform), and 5y (X) the joint distribution of all variables under the intervention î. In other words, these two distributions represent our interventional data distribution. With py.e(C'). we denote the distribution over adjacency matrices C under 7 0, where Ci estimate of variable X; conditioned on the parents according to C: (X,))) = log fo, (X;; C.,, O X=:). The second term of Equation 2 is an fr-reullarizer the edge probabilities. It acts prior, selecting the sparsest graph of those with similar likelihood estimates by removing redundant edges.
((r(()(()) Lc (X,) is the negative lieelioel
Prediction. Alternating between the distribution and graph fitting stages allows to fine-tune the neural networks to the most probable parent sets along the training. After training, we obtain a graph prediction by selecting the edges for which (T:;) and (8er) are greater than 0.5. The orientation parameters prevent loops between any two variables, since (0:s) can only be greater than 0.5 in one direction. Although the orientation parameters do not guarantee the absence of loops with more variable, we show that under certain conditions ENCO yet converges to the correct, acyclic graph.
3.3 LOW-VARIANCE GRADIENT ESTIMATORS FOR EDGE PARAMETERS
To update y and 0 based on Equation 2, we need to determine their gradients through the expectation (C)) where C is discrete variable. For this, we apply REINFORCE (Williams, 1992). For clarity of exposition, we limit the discussion here to the final results and provide the detailed derivations in Appendix A. For parameter Kij. we obtain the following gradient:
θ Ong A C= ((r)) (00))) Ex,c_,, Exx.xx((,))
(X,Ax,((,)) + Aegarse]
(3)
where Ex,c is summarizes for brevity the three expectations in Equation 2, excluding the edge X, + X,, from C. Further, Lx,- (Xş) denotes the negative logllikelihood for X,, if we include the edge X, + X, to the adjacency matrix C-aj. i.e., Cj = 1, and Lx, (X,) if we set C.y = O. The gradient in Equation 3 can be intuitively explained: if by the addition of the edge X; Xj, the
4

---

Published as a conference paper at ICLR 2022
log-likelihood estimate of Xj is improved by more than Aspanee increase the corresponding edge parameter "lys otherwise, we decrease it.
We derive the gradients for the orientation parameters 0 similarly. As mentioned before, we only take the gradients for 0,g when we perform an intervention on either X; or Xj. This leads us to:
① Ê = o' d'(0.) (B(x))(())))) ( o(1))) BIx,, (x, Lx,Ax, (X,)) ,8,,
(4)
p(xx,) (15i)) Bix, X,C- [xx (X,) -0000-xxxpp Lx, (X.]]
The probability of taking an intervention on X; is represented by pl p(xx;) (usually uniform across variables), and Fxxx,xXC. x.C yf the same expectation as before under the intervention on X,, When the oriented edge X, + X3 improves the log- likelihood of X, under intervention on X,, then the first part of the gradient increases 0zy. In contrast, when the true edge is X, + X,, the correlation between X; and X, learned from observational data would yield worse likelihood estimate of X; on interventional data on X, than without the edge X; + X; This is because pl p(X,) (x,,...) does not stay invariant under intervening on X;, The same dynamic holds for interventions on Xj. Lastly, for independent nodes, the expectation of the gradient is zero.
Based Equations 3 and 4, we obtain a tractable, unbiased gradient estimator by using Monte-Carlo sampling. Luckily, samples can be shared across variables, making training efficient. We first sample an intervention, a corresponding data batch, and K graphs from Py, (C)) (K usually between 20 and 100). We then evaluate the log likelihoods of all variables for these graphs on the batch, and estimate fx,+x, (X,) and A (X,) for all pairs of variables X; and X; by simply averaging the results for the two cases separately. Finally, the estimates are used to determine the gradients for y and A.
Low variance. Previous methods (Ke al., 2019; Bengio et al., 2020) relied on a different REINFORCE-lik estimator proposed by Bengio et al. (2020). Adjusting their estimator to our setting of the parameter Yiy; for instance, the gradient looks as follows:
EE [(((T))) -Cy)) -Xo(),
gi = o Ex
+ >Apparse
(5)
Ec [Cxxx]]
where Pú represents the gradient of Yt. Performing Monte- Carlo sampling for estimating the gradient leads to biased estimate which becomes asymptotically unbiased with in- creasing number of samples (Bengio et al., 2020). The di-
10"
ENCO (Ours)
wug.graa gradientann std of
W-REINFORCE
vision by the expectation of (X,) 3)x is done for variance reduction (Mahmood et al., 2014). Equation 5, however, is still sensitive to the proportion of sampled Ciyy being one or zero. A major benefit of gradient formulation in Equation 3, instead, is that it removes this noise by con- sidering the difference of the two independent Monte-Carlo estimates - (X,) and Ax,Ax, (X,). Hence, we can use smaller sample size than previous methods and attain 10 times lower standard deviation, as visualized in Figure 2.
10
Avg.
10²2
0
1000
2000
3000
4000
Number of graph samples
Figure ENCO estimates gradients of significantly lower variance com- pared to (Bengio et al., 2020).
3.4 CONVERGENCE GUARANTEES
Next, we discuss the conditions under which ENCO convergences to the correct causal graph. We show that not only does the global optimum of Equation 2 correspond to the true graph, but also that there exist no other local minima ENCO can converge to. We outline the derivation and proof of these conditions in Appendix B. and limit our discussion here to the main assumptions and implications.
To construct theoretical argument, we make the following assumptions. First, we assume that sparse interventions have been performed on all variables. Later, we show how to extend the algorithm to avoid this strong assumption. Further, given a CGM, we assume that its joint distribution (X) is Markovian with respect to the true graph G. In other words, the parent set pa(X,) reflects the inputs to the causal generation mechanism of X,. We assume that there exist no latent confounders in g. Also, we assume the neural networks in ENCO are sufficiently large and sufficient observational data is provided to model the conditional distributions of the CGM up to an arbitrary small error.
Under these assumptions, we produce the following conditions for convergence:
5

---

Published as a conference paper at ICLR 2022
Theorem 3.1. Given a causal graph 9 with variables X and conditional observational distributions p(X.l...) the proposed method ENCO will converge to the true, causal graph G. if the following conditions hold for all edges X, + -x in G:
1. For all possible sets of parents of X, excluding X, by adding X, the log-likelihood estimate of x is improved or unchanged under the intervention on X;:
Vpa(X;) S X-isj : Exx,
x (appp((X((((())))) X.) log (o, 2 0
(6)
2. For at least one set of nodes ). for which the probability to be sampled as parents of X is greater than O, Equation 6 must be strictly greater than zero.
3. The effect of Xx on X2 cannot be described by other variables up to Xsparse:
min (iC
Bn Bi, p(X; pa, X.) - log p( [pa]
>
Aaporse
where gpa, is the set of nodes excluding X; which, according to the ground truth graph. could have an edge to X3 without introducing cycle, and p1_, refers to the distribution over interventions pi PT(T) excluding the intervention on variable X--
Further, for all other pairs X,, X3 for which Xj is descendant of X,, condition I and 2 must hold.
Condition and 2 ensure that the orientations can be learned from interventions. Intuitively, ancestors and descendants in the graph have to be dependent when intervening the ancestors. This aligns with the technical interpretation Theorem 3.1 that the likelihood estimate of the child variable must improve when intervening and conditioning that the sparsity regularizer needs to be selected such that it chooses the sparsest graph among those graphs with equal joint distributions as the ground truth graph, without trading sparsity for worse distribution estimates. The specific condition in Theorem 3.1 ensures thereby that the set can be learned with gradient-based algorithm. We emphasize that this condition only gives an upper bound for Asparse when sufficiently large datasets are available, In practice, the graph can thus be recovered with sufficiently small sparsity regularizer and dependencies among variables under interventions. We provide more details for various settings and further intuition in Appendix B.
its ancestor variables. Condition 3 states intuitively
Interventions on fewer variables. It is straightforward extend ENCO to support interventions on fewer variables. Normally, in the graph fitting stage, we sample one intervention at a time. We can, thus, simply restrict the sampling only to the interventions that are possible (or provided in the dataset). In this case, we update the orientation parameters O,y of only those edges that connect to an intervened variable, either X; or Xj. as before. For all other orientation parameters, we extend the gradient estimator include interventions on all variables. Although this estimate is more noisy and does not have convergence guarantees, it can still be informative about the edge orientations.
Enforcing acyclicity When the conditions are violated, e.g. by limited data, cycles can occur in the prediction. Since ENCO learns the orientations as a separate parameter, we can remove cycles by finding the global order of variables O E SN, with N being the set of permutations, that maximizes the pairwise orientation probabilities: arg maxo II,, I T (fo..,,). ). This utilizes the learned ancestor-descendant relations, making the algorithm more robust to noise in single interventions.
3.5 HANDLING LATENT CONFOUNDERS
So far, we have assumed that all variables of the graph are observable and can be intervened on. A common issue in causal discovery the existence of latent confounders, i.e., an unobserved common cause of two or more variables introducing dependencies between each other. In the presence of latent confounders, structure learning methods may predict false positive edges. Interestingly, the context of ENCO latent confounders for two variables X,,X, cause unique pattern of learned parameters. When intervening on X; or Xj. having an edge between the two variables is disadvantagcous, as in the intervened graph X; and X, are (conditionally) independent. For interventions on all other variables, however, an edge can be beneficial as X, and X, are correlated.
Exploiting this, we extend ENCO to detect latent confounders. We focus on latent confounders between two variables that do not have any direct edges with each other, and assume that the confounder is not child of any other observed variable. For all other edges besides between X;

---

Published as a conference paper at ICLR 2022
Table 1: Comparing structure learning methods in terms of structural hamming distance (SHD) on common graph structures (lower is better), averaged over 25 graphs each. ENCO outperforms all baselines, and by enforcing acyclicity after training, can recover most graphs with minimal errors.
Graph type
bidiag
chain
collider
full
jungle
random
GIES
33.6
17.5
(42.9
216.5 (+15.2)
33.1 (#2.9)
(+1522)
IGSP
32.7 (=5.1)
253.8 23.7(12.3) (12.3))
35.9 (45.2)
(18.0))))
GES Orientating
14.8 (+2.6)
(=0.7
(+2.4)
282.8 (+4.2)
14.7
60.1 (+8.9)
SDI
(#2.6)
(=2.0)
16.1 (#2.4)
153.9 (=10.3)
(.2.3))
10.8 (#3.9)
DCDI
16.9 (12.0)
10.1 (11.1)
10.9 (13.6)
21.0 (14.8)
17.9
7.7 (13.2)
ENCO (ours)
2.2 (41.4)
(=1-)
(..
9.2 (43.4)
(+7.)))
4.6 (+1.9) 0.2 (1:0.5
ENCO-acyclic (ours)
0.0 (1:0.0)
0.0 (1:0.))
(.1.6))))
5.3 (12.3)
(.1.))
and X,, we can still rely on the guarantees in Section 3.4 since Equation 7 already includes the score every pair of variables on how likely they share a latent confounder using a function lc(.) that maximized in the scenario (D) (0) where 7aj is only updated with gradients
possibility of additional edges in such cases. After convergence,
mentioned above. For this, we define Tig = 7ij + 7ij
(1)
(O) all others. With this separation, define
from Equation 3 under interventions on X₁, and Tig
the following score function which is maximized by latent confounders:
Ic(X,, Xy) = (7.0
(1-0 (1-)))((--0 (1-O
(8)
To converge to the mentioned values, especially of 7/00, we need similar condition as in Equation 7: the improvement on the log-likelihood estimate gained by the edge X; + X, and conditioned on all other parents of Xj needs to be larger than Aquase on interventional data excluding X; and Xj- If this is not the case, the sparsity regularizer will instead remove the edge between X; and X, preventing any false predictions among observed variables. For all other pairs of variables, least one of the terms in Equation 8 converges to zero. Thus, we can detect latent confounders by checking whether the score function (e, X,) is greater than a threshold hyperparameter T E (0.0, 1.0). We discuss possible guarantees in Appendix B, and experimentally verify this approach in Section 4.5.
4 EXPERIMENTS
We evaluate ENCO on structure learning on synthetic datasets for systematic comparisons and real- world datasets for benchmarking against other methods in the literature. The experiments focus on graphs with categorical variables, and experiments on continuous data are included in Appendix D.5. Our code is publicly available at https: / /g thub. com/ //m/pphll ippe ippe/ENCO.
4.1 EXPERIMENTAL SETUP
Graphs and datasets. Given ground-truth truth causal graphical model, all methods are tasked to recover the original DAG from a set of observational and interventional data points for each variable. In case of synthetic graphs, we follow the setup of Ke al. (2019) and create the conditional distributions from neural networks. These networks take as input the categorical values of its variable's ' parents, and are initialized orthogonally to output a non-trivial distribution.
Baselines. We compare ENCO to GIES (Hauser & Bühlmann, 2012) and IGSP (Wang et al., 2017: Yang et al., 2018) as greedy score-based approaches, and DCDI (Brouillard et al., 2020) and SDI (Ke et al., 2019) as continuous optimization methods. Further, as common observational baseline, we apply GES (Chickering, 2002) on the observational data to obtain graph skeleton, and orient each edge by learning the skeleton on the corresponding interventional distribution. We perform separate hyperparameter search for all baselines, and use the same neural network setup for SDI, DCDI, and ENCO. Appendix C provides a detailed overview of the hyperparameters for all experiments.
4.2 CAUSAL STRUCTURE LEARNING ON COMMON GRAPH STRUCTURES
We first experiment on synthetic graphs. We pick six common graph structures and sample 5,000 observational data points and 200 per intervention. The graphs chaán and full represent the
7

---

Published as a conference paper at ICLR 2022
bidiag
collider
jungle
20
40
Method
20
SHD
SHD
DCDI
10
20
ENCO (Ours)
ENCO ENCO-acyclie (Ours)
0
4 810
15
20
25
4.6 810
15
20
25
4 6 810 Num. intervened variables
15
20 25
Num. intervened variables
Num. intervened variables
Figure 4: Experiments on graphs with interventions on fewer variables. Additional graphs are shown in Appendix D.2. ENCO outperforms DCDI on bidiag and jungle, even for very few interventions.
minimally- and maximally-connecee DAGs. The graph bi bidiag is a chain with 2-hop connections, and jungle a like-like graph. the collider graph, one node has all other nodes as parents. Finally, r random has a randomly sampled graph structure with a likelihood of 0.3 of two nodes being connected by a direct edge. For each graph structure, we generate 25 graphs with 25 nodes each, on which we report the average performance and standard deviation. Following common practice, we use structural hamming distance (SHD) as evaluation metric. SHD) counts the number of edges that need to be removed, added, or flipped in order to obtain the ground truth graph.
Table shows that the continuous optimization methods outperform the greedy search approaches on categorical variables. SDI works reasonably well on sparse graphs, but struggles with nodes that have many parents. DCDI can recover the collider and full graph to a better degree. yet degrades for sparse graphs. ENCO performs well all graph structures, outperforming all baselines. For sparse graphs, cycles can occur due to limited sample size. However, with enforcing acyclicity, ENCO-acyclic is able recover four out of six graphs with less than one experiments with various sample sizes in Appendix D.1. While other methods do not reliably recover the causal graph even for large sample sizes, ENCO attains low errors even with smaller sample sizes.
average. We further include
4.3 SCALABILITY
Next, we test ENCO on graphs with large sets of variables. We create zandom graphs ranging from N = 100 to N = 1,000 nodes with larger sample sizes. Every node has on average edges and maximum of 10 parents. The challenge of large graphs is that the number of possible edges grows quadratically and the number of DAGs super-exponentially. requiring efficient methods.
Structural Hamming Distance
ENCO (Ours)
SDI
DCDI
5000
4000
3000
2000
1000
We compare ENCO to the two best performing baselines from Table 1. SDI and DCDI. All meth- ods were given the same setup of neural networks and a maximum runtime which corresponds to 30 epochs for ENCO. We plot the SHD over graph size and runtime in Figure 3. ENCO recovers the causal graphs perfectly with no eтrors except for the 1.000-node graph, for which it misses only one out of million edges in 2 of 10 experiments. SDI and DCDI achieve considerably worse performance. This shows that ENCO can efficiently be applied to 1,000 variables while maintaining its convergence guarantees, underlining the benefit of its low-variance gradient estimators.
4030
0089
100 nodes
200 nodes (45mins)
400 nodes (2.5hra)
1000 rodes (Shrs)
(15 mins)
Graph size (Max. runtime)
Figure 3: Evaluating SDI, DCDI, and ENCO on large graphs in terms of SHD (lower is better). Dots represent single experiments, lines connect the averages. DCDI ran OOM for 1000 nodes.
4.4 INTERVENTIONS ON FEWER VARIABLES
We perform experiments the same datasets as in Section 4.2, but provide interventional data only for a randomly sampled subset of the 25 variables of each graph. We compare ENCO to DCDI, which supports partial intervention sets, and plot the SHD over the number of intervened variables in Figure 4. Despite ENCO's 's guarantees only holding for full interventions, it is still competitive and outperforms DCDI in most settings. Importantly, enforcing acyclicity has an even greater impact on fewer interventions as more orientations are trained on non-adjacent interventions (see Appendix B, for detailed discussion). We conclude that ENCO works competitively with partial interventions too.
8

---

Published as a conference paper at ICLR 2022
Table 3: Results on graphs from the BnLearn library measured in structural hamming distance (lower is better). Results are averaged over 5 seeds with standard deviations listed in Appendix C.5. Despite deterministic variables and rare events, ENCO can recover all graphs with almost no errors.
Dataset
cancer (5 nodes)
asia (8 nodes)
sachs (11 nodes)
child (20 nodes)
alarm (37 nodes)
diabetes (413 nodes)
pigs (441 nodes)
SDI
3.0
4.0
7.0
11.2
24.4
422.4
18.0
DCDI
4.0
5.0
5.4
8.4
30.0
ENCO (Ours)
0.0
0.0
0.0
0.0
1.0
2.0
0.0
4.5 DETECTING LATENT CONFOUNDERS
To test the detection of latent confounders, we create set of 25 random graphs with 5 additional latent con- founders. The dataset generated in the same way as before, except that we remove the latent variable from the input data and increase the observational and inter- ventional sample size (see Appendix C.3 for ablation studies). After training, we predict the existence latent confounder on any pair of variables X; and Xj 1,, is greater than т. We choose = = 0.4 but
Table 2: Results of ENCO on detecting la- tent confounders. The missed confounders do not affect other edge predictions.
Metries SHD
ENCO
0.0 (+-0.0)
Confounder recall
96.8% (#.9.5%)
Confounder precision
100.0% (+0.0%)
verify in Appendix C.3 that the method is not sensitive to the specific value of r. As shown in Table 2,
ENÇO detects more than 95% of the latent confounders without any false positives. What is more, the few mistakes do not affect the detection of all other edges, which recovered perfectly.
4.6 REAL-WORLD INSPIRED DATA
Finally, we evaluate ENCO on causal graphs from the Bayesian Network Repository (BnLearn) (Scutari, 2010). The repository contains graphs inspired by real-world applications that are used as benchmarks literature. In comparison to the synthetic graphs, the real-world graphs are sparser with maximum of parents per node and contain nodes with strongly peaked marginal distributions. They also include deterministic variables, making the task challenging even for small graphs.
We evaluate ENCO, SDI, and DCDI 7 graphs with increasing sizes, see Table 3. We observe that ENCO recovers almost all real-world causal graphs without errors, independent of their size. In contrast, SDI suffers from more mistakes as the graphs become larger. An exception is pigs (Scutari, 2010), which has maximum of 2 parents per node, and hence is easier to learn. The most challenging graph di abetes (Andreassen al., 1991) due to its large size and many deterministic variables. ENCO makes only two mistakes, showing that it can handle deterministic variables well. We discuss results on small sample sizes in Appendix C.5, observing similar trends. We conclude that ENCO can reliably perform structure learning on wide variety of settings, including deterministic variables.
CONCLUSION
We propose ENCO, an efficient causal structure learning method leveraging observational and interventional data. Compared to previous work, ENCO models the edge orientations separate parameters and uses an objective unconstrained with respect to acyclicity. This allows for easier optimization and low-variance gradient estimators while having convergence guarantees. As consequence, the algorithm can efficiently scale to graphs that are at least one order of magnitude larger graphs than what was possible. Experiments corroborate the capabilities of ENCO compared to the state-of-the-an on extensive array of settings graph sizes, sizes of observational and interventional data, latent confounding, as well as on both partial and full intervention sets.
Limitations. The convergence guarantees of ENCO require interventions all variables, although experiments on fewer interventions have shown promising results. Future work includes investigating guarantee extensions of ENCO to this setting. A second limitation is that the orientations missing transitivity: if X, > X2 and X₂ > X3. then X, > X₃ must hold. A potential direction is incorporating transitive relations for improving convergence speed and results fewer interventions.

---

Published as a conference paper at ICLR 2022
ETHICS STATEMENT
Causal structure learning algorithms such as the proposed method are mainly used to uncover and understand causal mechanisms from data. The knowledge of the underlying causal mechanisms can then be applied to decide on specific actions that influence variables or factors in desired way. For instance, by knowing that the environmental pollution in a city has impact the risk of cancer of its residents, one can try to reduce the pollution to decrease the risk of cancer. The applications of causal structure learning biology (Friedman et al., 2000; Sachs et al., 2005; Opgen-Rhein & Strimmer, 2007), epidemiology (Robins al., 2000; Vandenbroucke et al., 2016), and economics (Pearl, 2009; Hicks et al., 1980). We envision that our work can have positive impacts on those fields. One example we want to highlight the field of genomics. Recent advances have enabled perform gene knockdown experiments in large scale, providing large amounts of interventional data (Dixit et al., 2016; Macosko et al., 2015). Gaining insights into how specific genes and diseases interact can lead to the development of novel pharmaceutic methods for treating current diseases. Since the number of variables in those experiments is tremendous, efficient causal structure learning algorithms are needed. The proposed method constitutes first step towards this goal, and work foster future work for creating algorithms scaling beyond 10,000 variables.
ranging across many scientific disciplines, including computational
Since the possible applications are fairly vide-ranging there might be potential impacts we cannot forecast at the current time. This includes misuses of the method for unethical purposes. For instance, the method can be used to justify gender and race causes for irrelevant variables if the output is misinterpreted,, initial assumptions of the model are ignored, or the input data has been manipulated. Hence, the obligation use this method in correct way within ethical boundaries lies on the user. We emphasize this responsibility of the user in the license of our code.
REPRODUCIBILITY STATEMENT
To ensure reproducibility, have published the source code of the proposed method ENCO https: / github...... com/ ppl ippe Appe/ENCC. ENCO. The code includes instructions on how to download the datasets, and reproduce the experiments in Section 4 and additional experiments in Appendix D. Further, for all experiments of Section 4, we have included detailed overview in Appendix C of (a) the used data and its generation process, (b) all hyperparameters used for all methods, and (c) additional details on the results. All experiments have been repeated with 5 to 25 seeds to obtain stable, reproducible results. Appendix C.1.2 outlines the packages that have been used for running the baselines.
The computation resources deployed for all experiments are a 24-core CPU with a single NVIDIA RTX3090 GPU. All experiments can be reproduced computer with single GPU, and only the experiments on graphs larger than 100 variables require a GPU memory of about 24GB. The other experiments can be performed smaller GPUs well.
ACKNOWLEDGEMENTS
This work is financially supported by Qualcomm Technologies Inc., the University of Amsterdam and the allowance Top consortia for Knowledge and Innovation (TKIs) from the Netherlands Ministry of Economic Affairs and Climate Policy. We thank Pascal Mettes, Christina Winkler, and Sara Magliacane for their valuable comments and feedback on an initial draft of this work. We thank the anonymous reviewers for the reviews, suggestions, and engaging discussions which helped to improve this work further. Finally, we thank SURFsara for the support in using the Lisa Compute Cluster.
REFERENCES
Steen Andreassen, Roman Hovorka, Jonathan Benn, Kristian G. Olesen, and Ewart R. Carson. A model-based approach to insulin adjustment. In M. Stefanelli, A. Hasman, M. Fieschi, and J. Talmon (eds.), Lecture Notes in Medical Informaties, volume 44 of Lecture Notes in Medical Informatics, pp. 239,249, Germany, 1991. Springer. Conference date: 24-06-1991 Through 27-06-1991.
10

---

Published as a conference paper at ICLR 2022
Ingo A. Beinlich, Henri Jacques Suermondt, R. Martin Chavez, and Gregory F. Cooper. The ALARM Monitoring System: A Case Study with two Probabilistic Inference Techniques for Belief Networks. In Jim Hunter, John Cookson, and Jeremy Wyatt (eds.), AIME 89, pp. 247-256 Berlin, Heidelberg, 1989. Springer Berlin Heidelberg. ISBN 978-3-642-933377-7
Yoshua Bengio, Tristan Deleu, Nasim Rahaman, Nan Rosemary Ke, Sébastien Lachapelle, Olexa Bilaniuk, Anirudh Goyal, and Christopher J. Pal. A Meta-Transfer Objective for Learning to Disentangle Causal Mechanisms. In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020.
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexan- dre Drouin. Differentiable Causal Discovery from Interventional Data. In Hugo Larochelle, Marc' Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan- Lin (eds.), Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Process- ing Systems 2020, NeurIPS 2020, December 6-12. 2020, virtual, 2020.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507 2002.
Heinze-Deml Christina, Meinshausen Nicolai, and Peters Jonas. Invariant Causal Prediction for Nonlinear Models. Journal of Causal Inference, 6(2):1335, 2018.
Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. John Wiley and Sons, Ltd, 2005. ISBN 9780471748823.
A. Dixit, O. Parnas, B. Chen, C. P. Fulco, J.rerby-Aronn, N. D. Marjanovic, D. Dionne, T. Burks, R. Raychowdhury, B. Adamson, T. M. Norman, E. S. Lander, J. S. Weissman, N. Friedman, and A. Regev. Perturb-Seq: Dissecting Molecular Circuits with Scalable Single-Cell RNA Profiling of Pooled Genetic Screens. Cell, 1853311853-1866 2016.
Frederick Eberhardt. Almost Optimal Intervention Sets for Causal Discovery. In Proceedings of the Twenty-Foure Conference on Uncertainty in Artificial Intelligence, UAI'08, pp. 161-168, Arlington, Virginia, USA, 2008. AUAI Press. ISBN 0974903949.
Frederick Eberhardt, Clark Glymour, and Richard Scheines. On the Number of Experiments Suffi- cient and in the Worst Case Necessary to Identify All Causal Relations among N Variables. In Proceedings of the Twenty-Firss Conference Uncertainty in Artificial Intelligence, UAI'05, pp. 178-184, Arlington, Virginia, USA, 2005. AUAI Press. ISBN 0974903914.
Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe'er. er. Using Bayesian networks to analyze expression data. Journal of computational biology, (33-4):60-662) 2000.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of Causal Discovery Methods Based on Graphical Models. Frontiers in Genetics, 10:524, 2019. ISSN 1664-8021.
Olivier Goudet, Diviyan Kalainathan, Philippe Caillou, Isabelle Guyon, David Lopez-Paz, and Michèle Sebag. Causal generative neural networks. arXiv preprint arXiv:1711.0m 1.77911..866... 2017.
Ruocheng Guo. Lu Cheng, Jundong Li, Richard Hahn, and Huan Liu. Survey of Learning Causality with Data: Problems and Methods. ACM Comput. Surv., 53(4), 2020. ISSN 0360-0300.
Alain Hauser and Peter Bühlmann. Characterization and Greedy Learning of Interventional Markov Equivalence Classes of Directed Acyelic Graphs. Journal of Machine Learning Research, 13(1): 2409-2464, 2012. ISSN 1532-4435.
John Hicks et al. Causality in economics. Australian National University Press, 1980.
Kurt Horiik, Maxwell Stinchcombe, and Halbert White. Multilayer feedforward networks are universal approximators. Neural Networks, 2(5):359-385 1989. ISSN 0893-6080.
Chin-Wa Huang, David Krueger, Alexandre Lacoste, and Aaron C. Courville. Neural Autoregressive Flows. In Jennifer G. Dy and Andreas Krause (eds.). Proceedings of the 35th International Conference on Machine Learning ICML 2018, Stockholmsmässann Stockholm, Sweden, July 10-15, 2018, volume 80 of Proceedings of Machine Learning Research, pp. 2083-2092. PMLR, 2018
11

---

Published as a conference paper at ICLR 2022
Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo. Constraint- bonstraint-based Causal Discovery: Conflict Resolution with Answer Set Programming. In Nevin L. Zhang and Jin Tian (eds.), Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI 2014, Quebec City, Quebec. Canada, July 23-27, 2014, pp. 340-349. AUAI Press, 2014.
Oleg Ivanov, Michael Figurnov, and Dmitry P. Vetrov. Variational Autoencoder with Arbitrary Conditioning. In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019, 2019.
Amin Jaber, Murat Kocaoglu, Karthikeyan Shanmugam, and Elias Bareinboim. Causal Discovery from Soft Interventions with Unknown Targets: Characterization and Learning. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), Advances in Neural Information Process- ing Systems, volume 33, pp. 9551-9561.
Curran Associates, Inc., 2020.
Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. Causal Discovery Toolbox: Uncovering causal relationships in Python. J. Mach. Learn. Res., 21:37-1, 2020.
Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Chris Pal, and Yoshua Bengio. Learning neural causal models from unknown interventions. arXiv preprint arXiv:1910.01000 2019.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In Yoshua Bengio and Yann LeCun (eds.), 3rd International Conference on Learning Representations, ICLR 2015, San Diego CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
Murat Kocaoglu, Amin Jaber, Karthikeyan Shanmugam, and Elias Barcinboim. Characterization and Learning of Causal Graphs with Latent Variables from Soft Interventions. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Buc, E. Fox, and R. Gamett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Kevin B. Korb and Ann E. Nicholson. Bayesian artificial intelligence, second edition. CRC Press, Australia, 2010. ISBN 9781439815915.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien Gradient-Based Neural DAG Learning. In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020.
S. Lauritzen and D. Spiegelhalter. Local computations with probabilities graphical structures and their application to expert systems. Journal of the royal statistical society series b-methodological, 50:415 5484 1988.
Yang Li, Shoaib Akbar, and Junier Oliva. ACFlow: Flow Models for Arbitrary Conditional Likeli- hoods. In Proceedings of the 37th International Conference on Machine Learning, ICMI, 2020, 13-18 July 2020, Virtual Event, volume 119 of Proceedings of Machine Learning Research, pp. 5831-5841. PMLR, 2020.
R Duncan Luce. Individual choice behavior. John Wiley, Oxford, England, 1959.
E. Z. Macosko, A. Basu, R. Satija, J. Nemesh, K. Shekhar, M. Goldman, I. Tirosh, A. R. Bialas, N. Kamitaki, F. M. Martersteck, J. Trombetta, D. A. Weitz, J. R. Sanes, A. K. Shalek, A. Regev, and S. A. McCarroll. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell, 161(5): 2015.
Ashique Rupam Mahmood, Hado van Hasselt, and Richard S. Sutton. Weighted importance sampling for off- policy learning with linear function approximation. In Zoubin Ghahramani, Max Welling, Corinna Cortes, Neil D. Lawrence, and Kilian Q Weinberger (eds.), Advances in Neural Informa- tion Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-1 2014, Montreal, Quebec, Canada, pp. 3014-3022, 2014.
Christopher Meek. Graphical Models: Selecting causal and statistical models. PhD thesis, Carnegie Mellon University, 1997.
12

---

Published as a conference paper at ICLR 2022
Ricardo Pio Monti, Kun Zhang, and Aapo Hyvärinen. Causal Discovery with General Jon-Linear Relationships using Non-Lineat ICA. In Amir Globerson and Ricardo Silva (eds.), Proceedings of the Thirty-Fifth Conference Uncertainty in Artificial Intelligence, UAI 2019, Tel Aviv, Israel, July 22-25, 2019, volume 115 of Proceedings of Machine Learning Research, pp. 186-195. AUAT Press, 2019.
Joris M. Mooij, Sara Magliacane, and Tom Claassen. Joint Causal Inference from Multiple Contexts. Journal of Machine Learning Research, 21(99):1-108, 2020.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the Role of Sparsity and DAG Constraints for Learning Linear DAGs. In Advances in Neural Information Processing Systems, 2020.
Rainer Opgen-Rheee Rhein and Korbinian Strimmer. From correlation to causation networks: a simple approximate learning algorithm and its application to aigh-dimensional plant gene expression data. BMC Systems Biology, 1(1):77, 2007. ISSN 1752-0509.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang. Zachary Vetoto Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' d'Alché-Bus Emily B. Fox. and Roman Garnett (eds.). Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pp. 8024-800 0044--03,5250 2019.
Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. ISBN 0934613737
Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009. ISBN 052189560X.
Jonas Peters and Peters Bühlmann. Structural Intervention Distance (SID) for Evaluating Causal Graphs. Neural Computation, 27(3):771-799) 2015.
Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Statistical Sociery. Series (Statistical Methodology), pp. 947-1012, 2016.
Jonas Peters, Dominik Janzing, and Bernhard Schlkopf. Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press, 2017. ISBN 0262037319.
R. L. Plackett. The Analysis of Permutations. Journal of the Royal Statistical Society, Series (Applied Statistics), 24(2):193-202, 1975. ISSN 00359254, 14679876.
Danilo Jimenez Rezende and Shakir Mohamed. Variational Inference with Normalizing Flows. In Francis R. Bach and David M. Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, ICML. 2015, Lille, France, 6-11 July 2015, volume 37 of JMLR Workshop and Conference Proceedings, pp. 1530-1538. JMLR.org, 2015.
J.M. Robins, M.A. Hernan, and B Brumback. Marginal Structural Models and Causal Inference in Epidemiology. Epidemiology, 11(5):500-560, 2000.
K. Sachs, O. Perez, D. Pe'er, D. Lauffenburger, and G. Nolan. Causal Protein- Sieeiin-iinnaing Networks Derived from Multiparameter Single-Cell Data. Science, 308:523-529, 2005.
Bernhard Schölkopf, Dominik Janzing, Jonas Peters, Eleni Sgouritsa, Kun Zhang, and Joris M. Mooij. On causal and anticausal learning. In Proceedings of the 29th International Conference on Machine Learning, ICML 2012, Edinburgh, Scotland, UK, June 26 July 1. 2012. ieml.cc Omnipress, 2012.
Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEFF, 109(5):612-6) 664,, 2021.
13

---

Published as a conference paper at ICLR 2022
Marco Scutari. Learning Bayesian Networks with the bnleamn R Package. Journal of Statistical Software, 1-2,, 2010.
D. J. Spiegellhller and R. G. Cowell. Learning in probabilistic expert systems. Bayesian Statistics, 4: 447-466, 1992.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search, Second Edition. Adaptive computation and machine learning. MIT Press, 2000. ISBN 978-0-262-194400
Xiaohai Sun, Dominik Janzing, Bernhard Schölkopf. and Kenji Fukumizu. A bernel-based causal learning algorithm. In Zoubin Ghahramani (cd.), Machine Learning, Proceedings of the Twenty- Fourth International Conference (ICML 2007), Corvallis, Oregon, USA, June 20-24, 2007, volume 227 of ACM International Conference Proceeding Series, pp. 855 862. ACM, 2007.
Ioannis 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. ISSN
0885-6125.
Jan P Vandenbroucke, Alex Broadbent, and Neil Pearce. Causality and causal inference in epidemiol- ogy: the need for a pluralistic approach. International journal of epidemiology, 1766 2016.
Gherardo Varando. Learning DAGs without imposing acyclicity. arXiv preprint 2:00::2006.03005, 2020.
Yuhao Wang, Liam Solus, Karren D. Yang, and Caroline Uhler. Permutation-basee Causal Inference Algorithms with Interventions. In Isabelle Guyon, Ulrike von Luxburg. Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA, pp. 5822-5831, 2017.
Ronald J. Williams. Simple statistical soudient-following algorithms for connectionist reinforcement learning. In Machine Learning. pp. 256, 1992.
Karren D. Yang, Abigail Katoff, and Caroline Uhler. Characterizing and Learning Equivalence Classes of Causal DAGs under Interventions. In Jennifer G. Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning. ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of Proceedings of Machine Learning Research, pp. 5537-5546. PMLR, 2018.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG Structure Learning with Graph Neural Networks. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th Inter- national Conference on Machine Learning, ICML 2019, 9-15 June 2019. Long Beach, California, USA, volume 97 of Proceedings of Machine Learning Research, pp. 7154-7163. PMLR, 2019.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurlPS 2018, December 3-8, 2018, Montréal, Canada, pp. 1492-9503, 2018
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. Learning Sparse Nonparametric DAGs, In Silvia Chiappa and Roberto Calandra (eds.), The 23rd International Conference on Artificial Intelligence and Statistics, AISTATS 2020, 26-28 August 2020, Online [Palermo, Sicily, Itly], volume 108 of Proceedings of Machine Learning Research, pp. 3425-3425. PMLR, 2020.
Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal Discovery with Reinforcement Learning, In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020.
14

---

Published as a conference paper at ICLR 2022
SUPPLEMENTARY MATERIAL EFFICIENT NEURAL CAUSAL DISCOVERY WITHOUT ACYCLICITY CONSTRAINTS
TABLE OF CONTENTS
Gradient estimators
15
A.1 Low-variance gradient estimator for edge parameters
16
A.2 Low-variance gradient estimator for orientation parameters
19
A.3 Training loop
19
Conditions for converging to the true causal graph
B.I Assumptions
21
B.2 Convergence conditions
B.3 Convergence conditions for latent confounder detection
31
B.4 Convergence for partial intervention sets
................................
32
B.5 Example for non-faithful graphs
33
C Experimental details
35
C.I Common graph structure experiments
35
C.2 Scalability experiments
38
C.3 Latent confounder experiments
40
C.4 Interventions on fewer variables
C.5 Real-world inspired experiments
42
D Additional experiments
44
D.1 Effect of the sample size
44
D,2 Interventions on fewer variables
46
D.3 Ablation study on gradient estimators
46
D.4 Deterministic variables
47
D.5 Continuous data
49
D.6 Skeleton learning with observational baseline
49
D.7 Non-neural based data simulators
50
A GRADIENT ESTIMATORS
The following section describes in detail the derivation of the gradient estimators discussed in Section 3.3. We consider the problem of causal structure learning where parameterize the graph by edge existence parameters 7 and orientation parameters 9. Our objective is to optimize ry and such that we maximize the probability of interventional data, i.e. data generated from the true
15

---

Published as a conference paper at ICLR 2022
graphs under (arbitrary) interventions single variables. Thereby, the likelihood estimates have been trained on observational data only. Additionally, we want to ensure that the graph is as sparse possible to prevent unnecessary connections. Thus, an & regularizer is added top of the edge probabilities. The full objective can be written as follows:
N Σ (x(),)
N o (nyy (0,)))))))))
È = E (x)
Asparsc
where:
is the number variables in the causal graph (XI..., Xx):
Pr (I) is the distribution over interventions that are performed. This distribution can set as hyperparameter to weight certain interventions higher than others. In our experiments, we assume be uniform across interventions on variables;
(X) is the joint distribution of all variables under the intervention î;
(C)) is the distribution over adjacency matrices C, which we model as a product of independent edge probabilities: (C) I Iī=1,=; (71s)) ((0y)):
Lx(X,) is the negative log- ligelhod estimate of variable X; under sampled adjacency matrix C: (c(X,) log fo. (X,; ; C. @ X_:);
Asparse a hyperparameter representing the regularization weight.
Based on this objective, derive the gradient estimators for optimizing both edge existence and orientation parameters.
A.1 LOW-VARIANCE GRADIENT ESTIMATOR FOR EDGE PARAMETERS
In order optimize the edge parameters via SGD, need determine the gradient a É. Since L consists of sum of two terms, i.e., the log- likelihood estimate and the regularization, we can look at both parts separately. To prevent any confusion of index variables, we will use k, as indices for the parameter "Iu for which we determine the gradient, i.e., Trae
2%
a_ L. and i, as indices for sums.
As first step, we determine the gradients for the regularization term. Those can be found by taking derivative of the sigmoid:
N
((H))o((1))))) a(v)) (nna)) o A
(10)
Apparse
1=1=
Thus, is straight-forward to calculate for any edge parameter. In the following, we use ('...) ...) to abbreviate the derivate of the sigmoid: o' d'(7)) (1 (7w)))
For the log-likelihood term, we start by reorganizing the expectations to simplify the gradient 8 expression. The derivate term Jāu can be moved inside the two expectations over interventional data since those independent of the graph parameters. Thus, we can write:
2 Oryka
L'
Biippr (x)
E Oyya
(C))
Σ Lo(x.)
N
(11)
For readability, we denote E' to be the objective in Equation without the regularizer.
Next, take closer look at the derivate of the expectation over adjacency matrices. Note that have defined the adjacency matrix distribution ph, (C) =I= II Iي= (71))o(01)) with Cj representing the edge X, - Xj. Since parameter "ay only influences the likelihood of the edge X, + X; and no other edges, we can reduce the expectation to a single binary variable over which we need differentiate the expectation:
0 Oyk
N Σ
Epy,o(C) (C)
Σ (x.))
Epy,o(C uz)
EppolCha Oyss
(x,)))
(12)
16

---

Published as a conference paper at ICLR 2022
Intervention
(
Estimato ln Wallhoods and average
Data batch
Lx +x,+xx(X2
C.
17112
Sample grapha
Datermine grastents
Lx, Lxx (X2)
((
E
Xxx (X2)
Graph parameters
LxxssX((22)
C
2012
C
Figure 5: Visualizing the gradient calculation for the incoming edges of X2 in an example graph with three variables. The intervention is being performed on X1, and the data is used to calculate the log- likelihood estimates under the three randomly sampled graphs: Loi (X2), Lea (X2) and Lea(X2) (X2). Those terms are assigned to the Monte-Carlo estimators for Lx.+x, (X2) and Ax,4x, (X2). and finally used to determine the gradients for y and O. The same process is performed for X3 well.
where (CRi) = (rta) a(0s).. The first expectation over P-, e (C'+s)) ) is independent of 7ie as we have defined the adjacency matrix distribution be a product of independent edge probabilities.
The log-likelihood estimate of a variable, (x;),,, depends on the adjacency matrix column C, which represents the input connections to the node X,, All other edges have no influence on the log- oogliikel estimate of X, Hence, the parameter only influences Lo (X)). and thus we can reduce the sum inside the expectation to:
Σ Eo(x))
① Onma o(Cm) (oo,(i)))(((X))) [Xi()
E Oyka
P., a (Ct))
=
(13)
The REINFORCE trick is simple method to move the derivative of discrete distribution inside the expectation. Applied to our situation, we obtain:
(o(Xi) Ölog ((,)) Oyke
Onk Ep,,o [Xo(Xi) = Epy,0((2e))
(14)
This leaves us with cases in the expectation: Cw 0 and Cst In other words, we need to distinguish between samples of O' where we have the edge X₁ + X1, and where we do not have such an edge (Xx f X1). Thus, we can also write the expectation as a weighted sum of those two cases:
Ep, (x)) (logp),((Caa) Tym
ò @log (T)a)) a(OR))
o (Tw)) (0m)) (xx-xx(X)))
Orine
11-o o ((b))) (X+) Ølog (1 a(Ta) (fw))
Oyks
(15)
We use Xx (Xt) to denote the (expected) negative log likelihood for X₁ under adjacency matrices where we have an edge from X₃ X
Lx.-xx (Xt) E. Cu: (Cx(
(16)
Lx,Ax (Xt) Ep,, (C (Xi))
(17)
The final step is to solve the two derivative terms Equation 14. This done as follows:
a log (Tw)) (00a) 2 log log or (0ka)
a(T))
(18)
Oykl
Oy
@ log (1 (vva) σr (Oxz)) (Tm) (1 o (TT))) (0(0))
(19)
Oyk
((T)) (Baa)
17

---

Published as a conference paper at ICLR 2022
Gradient estimation variance y-parameters
0.3
Estimator W-REI NFORCE
Gradient estimate
0.2
ENCO
Y5.4
Y5.6
Vs.7
Parameter
Figure 6: Plotting the gradient estimate variance of ENCO for three variables of y compared to previous REINFORCE- approach by Bengio et al. (2020) on an example graph with 100. The gradients have been scaled to match in terms of averages. We can see clear reduction in variance with the gradient estimator of ENCO allowing us to use lower sample sizes.
Putting these results back in the original equation and adding the sparsity regularizer, we get:
0% -= o7ks) (1- ((1))) ((0;)) Bx.c., (xx,-x, Lx, (X₃) + Asparse| (20)
To align the result with the gradient in Section 3.3, we switch the index notation from kt, to j again. The expectation Ex,c-y a short form for the expectations Bi)p((1) (T) BB, Ex. is) From this expression, we can see that the gradients of "i are proportional to the difference of the expected negative log-likelihood of X3 with having edge between X, - X,, and the cases where X; A Xj. The sparsity regularizer thereby biases the difference towards the no-edge case. The value of Vij and 0g only scale the gradient, but do not influence its direction.
In order to train this objective on a dataset of interventional data, we can use Monte-Carlo sampling to obtain an unbiased gradient estimator. Note that the adjacency matrix samples to estimate Lx.+xX (X₂) and Ax,AX, (X₃) are not required to be the same. For efficiency, we instead sample K adjacency matrices from (C), evaluate the likelihood of batch X under all these graphs. Afterwards, we assign the evaluated samples to one of the two cases, depending on Cij being zero or one. This way, we can reuse the same graph samples for all edge parameters 7. We visualize the gradient calculation in Figure 5. In the cases where we perform intervention X,, we do not optimize Tij for this step and set the gradients to zero. The same holds for gradient steps where we do not have any samples for one of the two oog-likelihood estimates.
A.1.1 COMPARISON TO PREVIOUS GRADIENT ESTIMATORS
As discussed in Section 3, previous work similar structure learning methods (Bengio et al., 2020; Ke et al., 2019) relied on different estimator. In terms of derivation, the main difference is the continuation from Equation 14 on. In proposed method, we write the expectation the of two terms that can independently be approximated via Monte-Carlo sampling. In comparison, Bengio al. (2020) proposed directly apply a Monte-Carlo sampler to Equation 14, and apply an importance sampling weight to reduce the variance. This estimator is also used in the method SDI (Ke et al., 2019) which we have experimentally compared our method.
Figure compared the gradient estimator in terms of standard deviation. The gradient estimator of ENCO achieves a 10 times lower standard deviation compared to Bengio et al. (2020) making much more efficient. Since the estimator by Bengio et al. (2020) is biased and has a different mean, we have scaled both estimators to have the same mean. Specifically, we have applied ENCO to r random graphs from our experiments on synthetic graphs (see Section 4.2), and evaluated 64, 000 sampled adjacency matrices in terms of log-likelioood estimates. These 64, 000 samples are grouped into sets of K samples which we have used to estimate the gradients of 7. We evaluated different values of K, from K = 20 to K = 4000, and plotted the standard deviation of those estimates in Figure 2. We have also visualized three examples as violin plots in Figure 6 that demonstrate that despite both estimators having a similar mean, the variance of gradient estimates is much higher for Bengio et al. (2020).
18

---

Published as a conference paper at ICLR 2022
To verify that the improvement of ENCO is not just because of the gradient estimators, have performed an ablation study with ENCO deploying the gradient estimator of Bengio et al. (2020) in Appendix D.3.
A.2 LOW-VARIANCE GRADIENT ESTIMATOR FOR ORIENTATION PARAMETERS
To derive the gradients for the orientation parameters 6, we can mostly follow the same approach as for the edge existence parameters 7. However, we have to keep in mind the constraint Ost -0s which ensures that the orientation probability sums to one: ((0x) + a (0ux) 1.
To determine the gradient of the likelihood term, can separate the two gradients of 0ut and Oyk- This is because 0µt only influences the expectation over (X())), while 0ix concerns Lx(Xs). We can follow Equation 1I to Equation 20 of Section A. by swapping Out and Tea. For the derivative through the expectation, we obtain the following gradient:
① 00w p .()) [Xi() =0((0)) (Tma) (xx,-xx))xXXxxxxxx)xxxxx Cxexx (21)
Since we have the condition that Ou = -0u. the full gradient for On would therefore consist of the gradient above minus the gradient of Equation 21 with respect to 0zz. However, discussed in Section 3.3, the orientation of an edge cannot be learned from observational data in this framework. Hence, we only want to use the gradients of 0µr if we intervene on node Xg. which gives us the following gradient expression:
0 Ż = (()) ((s)) 00,,
[xx-xx(XX) Lx,Ax, (X²)]
(22)
plxx,) . o(35)) Exxx,
(X:) Cx,xX,(X,)
To align the equation with the one in Section 3.3, we swap the indices k. I with i, again. The first line represents cases where we have an intervention on the variable X;, while we have it over interventions on the variable X, the second line. The two terms are weighted based on the edge existence likelihood a ((as) and (15i) respectively, and the likelihood of performing an intervention on X, or X,. our experiments, we use a uniform probability across interventions on variables, but could design heuristics that selects the
emphasize that this is not strictly required. Moreover, intervention to update the parameters on with the aim of increasing computational efficiency, The gradient estimators presented in Equation 22 would still be valid in such a
We clarify that we do not consider the gradients of 0:j with respect to the edge regularizer. This is done for two reasons. Firstly, the orientation parameter models only the direction of the edge, not whether it exists or not. The regularizer would increase Oiz if the edge existence for the opposite direction would be greater than for the direction from X, to Xj. i.e.ii < 7ji. and decrease 0.j if we have Ty > Tys- However, the orientation should only model the causal direction of an edge. Hence, we do not gain any value from causal perspective when adding the regularizer to the gradient. Secondly, the regularizer would require us to take additional assumptions for guarantecing the discovery of the true graph upon convergence. In experiments with using a regularizer in the 0 0-gradient, we did not observe any difference to the experiments without the regularizer.
We note that the orientation parameters are considered be pairwise independent. In other words, Oij and 0g are considered independent parameters if // k. / and j / I. Global order distributions such as Plackett-Luce (Plackett, 1975; Luce, 1959) can be used to also incorporate transitive relations. However, those require high variance gradient estimators and struggled with chains in early experi- ments. The pairwise orientation parameters provide much easier optimization while still providing convergence guarantees for the full intervention setting.
A.3 TRAINING LOOP
Finally, we give an overview over the full training loop in Algorithm 1. The distribution over interventions p(T) is set to a uniform distribution for all our experiments. However, the distribution
19

---

Published as a conference paper at ICLR 2022
can also be replaced by heuristic which selects interventions to increase computational efficiency. To keep the convergence guarantees, p(T) would have to guarantee a zer- probability to pick any variable. In experiments, experienced that the Adam optimizer (Kingma & Ba, 2015) speeds up the convergence of the parameters y and while not interfering with the convergence guarantees in practice.
Algorithm 1: Learning algorithm of ENCO
Data: variables {X,,.., XN): J: observational dataset Dobs: interventional datasets Din (i) for sparse, perfect interventions on all variables; distribution over interventions p(I) Result: graph structure corresponding to the cuasal relations among variables Initialize 7 = 0
for number of epochs do
/+ Distrimutien sitting for iterations do
X ~ Dobs
for <- 1 to N do
M1 ~ Ber(or(rus) (B(1i)
log fø. (X;; M-i @ X-i)
6 < Adam @0 V L)
end
end
/+ Graph fitting
for
iterations do
/* Sample an intervention
i ~~ppl p(T) with intervention target Xx
X (fm
/+ Evaluate multiple graph samples for gradient estimator
for k = to do
((k) Ber(o(A) (
for i = N do
Law t log fo. (Xzi (X C( x_,)
end end
/+ Update paraneters to N do
for
for j 1 to N where j + and j & L do /* Considering edge X; + X,
DA ot (X,) ΣK d(t)
Cx,-x, (X₂)
Lx,AX, (X₂) ΣK (1-c ). (X,)
Σk_ (1-ch()
Yi -- Vij - o ('(iss) (0 [xx, (X₂) LxiAx, (X,) + Aspane)]
if then
0ij Oij d' (0.g) a(7))) [xx,- (X,) Lxxxxxx
0₂4 02j
end
end
end
end
end
return G= (V, E) with V = X and
E ((Xi, i E [1, N] and A and (vx)) >> 0.5 and a (0,))) > 0.5)
20

---

Published as a conference paper at ICLR 2022
B CONDITIONS FOR CONVERGING TO THE TRUE CAUSAL GRAPH
The following section gives an overview and proves the conditions under which ENCO converges to the correct causal graph given sufficient time and data. We emphasize that we provide conditions here for which no local optima exist, meaning that if ENCO converges, returns the correct causal graph. This is a stronger statement than showing that the global optimum corresponds to the true graph, since wradient-based algorithm can get stuck in local optimum. We will discuss the conditions for the global optimum in Appendix B.2.5.
To make the proof more accessible, we will first discuss the assumptions that needed for the guarantee, and then give sketch of the proof. The proof will first assume that we work in the data limit, i.e. have given sufficient data, such that can derive conditions that solely depend on the causal graphical model. In Appendix B.2.3, we extend the proof to the limited data setting.
B.1 ASSUMPTIONS
We given dataset of observational data from the joint distribution p( X). Additionally, we have N interventional datasets for N variables where in each intervention a different node intervened on (the intervention size for each dataset is therefore 1).
Assumption 1
Assumption A common assumption in causal structure learning is that the data distribution over all variables p(X) is Markovian and faithful with respect to the causal graph we are trying to model. This means that the graph represents the (conditional) independence relations between variables in the data, and (conditional) independence relations in the data reflect the edges in the graph. For ENCO, faithfulness is not strictly required. This is because we work with interventional data. Instead, we rely on the Markov property and assume that for all variables, the parent set pa(X;) reflects the inputs to the causal generation mechanism of X, This allows us to also handle deterministic variables.
For this proof, we assume that all variables of the graph are known and observable, and no latent confounders exist. Latent confounders can introduce dependencies between variables which are not reflected by the ground truth graph solely on the observed variables. We discuss the extension of latent confounders in Section 3.5 and Appendix B.3.
Assumption 3
ENCO relies on neural networks to determine the conditional data distributions p(),)...) Hence, for providing guarantee, assume that in the graph learning step the neural networks have been sufficiently trained such that they accurately model all possible conditional distribution p(X.)...). In practice, the neural networks might have a slight error. However, as long as enough data, network complexity, and training time is provided, it is fair to assume that the difference between the modeled distribution and the true conditional is smaller than an arbitrary constant c, based on the universal approximation theorem (Hornik et al., 1989). For the limited data setting, see Appendix B.2.3.
Assumption 4
Assumption 5 We are given a sufficiently large interventional dataset such that sampling data points from it models the exact interventional distribution under the true causal graph. This can be achieved by, for example, sampling directly from the causal graph, or having an infinitely large dataset. For the limited data setting, see Appendix B.2.3.
B.2 CONVERGENCE CONDITIONS
The proof of the convergence conditions consists of the following three main steps:
Step 1 We show under which conditions the orientation parameters 0,5 will converge to +oo, i.e. (0;;)) + X; an ancestor of X,. Similarly, if X; descendant of X,, the parameter Oz5 will converge to -o0, i.e. o (0.y) 0.
Step 2 Under the assumption that the orientation parameters have converged as in Step 1, we show that for edges in the true graph, "ih will converge to 1. Note that we need to take additional assumptions/conditions with respect Asparse here.
21

---

Published as a conference paper at ICLR 2022
Step 3 Once the parameters Vay and Ogy have converged for the edges in the ground truth graph, show that all other edges will be removed by the sparsity regularizer.
The following paragraphs provide more details for each step. Note that causal graphs that do not fulfill all parts of the convergence guarantee can still eventually be recovered. The reason is that the conditions listed in the theorems below ensure that there exists no local minima for 0 and ry to converge in. Although local minima exist, the optimization process might converge to the global minimum of the true causal graph.
Theorem B.1. Consider the edge X; + X3 in the true causal graph. The orientation parameter Osj converges to a (0:s) = if the following two conditions are fulfilled:
(1) for all possible sets of parents of Xj excluding X,, adding X; improves the log- log-likelihood estimate of Xj under the intervention on X;, or leaves it unchanged:
Vpa(X;) S X-is" EIxx X X llog ((x)) -oog (XxXx) 0
(23)
(2) there exists a set of nodes (i(X)), for which the probability to be sampled as parents of Xj is greater than O, and the following condition holds:
for
Xx) S X Eix, x [oogp(X (xx)) 1ogg(Xxlba(Xx).
- logp(X, pa(Xx))) 0
(24)
Proof. Based on the conditions in Equations 23 and 24, we need to show that the gradient of 0. negative in expectation, independent of other values of y and O. For readability, we define the following function:
T(x: = Eix, X.C_M (Cxe-x - (xxxx(
(25)
Hence, the gradient of Oig can be written as:
0 L= (p(Tx)) a(1))) X(X,, p(Txx) ((9))) (X))
(26)
08:3
Looking at the gradient of O.g in Equation 26, the conditions correspond to T(X,, X,) being smaller or equals to zero. Note that the sign is flipped because in T(X,, X,,, we have negative licelioos represented by Lx2 -xx Ex,+xx(XX (Xi). while in Equations 23 and 24, we have log-likelihoods. Further, the other factors of o' (0.g), (1ig) and p(lx) are all limited in the range of (0, 1) meaning that the sign of the gradient is solely determined by T(X;, and T ((X,,X)) If T(X,, X,) T(X,, X;) is smaller than zero, then the gradient of 0zg is negative, i.e. iee.ncreasiig θig-
First, we look when T(X,, X,) O. The condition in Equation 23 ensures that conditioning X; on a true parent X; when intervening on X; does not lead to a worse log-likelihood estimate than without. While this condition might seem natural, there are special cases where this condition does not hold for all variables (see Section B. B.2.4). The second condition, Equation 24, guarantees that there is at least one parent for which T T(2,, X,) is negative. Therefore, in expectation over all possible adjacency matrices (C), T(X,, X,) is smaller than zero if the two conditions hold.
To guarantee that the whole gradient of f.y is negative, we also need to show that for interventions on X,, T(X,. X,) can only be positive. When intervening on X,. X; and X; become independent as the edge X; - Xj is removed in the intervened graph. A distribution p(X.|Xx..>. relying on correlations between X; and X, from observational data cannot achieve better estimate than the same distribution when removing X,. This because the cross entropy is minimized when the sampled distribution, in this case p( X(),), is equal to the log- likelihooo estimator (Cover & Thomas, 2005):
Σ P(X: = 2z) = =x) S Σ p(X: = x;) log ( = =2 = zz)) 21,2,
(27)
The only situation where X; and X; can become conditionally dependent under interventions on X; if X; and X, share a collider Xx. and X, is being conditioned on the collider X₁ and X,. However, this requires that bæ has negative gradients, i.e. 6%z increasing, when intervening on Xx. This cannot be the case since under interventions on X3, X, and X₁ become conditionally independent, and the correlations learned from observational data cannot be transferred to the interventional setting. If X
22

---

Published as a conference paper at ICLR 2022
and X; again share collider, we can apply this argumentation recursively until node X, does not share a collider with X₁. The recursion will always come to an end as we have a finite set of nodes, and the causal graph is assumed to be acyelic.
Thus, if the conditions in Equations 23 and 24 hold for an edge X; + X5 in the causal graph, we can guarantee that with sufficient time and data, the corresponding orientation parameter 0iz will converge to o (0.g)) 1.
Theorem B.2. Consider a pair of variables X;, f which X; an ancestor of X; without direct edge in the true causal graph. Then, the orientation parameter of the edge X; + X, converges to a (00,g) ) = 1 if the same conditions as in Theorem B.1 1 hold for the pair of X3.
Proof. To show this theorem, we need consider two cases for pair of variables X; and X,; X; and X3 are conditionally independent under a sampled adjacency matrix C, or X; and X, are not independent. Both cases need to be considered for an intervention on X; with the log- likelihoood estimate of X,, and an intervention on X, with the log-likelihood estimate of X;
First, we discuss interventions on X₃ If under the sampled adjacency matrix C, X, is conditionally independent of X,, the difference in the log- liklliieehood estimates T(X,, X,) is zero in expectation. The variables can be independent if, for example, the parents of X; are all parents of the true causal graph. If X; is not conditionally independent of X; the conditions in Equations 23 and 24 from Theorem B. ensure that X, has, in expectation, a positive effect on the log-likelihood estimate of Xy. Thus, under interventions on X,, the gradient of 0t is smaller or equals to zero, iee.neresssss 0yy
Next, we consider interventions on Xy. If under the sampled adjacency matrix X, is conditionally independent of X,. the difference in the log- liglliiel estimates T(X;, X,) is zero. The variables can be independent if X; is conditioned on variables that d-separate X; and X, in the true causal graph. For instance, having the children of X; as parents of X, creates this scenario. However, for this scenario to take place, one or more orientation parameters of parent-child or ancestor-desseedant pairs must be incorrectly converged. In case of parent-child pair X,, Xx. Theorem B. shows that or (0xx) will converge to one removing any possibility of a reversed edge to be sampled. In case of an uncestor-descendant pair X,, X₁, we can apply a recursive argument: as X₁ d-separates X; and X,. X, must come before X, in the causal order. If for the gradient Ou. we have a similar scenario with X, being conditionally independent of X,. the same argument applies. This can be recursively applied until no more variables except direct children of X; can d-separate X; and X,. In that case, o (0a)) will converge to one, which leads to all other orientation parameters to converge to one as well. If X; is not conditionally independent of X,, we can rely back on the argumentation of Theorem B. when we have an edge X, - Xy: as in the intervened causal graph, X, and X, are independent, any correlation learned from observational data can only lead to worse log- ligeliilldoood estimate. In cases of colliders, we can rely on the recursive argument from before. Thus, under interventions on Xj. the gradient of 0,y must be smaller equals to zero in expectation, i. incceeesss 0yy
Therefore, we can conclude that o (Bis) converges to one for any ancestor-descendee pairs X,, X under the conditions in Theorem B.I.
Theorem B.3. Consider an edge X; + X; in the true causal graph. The parameter Yig converges to (14)) = 1 if the following condition holds:
min (1) [logp][[geeXXee X.) - log ((xl) >
Asparse
(28)
where (X;) is the set of nodes excluding X; which, according to the ground truth graph, could have an edge to X, without introducing a cycle, and pI_, (I) refers to the distribution over interventions (I) excluding the intervention on variable X3.
Proof. To show this convergence, we assume that the orientation parameters have converged cor- responding to Theorem B.1 and B.2. The parameter "a converges a (71y)) 1 if its gradient, Byyy e is negative independent of other values of y and orientation parameters 0 that are not included in Theorem B.1 and B.2. The gradient of Tiy includes an expectation over adjacency matrices (C)) Based on the converged 0 0-values, only need to consider sets of nodes as parents for X, that contain parents, ancestors, or (conditionally) independent nodes according to the ground truth graph. This sets of parents is represented by gpa, (X;). Among those remaining parent sets,
23

---

Published as a conference paper at ICLR 2022
need to ensure that for any such set, the gradient is negative. The condition in Equation 28 corresponds to the inequality -_ L < O since the term the left represents the log- -iglikklihood difference -x,-x, (X,) Ax,, (X,) in the gradients of Yi in Equation 20 with a flipped sign. For readability and better interpretation, Agparse has been moved on the right site of the inequality. This is possible as Asparse is independent of the two expectations in Equation 28. If the inequality holds for all parent sets pa), the gradient of "/j be guaranteed to be negative in expectation, independent of the other values of ry Since the distribution over parent sets pa) depends on other values of y, the condition in Equation 28 ensures that even for the parent set with the lowest log-likelihood difference, it is still larger than Asparse. If this condition holds, then the gradient of Tht zero independent of other values of 7.
- will be smaller than
The condition in Equation 28 introduces a dependency between convergence guarantees and the regularizer parameter Aspurse- The lower we set the regularization weight Axuase, the more edges can guarantee recover. If the regularization weight is too high, false negative edge predictions. If the regularization weight is set very low, we take a longer time to converge as it requires lower gradient variance or more update steps, and is more sensitive in limited data regime. Nonetheless, if sufficient computational resources and data is provided, any value of Apurre > can be used.
eventually obtain
Theorem B.4. Assume for all edges X; - X in the true causal graph, (00;)) and o ((ig)) have converged to one. Then, the likelihood of all other edges, i.e.)(Or)) (0w),, will converge to zero under the condition that Apprre >> O.
Proof. If all edges in the ground truth graph have converged, all other pairs of variables X,, Xx (conditionally) independent in the graph. This statement follows from the Markov property of the graph and excludes ancestor-descendant pairs X,, X,. The possibility of having edges from descendants to ancestors has been removed by the fact that the orientation parameters Oij have
converged according to Theorem B.2. Thus, for those cases, we already have the guarantee that (0,,)) (B.,) converges to zero.
For conditionally independent pair X1, XR, the difference of the log-likelihood estimate in the gradient of "/u, Lx,.xx (Xx) Lx,AX, (X&), zero in expectation since independent nodes do not share any information. Thus, the gradient remaining is:
i C d'(w)) o (0ix))
Asparse
(29)
Oykk
Since the gradient is positive independent of the values of Tue and Ouc Tue will decrease until it converges to (71x))= 0.
Hence, if nuk decreases for all pairs of (conditionally) independent variables X1, Xx in the ground truth graph, and a (0(x)) converged to zero for children and descendants, the product a (nw)) a (00)) will converge to zero for all edges not existing in the ground truth graph.
For graphs that fulfill all conditions in the Theorems B. and B.4, ENCO is guaranteed to converge given sufficient data and time. The conditions the theorems ensure that there exist no local minima or saddle points in the loss surface of the objective in Equation 2 with respect to 7 and O.
Summary We can summarize the conditions discussed above as follows. Given a causal graph G with variables XI,...,N,, and sparse interventions on all variables, the proposed method ENCO will converge to the true, causal graph G. if the following three conditions hold for all edges X; X; in the true causal graph G:
1. For all possible sets of parents of X, excluding X,, adding X; improves the log- loglliked estimate of X2 under the intervention on X,, or leaves it unchanged:
*fa(X,) C X-1 Bix, xx (aaap(XX X,) log p( X Ipa
(30)
2. There exists a set of nodes pa( (X;), ), for which the probability to be sampled as parents of X; is greater than O, and the following condition holds:
7pa(X,) C X Exx,,x [logr(X; (X,), X,) logg(Xxlpa(xx))) 0
(31)
24

---

Published as a conference paper at ICLR 2022
X1
X2
X3
Figure 7: Example graph for which we check the convergence conditions described in Appendix B.2. The conditional distributions are given in Equation 33 to 35.
3. The effect of X; on X3 cannot be described by other variables up to Asuaree
min E i~pr ,00 (x) log (X, X;) log p (xppp( paCgpu
Apparse
(32)
where gpa,(X) is the set of nodes excluding X, which, according to the ground truth graph, could have edge Xj without introducing cycle.
Further, for all other pairs X,, Xj for which Xj is descendant of X,, conditions (1) and (2) need to hold as well.
B.2.1 EXAMPLE FOR CHECKING CONVERGENCE CONDITIONS
In the following, we will provide a walkthrough for how the conditions above can be checked on a simple example graph. For further details the precise calculations, provide Jupyter Notebook that contains all calculations in this example'.
Suppose have graph with 3 binary variables, X1, X2, X3, with the causal graph being X1 X2 + X3, i..., small chain. For simplicity, let us assume that the true, conditional distributions are the following:
((X1) = Bern(..7)
(33)
p(X)(x)) = [X1
with prob. 0.6
(34)
[X,@] with prob. 0.4
[X2
with prob. 0.2
X(xx(x2) =
(35)
(x2 @
with prob. 0.8
In other words, X₂ equals to the value of X, with a probability 0.6, and the opposite binary value otherwise. Similarly, X3 is equals to the value of X2 with probability of 0.2, and the opposite binary value with a probability of 0.8. Therefore, the sample with the highest probability in this joint distribution would be X1 1, 1, X3 O. Further, assume that all interventions replace the respective conditional distribution by uniform distribution, i.., Pix, (X,) = Bern(0.5). Next, will check the conditions for the edges in G, ie., X1 X2 and X₂ X3, and the remaining ancestor-descendamt pair X1,X3.
Edge X1 + X2:
Condition 1: the possible parent sets that exclude X, and Xa X3 and the empty set. For the empty set, we get:
Exx, x (logp(X2/Xx) logp( 0.023
For conditioning on X3. we obtain:
E Exx, x
logr(X2|XX)) ~ 0.015 2 0
Since both values are greater than zero, condition is fulfilled for X; - X2.
Condition 2: is already fulfilled by the equations in condition since all parent sets have difference greater than zero.
'The calculations can be found in the notebook called convergence_guuraateee. ipyab, https://gt thub. thwb.com/ph..... com/ph1 1pynb.
25

---

Published as a conference paper at ICLR 2022
Condition 3: the set gpa, (X2) is the empty set since X3 is descendant of X2, and other nodes exist in the graph. Thus, the parent set minimizing the expression on the left can only be the empty set, and we calculate it as follows:
E i~pr 2 (x) -og ~ 1/2 0.023 1/2 0.017 = 0.020 > Asparse
1x₃
Ix,
with assuming (I) being the uniform distribution, and excluding Ix, since we do not update "/2 in this case. Hence, as long as Aqurse is smaller than 0.02, the condition is fulfilled.
Edge X2 + X3:
Condition 1: the possible parent sets that exclude X2 and X3 are X1 and the empty set. For the empty set, get:
Exx, x (logp(Xs) - log logg((X)) 0.194 2
For conditioning on X1, we obtain:
E Exx, x [logp(X3 (x3 X2, X1) logg 0.200 2
Since both values are greater than zero, condition is fulfilled for X₂ - X3-
Condition 2: is already fulfilled by the equations in condition since all parent sets have difference greater than zero.
Condition 3: the set gpaz(X3) contains the variable X₁ since
introduce edge X1 - X3 without introducing acyclicity in the true, causal graph. Thus, we need to compare parent sets for finding the minimum of the left-side term: X, and the empty set. First, we consider the empty set:
E
1/2-
1/2-0.194 0.194 Aqparse
1x
Ix2
Again, we exclude Lx, since we do not update 7/23 in this case. The second case considers Xi as additional parent set pa:
[logr Naa X1) log (xx R 1/2. 0.186 1/2 0.200 E imps ,(n) (p)
fx
1xg
0.193
Asparso
The minimum of both values is 0.193. Hence, the edge X₂ + X3 can be recovered if 0.193 Asparse-
Ancestor-descendant pair X1, Xz:
Condition 1: the possible parent sets that exclude X, and Xg are X2 and the empty set. For the empty set, we get:
Exx, x [logp(Xs/Xx) logr(Xa)) 0.008
For conditioning on X2, we obtain:
E Exx, x [log Nogp(X3/Xx X2) log X2)) 0 0
The difference is zero because X₃ is independent of X1 when conditioned X2: p(Xy] X2) = p(Xx Since both values are greater or equals to zero, condition 1 is fulfilled for the pair X1, X1,
Condition 2: from condition 1, can see that the parent set of pa(Xa) being the empty set the only option that fulfills the condition being greater than zero. Since we start the optimization process with an initialization that assigns non-zero probability to all possible parent sets, follows that pa(X;) being the empty set has probability greater than zero throughout the optimization process. Hence, condition 2 is fulfilled well.
26

---

Published as a conference paper at ICLR 2022
Summary: in conclusion, for the discussed example, we can guarantee that ENCO converges to the correct causal graph if Aspase < 0.02. To experimentally verify this results, we applied ENCO on this graph with two hyperparameter settings for the sparsity regularizer: Axparse 0.019 and Xygarse = 0.021. We considered a very large sample size, more specifically 10k per intervention and 100k observational samples, to simulate the data limit regime. For Aupanse 0.019, ENCO was able to recover the graph without errors while for Apaasse = 0.021, the edge X1 + X2 was, as expected, missed. This verifies the theoretical result above with respect Asparse. Note that if the condition is not fulfilled by selecting too large sparsity regularizer, this does not necessarily mean that ENCO will not be able to recover the graph. This is because we consider the 'worst-case' parent set in condition 3, while this case might not be in the true causal graph to which the other edges converge.
B.2.2 INTUITION BEHIND CONDITION AND 2
As mentioned in the text, condition and 2 of Theorem 3.1 ensure that the orientation probabilities cannot converge to any local optima. Since the conditions explicitly involve the data distributions and implicitly the gradient estimators, we provide below an assumption from data generation mechanism perspective as an alternative, that ensures condition 1 and 2 to be satisfied.
Firstly, we assume that ancestors and descendants are not independent under interventions on the ancestors. Note that there can exist graphs where the ancestors are independent of descendants, for instance in linear Gaussian setting when the ancestor has weight of zero on the descendant. However, those graphs, violating faithfulness, are impossible to find for any causal discovery method since the variables are independent under any setting. In terms of condition I and 2, would imply that the inequality is always zero.
Next, we show that under the previous assumption, local optima of the orientation probabilities can only occur in the following structure: for an edge X; + Xj. there exist one or more parent(s) of Xj sharing a common confounder X₁ with X,, where X₁ is not a direct parent of Xj. An example of this structure is the following: Xi - X2, X3; X2, X3 - X4 where the orientations of the edges X₂ + X4 and X3 + X1 could have a local optimum. This statement can be proven as follows by using the do-calculus (Pearl, 2009). Suppose a graph that includes the three variables X1 X2. X3 with X, - X2, X3 - X2, and X2 having no parents besides X, and X3. X; and X3 do not share a confounder, then, from do-caleulus, we know that p(X² do X1 =1)) p(X2| = 21) and p(X, (x3 = ra)) = p(X, X3- zz). Furthermore, since the conditional entropy of a variable can only be smaller or equals to the marginal, i.e. (X)) Z 11 (X]Y), estimating X2 under interventions on X1 can only be improved by conditioning on X1, and similarly for X3. Thus, condition 1 strictly fulfilled when parents do not share a confounder under the previous assumption of no independence in all possible settings. Now, consider the situation where X1 and X3 share a common confounder. Then, from do-calculus, we can state that there can exist a parameterization of the conditional distributions for which p( X2| X; = 21)) * p(X2 X2 21). Under this setting, we cannot guarantee that condition is always fulfilled. However, whether this parent-confounder structure above actually leads to local optimum not depends on the distributions, which condition 1 models. Intuitively, this requires the mutual information between the two or more parents to be very high, and the initial edge probabilities of those edges to be very low. Further, as the results show, this combination of events is not very common in practice, meaning that as long as the ancestor and descendant are not independent under the interventions, we usually converge graph with the correct orientation.
Besides, if the confounder X, of X, and X3 is parent of X2, then the local optimum would disappear with learning that edge since p(Xgldo(X, = x,) = p(Xx[X, = Iι, Xz). In conclusion, for many of the graph structures like chain, bidiag, collider, full and jungle, this shows that there does not exist any local optima for the orientation probabilities. Only for the certain structures of confounded parents, there may exist local optima that depend on the specific distribution parameterization.
B.2.3 LIMITED DATA REGIME
Assumption (3) and (4) are taken with respect to the data limit such that the conditions derived in the next section solely depend on the given causal graphical model. However, in practice, often have a limited data set. The proof presented for the data limit is straightforward to extend to this setting with the following modification:
27

---

Published as a conference paper at ICLR 2022
X1
X2
X1
X2
X3
X3
(a) Fork example
(b) Fully connected example
Figure 8: Causal graph structures for which, under specific parameterization of the conditional distributions, the conditions for guaranteeing convergence can be violated.
The conditional distributions p(.)..) are replaced by the conditional distributions that follow from the given, observational data.
The expectations over the interventional data (X) is replaced by the joint distribution over samples given for the intervention! î.
Theorem B.1 and B.2 for the edge X; - Xj extended as follows:
(3) For all possible sets of parents of X, excluding Xj, adding X; does not improve the log-likelihodd estimate of Xi under the intervention on Xj. or leaves it unchanged:
vpa(X.) S C--ss Eixx, xX (x,pp(X//a(xx))...
log pog)((X.pa(Xx))] <
(36)
This condition is the inverse statement of Equation 23, in the sense that consider interventions on the child/descendant X,. In the data limit, this naturally follows from Equation 23 and Equation 24, but in the limited data regime, we might have violations of Equation 36 due to biases in our samples. Violations of Equation 36 are the cause of ENCO predicting cyclic graphs as seen in Section 4.2.
Finally, Theorem B.4 does not necessarily hold anymore since noise in our data can lead to overestimation of edges. Thus, we add the following condition:
(1) For all pairs of variables X5, X3 for which there exists no direct causal relation in the true causal graph, and Xj not being the ancestor of Xi, the following condition has to hold:
min (X,))xx)
E (i) BB [log (ogg((X)lpaXXXX lo) (Xppp((XXp))))))) pa)] < Ayuarse
(37)
where (X;; is the set of nodes excluding X; which, according to the ground truth graph, could have an edge to X, without introducing cycle.
This condition ensures that no correlations due to sample biases introduce additional edges in the causal graphs
If the conditions discussed above hold with respect the given observational and interventional dataset, we can guarantee that ENCO will converge to the true causal graph given sufficient time.
B.2.4 GRAPHS WITH LIMITED GUARANTRES
Most common causal graphs fulfill the conditions mentioned above, as long as a small enough value for Axparse is chosen. Still, there are situations where we cannot guarantee that ENCO convergences the correct causal graph independent of the chosen value of Xyparse Here, we want to discuss two scenarios visualized in Figure under which the guarantees fail. Still, want emphasize that despite graphs not fulfilling the conditions, ENCO might still converge to the correct DAG for those as the guarantee conditions assume the worst-case scenarios for 0 and ry in all situations.
The first example we discuss is based on a fork structure where we have three binary variables, [X1,X2, Xz), and the edges X, - X3 and X₂ - X3 (see Figure 8a). The parameterization look at is a (noisy) XOR-gate for X3 with its two input variables X,, X2 being independent of each other and uniformly distributed. The conditional probability distribution p( X3. X1, X2) can be
28

---

Published as a conference paper at ICLR 2022
summarized in the following probability function:
if 11 = 0, X2
if X1 = 11, =
x
(38)
c
X1 = -0,X² = =1
1 xf = 1, =
In other words, if 1 A X2, X3 is equals 1 with a likelihood of 1 e. If = X2, X3 is equals 1 with likelihood of c. The issue that this probability table creates the following. Knowing only one out of the two variables does not improve the log likelihood estimate for the output. This is because X2 and Xę are independent of each other, and p (1 = p (x)) is a uniform distribution. Hence, the worst-case parent set in Equation 6 would be the empty set, and leads to an expected difference log-likelihood difference of zero. As Asparse is required to be greater than zero for Theorem B.4, we cannot fulfill the condition for that graph. This means that an empty graph without any edges is a local minimum to which ENCO could converge. Yet, when the edge probabilities are non-zero, will sample adjacency matrices with both input variables being a parent of X3 with non-zero probability. Hence, the log- lig-liiel difference for X2 and X₂ to X3 is unequal zero. Further, this graph is still often correctly discovered despite ENCO not having a convergence guarantee for it. We have conducted experiments on this graph with e = (0.1, 0.2, 0.3, 0.4, 0.45] using a sparsity regularizer of >sparse 1e-4, and in all cases, ENCO converged to the correct, acyclic graph. Note that values close to 0.5 for e are most challenging, because the difference between the true conditional and marginal distribution goes against zero.
The second example we want to discuss aims at graphs that violate the condition in Theorem B. more specifically Equation 23. The graph we consider is fully connected graph with three variables X1, X3 (see Figure 8b). The scenario can be described as follows: if knowing X2 informs the log-likelihood estimate of X3 more about X, than about X₂ itself, an intervention on X2 and sampled graph with the edge X2 + X3 could lead to a worse likelihood estimate of X3 than without the edge. For this scenario to happen, p(Xg|_) must be close to deterministic. Additionally, p(Xx, 11, X2) must be much less reliant on X2 than on X1, such as in the following probability density:
C1
if X1 = 0, = 0
1 G
if X1 = 1, 0
p(X, = 1[x,,,,),,,
(39)
62 if X1 O,
1 (1-
if = 1, =1
The two variables EI, 62 represent small constants close zero. In this case, the graph can violate the condition in Equation 23 since intervening on X₂ breaks the dependency between X1 and X₂. The conditional distribution p(Xs| (22) learned from observational data relies on the dependency between X1 and X2 which can make it to a worse estimate than p(Xs). Note that if the edge X1 + X3 is learned by ENCO though, this will not constitute problem anymore since with conditioning on X,, i.e. p( X3 X2, X1). the edge X₂ - X3 will gain a gradient towards the correct graph. Thus, when y and 0 are not initialized with the worst-case values, the graph with both X1 and X₂ as parents of X3 can be sampled and provides gradients in the correct direction. Further, we did not observe any of these situations in the synthetic and real-world graphs we experimented on.
B.2.5 CONDITIONS FOR THE GLOBAL. OPTIMUM
So far, the discussion focused on proving that the optimization space does not contain any local optima with respect to the graph parameters y and 0 besides the global optimum. If these conditions are violated, ENCO might still converge to the correct solution, since we are not guaranteed to find and get stuck in one of these local optima. Thus, in this section, we provide conditions under which the ground truth graph is the global optimum of the objective in Equation 2. Graphs that fulfill these conditions are very likely be correctly identified by ENCO, but with a suboptimal choice of hyperparameters, initial starting conditions etc., we could return an incorrect graph.
The conditions and proof follow similar structure to those before for the local optima. We first discuss when we can guarantee that the global optima has the same orientation of edges, and then when we also find the correct parent set of the remaining variables.
29

---

Published as a conference paper at ICLR 2022
Theorem B.5. For every pair of variables X;, X3 where X, is parent of X,, the graph G that optimizes objective in Equation 2 models the orientation X; + X5. if there exists an edge between X; and Xj in G, under the following conditions:
Xx and Xj are not independent under observational data.
. Under interventions on X,, X; and X3 are not independent given the true parent set of X
Proof. If x, and X; are independent under observational data, the observational distributions would not identify any correlation among those two variables. Hence, transferring them for any graph to interventional data would have ((.) p(X,|X,, thus making the objective invariant the orientation of the edge, and removing any edge between X; and X; for sparsity.
If X; and X, are dependent, we can prove the statement by showing that modeling the orientation X, + X, will strictly lead to worse estimate under the intervention on X, since the orientation parameters are optimized by comparing the interventions of the two adjacent variables. Under interventions on X; the causal mechanism p( X, /pp Xj)), with (Xx)) being the parent set of the ground truth graph including X,, remains invariant under interventions on X;, and is strictly better than pXx, ppa(X;) (X) for estimating X, due the direct causal relation. Under interventions on X the causal mechanism p pX,, leads to a strictly worse estimate as discussed in Theorem B. 1, since the dependency between X; and X; does not exist in the interventional regime. Hence, the inverse orientation the edge X, + X,, i.e. X3- + X, cannot be part the global optimum. O
Theorem B.6. For every pair of variables X, where X; is ancestor but not direct parents of x, the graph G that optimizes objective in Equation 2 does not include the edge X, X, if the conditions in Theorem B.5 hold.
Proof. To show this statement, we need to consider different independence relations between X; and Xj. First, if X; and X, are independent in the observational dataset given any conditional set, the edge will be removed since any edge between two independent variables is removed for any Asparse > O. The same holds if X; and Xj are independent for interventions on X, and Xy-
If they dependent, p( Xx,.X..... transfers from observational to interventional data on X, since on interventions on X,, the causal mechanism of X, is not changed. Further, when intervening on X,, X; and Xj become independent such that any mechanism p(X,|_X,.... cannot transfer except if X, and X, are independent under interventions on Xj. In this case, the edge will be again removed by the sparsity regularizer. This shows that for any setting, the orientation of the edge X; + X; cannot lead to better estimate than X; > X,, and in case of independence, the edge X, - X, will be removed as well.
follow a similar argument in Theorem B.2. The causal mechanism
Theorem B.7. The graph G that optimizes objective in Equation 2 models the same parent sets for each variable under the following conditions:
The conditions in Theorem B.5 hold.
For any variable X; with its true parent set pa(X.), there does not exist smaller parent set pia X X,, descendants (X,) which approximates the log-likelihood of X; up to Xxparze (Ipa( X,) - (pal) on average.
The regularization parameter Xxparie is greater than zero.
Proof. If the orientations for all edges are according to the ground truth graph the global optimum following Theorem B.5 and B.6, the parent set for a variable X; is limited to those variables which are not descendants of X;. From the ground truth graph, we know that, conditioned on the true parent set pa(X,), ). X; is independent of all other non-descendant variables. Thus, the log-likelihood estimate of X,, i.e. the left part of the objective in Equation 2, is optimized by p(X, (XX,pa(XX)). To show that this is also the global optimum when combining with the regularizer, those parent sets of X, which obtain a lower penalty, i.c. smaller parent sets. The difference between two parent sets pa( (X;) and pa in terms of the regularizer corresponds Aguase ((pa(X.) pal).
need consider
30

---

Published as a conference paper at ICLR 2022
Xi
XI
X1
X
X1
X2
X3
(a) Independent children
(b) Ancestor-doscendane pair
Figure Visualization of the different graph structures we need consider in the guarantee discussion of latent confounder detection. Latent variables X₁ are shown in white, all other variables are observed. (a) The two children of X1, X1 and X2, are independent of each other. (b) X1 an ancestor of X3. and the two variables have shared latent confounder X
Thus, there exists no parent set for which this difference is greater than the penalty for the worse log-likelihood estimate, the true parent set pa(X,) constitutes the global optimum.
B.3 CONVERGENCE CONDITIONS FOR LATENT CONFOUNDER DETECTION
In Section 3.5, we have discussed that ENCO be extended to graph with latent confounders. For this, we have to record the gradients of Ki for the interventional data on X; and all other
interventional data separately. We define "i - Yi + +7g where Tij only updated with gradients
(I) (0)
(0)
from Equation 3 under interventions on X,, and Tif on all others. The score detect latent confounders is:
le lc(
(0) 115
(1-0(11)) (
(1-o
(35 ()
(40)
In this section, we show under which conditions the score Ic( X,) converges one if X, and X, share a latent confounder. We restrict our discussion to latent confounders between two variables that do have any direct edges with each other, and assume that the confounder is not child of any other observed variable, We assume that the causal graph based on the observed variable fulfills all conditions of Theorem B.1 B.4 in Section B.2, meaning that without the latent confounders, the graph could have been recovered without errors. Under those conditions, we can also show that the graph among observed variables with latent confounders is also correctly recovered. This is since the latent confounders only affect Theorem B.4: if X, and X, share latent confounder, they are not conditionally independent given their observed parents. Thus, we can rely on the fact that edges in the true causal graph will be found according Theorem B.1 to B.4. and the edges with latent confounders do not fulfill Theorem B.4.
For all pairs of variables that do not share a latent confounder, Ic X2) converges zero. The edges that are removed in Theorem B.4 converge to (7) (0)) = O which sets le(X,, Xj) K to zero. For edges that have been recovered, state in Equation 24 that the gradient for interventional data must be negative for interventions on the parent. Hence, (7,)1 converges to one which brings Ic)
(2),
to zero again.
For variables that share a latent confounder, we distinguish between two cases that are visualized in Figure 9. In the first case, we assume that X; and X, independent in the true causal graph excluding the latent confounder. This means that an intervention on X, does not cause any change in X,. and vice versa. The second case describes the situation where X; is an ancestor of X₃- The case of X₂ being a parent of X, has been excluded in earlier assumptions as in those cases, we cannot separate the causal effect of X, on X, based on its causal relation and the latent confounder.
In case that the two children of the latent confounder are not an ancestor-descendant pair, we can provide guarantee under the following conditions.
Theorem B.8. Consider a pair of variables X,, X, that share latent confounder Xt. Assume that X, and X2 are conditionally independent given the latent confounder and their observed parents. Further, all other edges in the causal graph have been recovered under the conditions of Theorem B.
31

---

Published as a conference paper at ICLR 2022
to B.4. The confounder score Ic(X, X,) converges one if the following conditions hold:
B Bimpix, (1) (x) ( (a,) logg(X,
Asparse
(41)
Bi-dix, (1) (x) log p( X |pa( (X;), X)) - -log
(esp(X,lpuee -logp(X; X:))] Apparee
(42)
Proof. We need to show that under the two conditions above, or (nl)) ) and o735)) (0)) are guaranteed to
converge to one while ( ( and (?7])²)) converge to zero. The distribution pt.x, (I) represents the distribution over interventions excluding the ones performed on the variable Xx. The two conditions
resemble Equation 28 with the difference that the intervention on the potential parent variable
excluded, and the parent set is the true parent set of the correct causal graph. This is because all other
edges have been correctly recovered, and the two conditions are concerning (>:1))). If the condition
in Equation 41 holds, it corresponds to a negative gradient in (0) following the argumentation in 7sj
Theorem B.3. The same holds for (0) Therefore, ((T5)) and (77 (2)))) are guaranteed to converge to one under the conditions given in Theorem B.8.
For the interventional parameters Vij and ?ja we show that the gradient can only be positive, i.e. decceeceai ) and ?/4e Under the intervention on X,, X; and X, become independent since we assume perfect interventions. In this case, the log-likelihood estimate of X, cannot be improved by conditioning on X,, Hence, the difference Ex,,xX((,)) zero. When further considering the sparsity regularizer Aqurse the gradient of Yij under interventions
Lx, (X,) is greater or equal to
(f) on X; can only be positive, i.e.eerreasing "ij The same argument holds for 7ij Thus, we can conclude that (11))) and (n!!1)) converge to zero.
If the conditions of Theorem B.8 are not fulfilled, ( results in the score lc(X, (. XX,)
(0) and o 10)) might converge to zero. This being zero, but also (7ay) converging to zero. Hence, we do not get any false positive edge predictions as we have seen in the experiments of Section 4.5.
For the second case where X, is an ancestor of Xj, we cannot give such a guarantee because of Theorem B.2, Theorem B.2 states that (f0,) converges to one for ancestor-descendaat pairs. However, a (f5:) is factor in the gradients of 7/j- This means that if ((3z) converges to zero according to Theorem B.2, we cannot guarantee that T/e converges to the desired value since its gradient becomes zero. Nevertheless, 59.2% of the latent confounders in experiments of Section 4.5 were on ancestor-descendart pairs. ENCO detects a majority of those confounders, showing that ENCO still works such confounders despite not having guarantees. Further, show in Section C.3 that the confounder scores Ic(X,, Xÿ) indeed converge to one for detected confounders, and zero for all other edges.
B.4 CONVERGENCE FOR PARTIAL INTERVENTION SETS
In Section 4.4, we have experimentally shown that ENCO works on partial intervention sets as well. Here, we will discuss convergence guarantees in the situation when interventions are not provided on all variables.
We start with discussing the case where, for a graph with N variables, we are given samples of interventions on N - variables. In this case, we can rely on the previous convergence guarantees discussed in Appendix B.2 with minor modifications. Specifically, for the variable X; on which we do not have interventions, the orientation parameters O, are only updated by interventions other variables. Hence, for this variable X,, the following conditions need to hold instead of Theorem B:
For all possible sets of parents of X, excluding Xj, adding X, does not improve the log-likelihood estimate of X, under the intervention on Xj, or leaves it unchanged:
Ypa(X.) C X Eix, xx (oog((X/pp(X))) X,) looglog lopp((X]paa()))] <
(43)
For at least one parent set (X,),), which has a probability greater than zero to be sampled, this inequality is strictly smaller than zero.
32

---

Published as a conference paper at ICLR 2022
This condition ensures that 8. converges to the correct values, where X, is the parent ancestor of Xy. Thus, in conclusion, we can provide convergence guarantees
interventions are provided.
Next, we can consider the case of having N interventions. With the conditions above, can ensure that the all orientation parameters are learned, excluding 0,z where X; and X; are the variables for which we have not obtained interventions. In this case, we cannot give strict convergence guarantees for the edge X; <> X,, especially when X, and Xj have a direct causal relationship. If X; is the child of X,, we might obtain the edge X; + X; which violates the assumptions in the second and third step of the proof. Therefore, we cannot give guarantees of correctness for incoming/outgoing edges of X; and Xj. and might make incorrect predictions of edges between these two variables.
When taking the next step to having M interventions provided, ENCO can create more incorrect predictions. For the variables for which interventions are provided, we can use the same convergence guarantees (Theorem B. BB.4). since all conditions are independent across variables. For variables without interventions, we cannot rely on those. While we have observed that learning the missing 0's from other interventions give reasonable results, we see a degradation of performance the further the distance is between node and the closest intervened variable. As an example, suppose we have chain with 5 variables, i.c. X1 + X2 - X3 + X1 Xs. and we are provided with an intervention on X1 only. This allows to learn the orientation between X1 and X2. The orientation between X2 and X3 is often learned correctly as well because adding the edge X2 -> X3 instead of X3 - X2 gives greater decrease in overall log-likelihood, since part of the information from X3 to predict X₂ is already included in 1 However, the further we go away from X1, the less information is shared between the child and the intervened variable. Moreover, the likelihood of mistake occurring due to limited data further increases. This is why the orientation of the edge X4 + Xs is not always learned correctly, which can also cause false positive edges.
Many scenarios of predicting false positive edges can, in theory, be solved by providing an undirected skeleton of the graph, for example, obtained from observational data. Still, one of the cornerstones of ENCO is that it does not assume faithfulness. Without faithfulness or any other assumption on the functional form of the causal mechanisms, the correct undirected graph cannot be recovered by any method. One of the future directions will be to include faithfulness in ENCO to solve the scenarios mentioned above, although this would imply that we might not be able to recover edges of deterministic variables anymore.
B.5 EXAMPLE FOR NON-FAITHFUU GRAPIS
Below we give an example of a distribution which is not faithful with respect to its graph structure,
but can yet be found by ENCO. Suppose we have a chain of three variables, X; - X2 + X3. For simplicity, we assume here that all the variables are binary, but the argument can similarly hold for any categorical data. The distribution p(X1) is an arbitrary function with 0 < p(X, = 1) < 1, and the other two conditionals
deterministic functions: XX) 6|X₂ = X1), p(Xx/X22) 5 5[X3- = X2). This joint distribution is not faithful to the graph, since X₃ is independent of X₂ given X1, which is not implied by the graph. We will now focus our discussion
the edge X; X3 to show that, despite the independence, the proposed method ENCO can identify the true parent set of Xa. The first step of ENCO is fit the neural networks to the observational distributions, which include p(X3), pl X Xx)X,),, ), p( X2), p(Xx| X2)X2) X2).
Now, the update of the edge parameter hi under interventions on X2 can be summarized as follows, where we marginalize out over graph samples:
3c
= - (fru)) Ex p(X22 + Xx) 1----lgge log + log p(Xx
0713
Gradients for graph sumples with edge X X
(44)
p(X2 f Nz)
[- logg xX X1) + log p( Xx))
Gradients for graph samples without edge X2
where the samples X are sampled from the graph under interventions on X2. Intuitively, the gradient points towards increasing the probability of the edge X₁ + X3 if adding X, to the conditionals improves the log-likeli likellood estimate of X3 in both graph samples, ie. where we have the edge X2 X3 or not. Thus, to show that the gradient points towards decreasing the probability of the
33

---

Published as a conference paper at ICLR 2022
edge X1 - X3, need to show that both of the above log- oikeleiidhood differences greater than zero (note that positive gradients lead to decrease since we minimize the objective).
In the first difference, Ex I- log p( X1, X2) + ((X3 (2)), we know that p(X3 X2) is the optimal estimator, since the ground truth data is generated via this conditional. Thus, independent of what function the network has learned for p(X3 X2), the difference above can be only greater or equals to zero. Note that for (xx3 (2), there are value combinations that have never been observed, i.e. Xi X2, such that in practice we can't I guarantee specific distribution to be learned for such conditionals. Still, this does not constitute a problem for finding the graph as shown above.
The second difference, Ex | p X1) + log p( (x))]. can be similarly reasoned. Since X1 is independent of X3 under interventions on X2. p(X33 [x1) cannot lead to a better estimator than p(Xs). even when trained on the interventional data. However, since p(Xg|_X (x was trained on observational data, where there exist strong correlation between X1 and X3. this estimator must be strictly worse than p (X)). Hence, this difference will be strictly positive.
To summarize, under interventions on X2, the edge X1 + X3 will be trained towards decreasing its probability. Further, under interventions on X1, the effect of X on X3 can be fully expressed by conditioning X2, making this gradient going to zero when the edge probability X₂ - X3 goes towards one. For the edge X2 ~ X3 itself, the same reasoning as here can be followed such that independent of whether the edge X1 + X3 is included in the graph not, conditioning X3 X2 can only lead to an improvement in its estimator. Therefore, ENCO is able to find the correct graph despite it not being faithful.
34

---

Published as a conference paper at ICLR 2022
bidiag
chain
collider
full
jungle
random
Figure 10: Visualization of the common graph structures for graphs with 8 nodes. The graphs used in the experiments had 25 nodes. Note that the graph random more densely connected for larger graphs, as the number of possible edges scales quadractically with the number of nodes.
C EXPERIMENTAL DETAILS
The following section gives overview of the hyperparameters used across experiments. Addition- ally, we discuss details of the graph generation and the learning process of different algorithms.
C.1 COMMON GRAPH STRUCTURE EXPERIMENTS
C.1.1 DATASETS
Graph generation The six common graph structures we have used for the experiments in Table visualized in Figure 10. In the graph bidiag, variable X; has X;-2 and X;-I as parents, and consequently X,+1 and X1+2 as children. Hence, this graph represents a chain with 2-hop connections. The graph chain is bidiag ag with single hop, meaning that X; is the parent of X,-1 but not X₁-2. In the graph collider the variable X has all other variables, X_,, as parents. In the graph full. the parents of variable X; all previous variables: pa( (X)) (X1, X2.... X,- Hence, it the densest connected graph possible. The graph jungle 1 represents a binary tree where node is also connected to its parent's parent. Finally, the graph random follows randomly sampled adjacency matrix. For every possible pair of variables X,, we sample an edge with a likelihood of 0.3. To determine the orientation of the edges, we assume the causal ordering of X; X In other words, if we have an edge between X; and X,. it is oriented X; X, is else X, + X;-
Conditional distributions In the generated graphs, we categorical variables that each have 10 categories. To model the ground-truth conditional distributions, we use randomly initialized neural networks. Specifically, we use MLPs of two layers where the categorical inputs are represented by embedding vectors. For a variable X, with M parents, we stack the M embedding vectors to form the input to the following MLPs. Each embedding has dimensionality of 4, hence the input size to the first linear layer is 4M. The hidden size of the layers is 48, and we use LeakyReLU activation function in between the two linear layers. Finally, softmax activation function is used on
35

---

Published as a conference paper at ICLR 2022
Table 4: Hyperparameter overview for the simulated graphs dataset experiments presented in Table 1.
Hyperparameters
SDI
ENCO
Sparsity regularizer Xxparse
{0.01, 0.02, 0.05, 0.1, 0.2)
{0.002, 0.004, 0.01)
DAG regularizer
(0.2, 0.5, 1.0, 2.0, 5.0)
Distribution model
2 layers, hidden size 64, Leaky 0.1)
Batch size
128
Learning rate model
(2e-3, 5e-3, 2e-2, 5e-2)
Weight decay model
(le-5, 1e-4, le-5)
Distribution fitting iterations F
1000
Graph fitting iterations G
100
Graph samples
100
Epochs
50
30 (5e-3, 2e-2, 5e-2)
Learning rate 7
(5e-3, 2e-2, 5e-2)
Learning rate
(5e-3, 2e-2, 5e-2, le-1)
the output to obtain a distribution over the 10 categories. The MLP and hyperparameters have been chosen based on the design networks used in ENCO, SDI (Ke et al., 2019) and DCDI (Brouillard et al., 2020). For the initialization of the networks, we follow Ke et al. (2019) and use the orthogonal initialize with gain 2.5. The biases are thereby initialized uniformly between -0.5 and 0.5. This way, we obtain non-trivial, random distributions. Experiments with different synthetic distributions are provided in Appendix D.7.
C.1.2 METHODS AND HYPERPARAMETEN
Baseline implementation We used existing implementations to run the baselines GIES (Hauser Bühlmann, 2012), IGSP (Wang et al., 2017), GES (Chickering, 2002) and DCDI (Brouillard et al., 2020). For GIES, we used the implementation from the R package poalg?. To run categorical data, used the GaussLOpenInt Score score function. For IGSP, used the implementation of the python package causal dag' As IGSP uses conditional independence tests in its score function, cast the categorical data into continuous space first and experiment with different kernel-based independence tests. Due to its long runtime for large dataset sizes, we limit the interventional and observational data set size 25k. Larger dataset sizes did not show any significant improvements. For details on the observational GES experiments, see Section D.6. Finally, we have used the original python code for DCDI published by the authors*. We have added the same neural networks used by ENCO into the framework to perform structure learning on categorical data. Bugs in the original code were corrected the best of our knowledge. Since SDI (Ke et al., 2019) has a similar learning structure as ENCO, we have implemented it in the same code base as ENCO. This allows us to compare the learning algorithms under exact same perquisites. Further, all methods with neural networks used the deep learning framework PyTorch (Paszke et al., 2019) which ensures fair run time comparison across methods.
To ensure a fair comparison, we performed hyperparameter search for all methods. The hyperparameter search was performed on hold-out set of graphs containing two of each graph structure.
Hyperparameters
GIES We performed hyperparameter search over the regularizer values A
(0.01, 0.02, 0.05, 011, 0.2, 0.5, 1,2,5, 5, 10, 20, 50, 100, 200). The value obtaining the best results in terms of structural hamming distance (SHD) was À = 20. The average run time of GIES was 2mins per graph.
hhtrps://rraa..-ppojjett httppsr//ggthbb..co////ee/eaaaaa
htttp:://githbb..omm//aaaaaeeee
36

---

Published as a conference paper at ICLR 2022
bidiag
chain
collider
1.0
1.0
1.00
Metrics
0.98
recall
0.8
0.8
0.96
precision
0.6
0.94
0.6
0.92
0.4
10
20
30
10
20
30
10
20
30
Epoch
Epoch
Epoch
full
jungle
random
1.000
1.0
1.00
0.985
0.9
0.95
0.990
0.8
0.90
0.985
0.7
0.980
0.85
Epoch
Epoch
Epoch
Figure 11: Learning curves of ENCO in terms of recall and precision on edge predictions for synthetic graph structures. The orientations for the ground-truth edges not plotted they have usually been correctly learned after the first epoch except for the graph full 1. Overall, we see that the edge recall starts very high for all graphs, and the precision catches up the epochs. This is in line the steps in the convergence proof in Section B.
IGSP We experimented with two different conditional independence tests, kci and hsic, c, and different cutoff values a (1e-5, 1e-4, 1e-3, 1e-2). The best hyperparameter setting was kci with a = 1e-3. The average run time of IGSP was 13mins.
SDI We focused the hyperparameter search for SDI on its two regularizers, Asparse and ÂDAG. as well as its learning rate for y. The other hyperparameters with respect to the neural networks were kept the same as ENCO for fair comparison. We show all details of the hyperparameter search in Table 4. The best combination of regularizers found was Asparse = 0.02 and ÀDAG 0.5. Lower values of Asparse lead to more false positives, especially in sparse graphs, while a lower value of ADAG caused many two-variable loops. Compared to the reported hyperparameter by Ke al. (2019) (ADAG 0.5, Aquase 0.1), we found lower sparsity regularizer to work better. This is likely because of testing SDI on larger graphs. In contrast ENCO, SDI needed a lower learning rate for ? due its higher variance gradient estimators. To compensate for it, we ran it for 50 instead of 30 epochs. In general, SDI achieved lower scores than in the original experiments by Ke al. (2019) which was because of the larger graph size and smaller dataset size. The average run time of SDI was 4mins per graph.
DCDI
The most crucial hyperparameter in DCDI was the initialization of the constraint factor Пg. We experimented with range of μu f1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5) and found =1-9 to work best. This close to the reported value of 1e-8 by Brouillard et al. (2020). Higher values lead to empty graphs, while lower values slowed down the optimization. Ad- ditionally, we search over the regularizer hyperparameter À E (1e-3, 1e-2, 1e-1, 1.0, 10.0) where found À 0.1, which the same value used by Brouillard et al. (2020). We stop the search after the Lagrangian constraint is below 1c-7, following Brouillard et al. (2020), 50k iterations have been used which was sufficient to converge on all graphs. We have experimented with using weight decay to prevent overfitting, but did not find it to provide any gain in terms of structure learning performance. The average run time of DCDI was 16 minutes.
Ho = 1e-9
ENCO
We outline the hyperparameters of ENCO in Table 4. As discussed in Section 3.4, the most crucial hyperparameter in ENCO the sparsity regularizer Xyparee- The larger it is, the faster ENCO converges, but the same time might miss edges the ground-truth graph.
37

---

Published as a conference paper at ICLR 2022
Table 5: Extension of Table with the metric structural intervention distance (SID) (lower is better), averaged over 25 graphs each. The conclusion is the same as for SHD, namely that ENCO outperforms all baselines, while the acyclic heuristic has an even greater impact.
Gruph type
bidiag
chain
collider
full
jungle
random
GIES
160.0 (=60.1)
221.2 3)
(3.3
527. 9 (1:557)
441.4
(133.7)))
IGSP
423.3 (#-48.2)
240.1 (=78.8)
(+:1-4)
554.8(+226.))
394.8 (+73.5)
(21.88))
SDI
243.7 (=:55.1)
70.2 (4-46.8)
(4:00)
537.1 (4:29.2)
180.0 (4:56.4)
(4:6777
DCDI
369.3 (#:47.5)
233.4 (#24.8)
10.9 (#3.6)
(83.6()
339.4 (4+59.1)
158.6 (4-69.1)
ENCO (ours)
28.2 (+:8.8)
(1:1.)))
(=133))
63.2 (4:22.6
16.3 (49.8)
77.6 2)
ENCO-acyclic (ours)
(++0.)
0.0 (+900)
7.4 (:183))))
(=..9
(=8.1)
5.3 (1:1.88
Lower values allow the detection of more edges for the price of longer training times. We have found that for the graph structures given, only the graph full was sensitive to the value of Asaare where Asparse = 0.002 and Asparse = 0.004 performed almost equally well. In comparison to SDL, ENCO can make use of larger learning rates due to lower variance gradient estimators. Especially for 0, we have noticed that high learning rates are beneficial. This is in line with our theoretical guarantees which require the orientation parameters to converge first. We use the Adam optimizer for y and O with the 3-hyperparameters (0.9, 0.9) and (0.9, 0.999) respectively. A lower 32 hyperparameter for y allows the second momentum to adapt faster to change of gradient scale which happens for initial false positive predictions. The average run time of ENCO was 2mins per graph. The algorithm could be sped up even more by reducing the number of graph samples K and model fitting iterations. However, for graphs of larger than 100 nodes, = 100 and longer model fitting times showed to be beneficial. The learning curves in terms of recall and precision shown in Figure 10.
Enforcing acyclicity ENCO is guaranteed to converge to acyelic graphs in the data limit; arguably, an assumption that does not always hold. In the presence of cycles, which can occur especially when low data is available, a simple heuristic is to keep the graph, which maximizes the orientation probabilities. Specifically,
aim to find the order O SN. where Sx represents the set of all permutations from 1 to the number of variables N, for which we maximize the following objective:
INN
N
Ô = arg max
II =1
II
(00..,,)
(45)
O
For small cycles it is easy to do this exhaustively by checking all permutations. For larger cycles, apply a simple greedy search that works just as well. Once the order Ô has been found, we remove all edges X; + X, where comes after j in Ó. This guarantees to result in an acyelic graph.
The intuition behind this heuristic is the following. Cycles are often caused by a single orientation pair being incorrect due to noise in the interventional data. For example, in a chain X; + X2 + X3 + X4 + Xs, it can happen that the orientation parameter 014 is incorrectly learned as orientation of the edge between X2 and X, as X,, + X₁ if the interventional data on X, does not show the independence of X, and X4. However, most other orientation parameters, e.g. 012, 013, 624, Bz4, etc., have been likely learned correctly. Thus, it is easy to spot that 014 an outlier, and this what the simple heuristic above implements.
Results Besides the structural hamming distance, a common alternative metric is structural inter- vention distance (SID) (Peters& & Bühlmann, 2015). In contrast to SHD, SID quantifies the closeness between two DAGs in terms of their corresponding causal inference statements. Hence, it is suited for comparing causal discovery methods. The results of the experiments the synthetic graphs in terms of SID are shown in Table 5, and show a similar trend as before, namely that ENCO is outperforming all baselines.
C.2 SCALABILITY EXPERIMENTS
Graph generation For generating the graphs, we use the same strategy as for the graph x random in the previous experiments. The probability of sampling edge is set to 8/N, meaning that on
38

---

Published as a conference paper at ICLR 2022
random with nodes
1.0
0.9
0.8
0.7
Metrics
0.8
recall
precision
5
10
15
20
25
30
Epoch
(a) Example graph 100 nodes
(b) Learning curve 1, 000 nodes
Figure 12: (a) Example graph of 100 variables (best viewed electronically). Every node has on average edges and maximum of 10 parents. (b) Plotting recall and precision of the edge predictions for the training on graphs with 1, 000 nodes. The small standard deviation across graphs shows that ENCO can reliably recover large graphs.
average, every node has in- and outgoing edges. We limit the number of parents to 10 per node since, otherwise, cannot guarantee that the randomly sampled distributions take all parents faithfully into account. This is also in line with the real-world inspired graphs of the Bnl Bnre repository, which have maximum of parents. To give intuition on the complexity of such graphs, we show an example graph of 100 nodes in Figure 12a. Accordingly to the number of variables, we have increased the data size to 4096 samples per intervention, and I 100k observational samples. We did not apply the order heuristic on the predictions, since ENCO was able to recover acyclic graphs by itself with the given data.
y-frrezing stage For ENCO, one challenge of large graphs is that the orientation parameters O are updated very sparsely. The gradients for 0.j require data from an intervention on one of its adjacent nodes X; Xj. which we evaluate less frequently with increasing N we iterate over interventions on all N nodes. Hence, we require more iterations/epochs just for training the orientation parameters while wasting a lot of computational resources. To accelerate training of large graphs, we freeze in every second graph fitting stage. Updating only allows us to use the same graph sample C_ for both -x,-x, (X,) and Lx,4x, (X,) since the log-likelihood estimate of X; only needs to be evaluated for 8.y. With this gradient estimator, we experience that as little as 4 graph samples are gradient
sufficient to obtain reasonable gradient variance. Hence, it is possible perform updates of 0 in the same computation time. Note that this is estimator not efficient when training y as require different C-ij samples for every i. In experiments, step with this pure f-training stage. We want to emphasize that this approach can also be used for small graphs obtaining similar results in Table 1. However, is less needed because the orientation parameters are more frequently updated in the first place. Such an approach is not possible for the baselines, SDI and DCDI, because they do not model the orientation separate variable.
alternate the standard graph fitting
Hyperparameters To experiment with large graphs, we mostly keep to the same hyperparameters as reported in Section C.I. However, all methods showed gain by small hyperparameter search. For SDI and ENCO, we increase the number of distribution fitting iterations as the neural networks need to model larger set of possible parents. We also increase the learning rate of r to 2e-2. However, while SDI reaches better performance with the increased learning rate at epoch 30, it showed to perform worse when training for longer. This indicates that high learning rates can lead to local minima in SDI. Additionally, convergence speed for ENCO while SDI did not improve with a higher sparsity regularizer. Table shows hyperparameter overview of ENCO large-scale graphs, and Figure 12b the learning curve on graphs of 1, 000 nodes.
noticed that slightly higher sparsity regularizer improved
For DCDI, noticed that the hyperparameters around the Lagrangian constraint needed to be carefully fine-tuned. The Lagrangian constraint can reach values larger than possible to represent with double, and starts with 8e216 for graphs of 0000 nodes. Following Brouillard et al. (2020),
39

---

Published as a conference paper at ICLR 2022
normalize the constraint by the value after initialization, which gives more reasonable value to start learning. We performed another hyperparameter search on Ho, noticed however that it did not have a major impact. In the run time of ENCO, DCDI just starts to increase the weighting factor of the augmented Lagrangian while the DAG constraint is lower than 1e-10 for the smallest graph. The best value found was f40 1e-7.
Table 6: Hyperparameter overview of ENCO for the scalability experiments presented in Table 7. Numbers in the center represent that we use the same value for all graphs.
Hyperparameters
100 nodes
200 nodes
400 nodes
1000 nodes
Distribution model
2 layers, hidden size 64, LeakyReLU(a = 0.1)
Batch size - model
128
128
64 128
Learning rate- model
5c-3
Distribution fitting iterations F
2000
2000
4000
4000
Graph fitting iterations G
100
Graph samples
100
Learning rate- Y
2e-2
Learning rate-
lc-1
O-training iterations O-training graph samples Sparsity regularizer Asparse
1000
1000
2000
2000
2+2
(0.002, 0.004, 0.01)
Number of epochs (max.)
30
Results For clarity, we report the results of all methods below. The exact values might not be easily readable in Figure due to large differences in performance.
Table 7: Results for graphs between 100 and 1000 nodes. We report the average and standard deviation of the structural hamming distance (SHD) over 10 randomly sampled graphs. Thee maximum runtime of ENCO was measured on an NVIDIA RTX3090. Baselines were executed on the same hardware.
Graph size
100 15mins
200
400
1000
Max Runtime
45mins
2.5hrs
9hrs
DCDI (Brouillard et al., 2020)
583.1 (+21.8)
1399.0 (+-67.5
4761.2
(4:303.4)
OOM
SDI (Ke et al., 2019)
35.7 (12.9)
100.9 (17.6)
356.3 (111.5
1314.4 (1158.5)
ENCO (Ours)
0.0 (+0.0)
0.0 (++0.)
0.0 (+:000)
0.2 (+:0.4))
C.3 LATENT CONFOUNDER EXPERIMENTS
Graph generation The graphs used for testing the latent confounding strategy are based on the random graphs from Section 4.2. We use graphs of 25 nodes, and add 5 extra nodes that represent latent confounders. Each latent confounder X; connected to two randomly sampled nodes X,, X, that do not have a direct connection. However, X; and X, can be an ancestor-descendamt pair and have any other (shared) parent set (see Figure 13a). In the adjacency matrix, we add the edges X₁ + X; and X: Xj. and perform the data generation as for the previous graphs. After data generation, we remove the latent confounders from both observational and interventional data. The task is to learn the graph structure of the remaining 25 observable variables, as well as detecting whether there exists a latent confounder between any pair of variables. We use the same setup in terms of dataset size as before for the observational samples, namely 5k, but increased the samples per intervention to 512. Little interventional data showed to cause high variance the interventional gradients, N2 (7) which is why more false positives occured. The results in for the limited data with 200 interventions, and results in the data limit, í.e. for 10k interventional samples and 100k observational samples, shown in Table 8.
40

---

Published as a conference paper at ICLR 2022
Table Results of ENCO on detecting latent confounders averaged over 25 graphs with 25 nodes in the data limit (10k samples per intervention, 100k observational samples) and limited data (200 samples per intervention, 5k observational samples). In the data limit, only false negative predictions of latent confounders occured which did not affect other edge predictions. With little interventional data, more false positives occur reducing the precision.
Metrics
ENCO
ENCO
ENCO (200 interv./Sk observ.)
(10k interv./100k observ.)
(512 interv./5k observ.)
SHD Confounder recall
0.0 (=0.0)
1.24 (#1.33)
4.12 (=1.86) 92.0% (=11.5%)
96.8% (49.5%)
93.6% (# 13.8%)
Confounder precision
100.0% (+0.0%)
96.5% (#7.1%)
83.8% ( 16.4%)
Confounder score Ic(X, X)
0.6
0.4
0.4
True confounders
Xi
0.2
Max. other pairs
0.0
X
Epoch
(a) Latent confounder example
(b) Threshold sensitivity
Figure 13: Left: Example of latent confounder scenario, where X₁ is not observed and introduces dependency between X; and X, on observational data. The dots on the left and right represent even- tual (observed) parents of X, and X,. Right: Plotting the average score (e(X,, Xj) for confounders X; X₁ X; in the true causal graph (orange) and maximum score of any other node pair (blue). The plot shows the detection latent confounders in ENCO not sensitive to the specific value of т.
Baselines None our previous continuous optimization baselines, i.e., SDI and DCDI, are able to deal with latent confounders. To the best of our knowledge, other methods that are able handle latent confounders commonly take assumptions that do not hold in our experimental setup. Further, most methods are able deal with latent confounders in the sense that they obtain the correct results despite latent confounding being present. However, in our case, we explicitly predict latent confounders which is different task by itself.
Hyperparameters To show that we can perform latent confounder detection without specific hyperparameters, we use the same hyperparameters structures (see Appendix C.3). record 7 and Tij separately, we use separate first and second order momentum parameters in the Adam optimizer. We plot in Figure 13b the latent confounder scores Ic(X,, X2) calculated based on Equation 8. We see that the score converges close to 1 for pairs with a latent confounder, and for all other, it converges to O. This verifies the score function discussed in Section 3.5, and also shows that the method is not sensitive to the threshold hyperparameter
for the experiment on the previous graph (0)
motivation of
We choose recorded for any other pair at early stages of training.
0.4 which slightly higher than the highest value
C.4 INTERVENTIONS ON FEWER VARIABLES
Datasets We perform the experiments of interventions on fewer variables on the same graphs and datasets used for the initial synthetic graphs (see Section C. To simulate having interventions fewer variables, we randomly sample a subset variables for which we include the interventional data, and remove for all others. The sampled variables are the same for both ENCO and DCDI, and
41

---

Published as a conference paper at ICLR 2022
differ across graphs. The dataset size is the same before, namely 200 samples per intervention and 5k observational datasets.
ENCO for partial intervention sets While the theoretical guarantees for convergence to an acyclic graph apply when interventions on all variables are possible, it is straightforward to extend the ENCO algorithm support partial interventions as well. Normally, in the graph fitting stage, sample one intervention at time, We can, thus, simply restrict the sampling only to the interventions that are possible (or provided in the dataset). In this case, we update the orientation parameters 6.z of only those edges that connect to intervened variable, either X; or Xj, as before. All other orientation parameters would remain unchanged throughout the training, since their gradients rely on interventions missing from the dataset. Instead, extend the gradient estimator in Equation 4 to not be exclusive to adjacent interventions, but include interventions on all variables. Specifically, for the orientation parameter Byy without any interventions on X; or Xj, we use the following gradient estimator:
0
Ż - d'(0)) 00
Exxg,x.C_1 ((0))) (n(ss) (X,,)x((X,) - Lx, (X,)
(46)
((yi)) E Exx, x,c- (x,
Ax xx (X:))
where we have an intervention on an arbitrary variable X with + i, + j. This still represents an unbiased gradient estimator since in the derivation of the estimator, excluded interventions other variables only to reduce noise.
ENCO has been designed under the assumption that interventional data is provided. When have interventional data on only a very small subset of variables, we might not optimally use the information that is provided by the observational data. To overcome this issue, we can run causal discovery method that solely work on observational data and return an undirected graph. This skeleton can be used prior, and prevents false positive edges between conditionally independent variables.
Hyperparameters We reuse the hyperparameters of the experiments on the synthetic graph except that we use a slighly smaller sparsity regularizer, Aspanse 0.002, and weight decay of 4e-5. For the orientation parameters without adjacent intervention, we use learning rate of 0.1 Irg which is le-3 for this experiment. For DCDI, we observed that a higher regularization term of À = 1.0 obtained best performance. All other hyperparameters are the same as in Section C.I.
Results For additional experimental results, see Section D.2.
C.5 REAL-WORLD INSPIRED EXPERIMENTS
Datasets We perform experiments on a collection of causal graphs from the Bayesian Network Repository (Bnleear)) (Scutari, 2010). The repository contains graphs inspired by real-world appli- cations that are used benchmarks in literature. We chose the graphs to reflect a variety of sizes and different challenges (rare events, deterministic variables, etc.). The chosen graphs are cancer (Korb & Nicholson, 2010), earthquake (Korb & Nicholson, 2010), asia (Lauritzen & Spiegelhalter, 1988), sachs (Sachs et al., 2005), child (Spiegelhalter & Cowell, 1992), alarm (Beinlich et al., 1989), diabetes (Andreassen et al., 1991), and pigs (Scutari, 2010). The graphs have been downloaded from the BnLearn website". For the small graphs, have used dataset size of 50k observational samples and 512 samples per intervention. This is a larger dataset size than for the synthetic graph because many edges in the real-world graphs have very small causal effects that cannot be recovered from limited data, and the goal of the experiment was to show that the convergence conditions also hold on real-world graphs. Hence, we need more observational and interventional samples. The results with a smaller dataset size, i.e. 5k observational and 200 interventional samples as before, are shown in Table For the large graphs, we follow the dataset size for the scalability experiments (see Section C.2).
We reuse most of the hyperparameters of the previous experiments. For all graphs less than 100 nodes, we use the hyperparameters of Appendix C.I, i.e. the synthetic graphs of
Hyperparameters
'https://www //t::///wwwwww//... www.bnnweaaa.. Teaeear.commm www.ea...ooooeeeeeeeeeoooeeeeeee tory/
42

---

Published as a conference paper at ICLR 2022
Table 9: Results on graphs from the BnLearn library measured in structural hamming distance (lower is better). Results are averaged over 5 seeds with standard deviations.
Dutaset
cuncer (5 nodes)
eurthquake (5 nodes)
sachs
child (20 nodes)
alarm (37 nodes)
diabetes (413 nodes)
pigs
(8 nodes)
(11 nodes)
(441 nodes)
SDI
3.0 (++000
0.4 (+:05)
4.0 (+000)
7.0 (+0.0)
11.2 (+0.4)
24.4 (+1177)
422.4 (#8.7)
18.0 (++1.))
DCDI
4.0 (+:00)
2.0 (+0.0)
5.0 (+:0.0
5.4 (4:2.)
8.4 (+0.7)
30.0 (1:4.2)
ENCO
0.0 (+000)
0.0 (#:0.0)
0.0 (#000)
0.0 (=0.0)
0.0 (1:0.0)
1.0 (#0.0)
2.0 (+0.0)
0.0 (#:0.0)
Table 10: Results on graphs from the BnLearn library measured in structural hamming distance (lower is better), using 5k observational and 200 interventional samples.
Dataset
cancer (5 nodes)
earthquake (5 nodes)
asia (8 nodes)
sachs
child
alarm (37 nodes)
(11 nodes)
(20 nodes)
SDI
3.0 (+0.0)
0.4 (+0.5)
4.6 (1055)
8.4 (1..5)
12.4 (=0.9)
26.6 (#:1.1)
DCDI
4.0 (10.0)
2.0 (10.0)
4.4 (=0.7)
7.2 (12.4)
9.8 (10.6)
31.4 (+0.7)
ENCO
1.2 (+0.4)
0.4 (+0.9)
1.4 (+0.5)
0.4 (+0.5)
0.8 (=1.1)
11.4 (+1.8)
25 nodes. For all graphs larger than 100 nodes, we the hyperparameters of Appendix C.2, i.e. the large-scale graphs. One exception is that we allow the fine-tuning of the regularizer parameter for both sets. For ENCO, we used a slightly smaller regularizer. Asparse 0.002, for the small graphs and a larger one, Asparse 0.02, for the large graphs. Due to the large amount of deterministic variables, ENCO tends to predict more false positives in the beginning before removing them one by one. For SDI, we also found a smaller regularizer, Xaparse = 0.01, to work best for the small graphs. However, line with the results of Ke et (2019), SDI was not able to detect all edges. Even lower regularizers showed to perform considerably worse the child dataset, while minor improvements were made on the small graphs. Hence, we settled for Asparse 0.01. In terms of run time, both methods used 100 epochs for the small graphs and 50 for the large graphs
Results The results including standard deviations can be found in Table 9. The low standard deviation for ENCO shows that the approach is stable across seeds, even for large graphs. SDI has zero standard deviation for a few graphs. In those cases, SDI converged to the same graph across seeds, but not necessarily the correct graph. We have also applied DCDI (Brouillard et al., 2020) to the real-world datasets and report the results in Table 9 and 10. DCDI performs relatively similar to SDI, making few more mistakes on the very small graphs (< 10 nodes) while being slightly better on sachs and ch d. Nonetheless, ENCO outperforms DCDI on all graphs. We do not report results of DCDI on the largest graphs, diabetes and pigs, because ran out memory for diabetes (larger number of max. categories per variable) and did not converge within the same time limitations as SDI and ENCO (see Section 4.3 for a comparison on scalability).
43

---

Published as a conference paper at ICLR 2022
D ADDITIONAL EXPERIMENTS
In this section, we show additional experiments performed as ablation studies of ENCO. First, we discuss further experiments We then discuss the effect of using our gradient estimators proposed in Section 3.4 compared to Bengio et al. (2020). Next, we show experiments on synthetic graphs with deterministic variables violating faithfulness, and experiments on continuous data with Normalizing Flows. Finally, we discuss experiments with different causal mechanism functions for generating synthetic, conditional categorical distributions besides neural networks.
D. EFFECT OF THE SAMPLE SIZE
The number of samples provided as observational and interventional data crucial for causal structure learning methods since the more data we have, the better we can estimate the underlying causal mechanisms. To gain further insights in the effect of the sample size on ENCO and the compared baselines, we repeat the experiments of Section 4.2 with different sample sizes.
Large sample size First, we use very large sample sizes to find the upper bound performance level that we can expect from each method. For this, we sample 100k observational samples per graph, and 10k samples per intervention. We observed that this is sufficient to model most conditional probabilities up to a negligible error. The results are shown in Table 11. We find that, in line with the theoretical guarantees, ENCO can reliably recover most graphs, only making 0.3 mistakes average on the full graph. Of the baselines, only DCDI is able to recover the collider graph without errors since its edges can be independently orientated. For all other graphs, DCDI converges to acyclic graphs, but incorrectly orients some edges and predicts false positive edges, while being times slower than ENCO on the same hardware. All other baselines show improved SHD scores than in Table as well, but are not able to match ENCO's performance. This shows that, even in the data limit, ENCO achieves notably better results than concurrent methods.
Next, we consider situations where data is very limited. Thereby, we consider two data sample axes: observational and interventional data.
Limited interventional data sample sizes We repeat the experiments of Table I for ENCO while limiting the sample size per intervention to 20, 50, and 100 (200 before). The observational dataset size of 5000 samples is thereby kept constant. We plot the performance for all graph structures in Figure 14. Overall, the decrease of performance with lower interventional sample size is consistent across graph structures. With only 20 samples per intervention, it becomes especially hard to reason about variables with many parents, since the variable's distribution is determined by many other parents as well. Yet, for four out of the six graphs, we obtain an SHD of less than with 100 interventional samples, and less than 6 when only 20 samples are available. In conclusion, ENCO works well with little interventional data most variables have a small parent set.
Limited observational data sample sizes Similarly as above, we repeat the experiments of Table 1 for ENCO but limit the observational sample size to 1000 and 2000 (5000 before) while keeping 200 samples per interventions. Observational data is important in ENCO for learning the conditional distributions. For variables with many parents, this becomes more difficult when fewer samples are available, because the input space grows exponentially with the number of parents. Thus, we
Table 11: Repeating experiments of Table with large sample sizes (10k samples per intervention, 100k observational samples). In line with the theoretical guarantees, ENCO can reliably recover five out of the six graph structures without errors.
Graph type
biding
chain
collider
full
jungle
random
GIES
47.4 (+5.2)
22.3 (+:3.5
13.3 (+3.0)
152.7 (+12.0)
53.9 (+8.9)
86.1 (#:12.0)
IGSP
33.0 (+4.2)
12.0 (+1.9)
23.4 (+2.2)
264.6 (+7.4)
38.6 (+5.77
76.3 (+7.7)
SDI
2.1 (+1.5)
0.8 (+0.9)
14.7 (+-10)
121.6 (1:8.4)
(+11))))
1.8 (+119)
DCDI
3.7 (+1.5)
4.0 (1:13)
(10.0)
2.8 (+2.1)
1.2 (11.5)
2.2 (+1.5)
ENCO (Ours)
0.0 (=000)
0.0 (+H000)
0.0 (++00)
0.3 (+-09)
0.0 (+00.)
0.0 (+0.0)
44

---

Published as a conference paper at ICLR 2022
bidiag
chain
collider
Method
10
SHD
ENCO-scyclic
100
200
100
200
100
200
Samples per intervennnnonn
Samples per intervention
Samples per intervention
jungle
10
SHD
0
20 100 200
20
100
200
20
100
200
Samples per intervention
Samples per intervention
Samples per intervention
Figure 14: Results of ENCO for different graph structures under limited interventional data sample size. Note the different scale of the y-axis for the six graphs. While the general trend is the same for all graphs, i.e. decreasing performance with fewer samples, the order heuristic reduce the SHI) error by a considerable margin for most graphs.
bidlag
chaln
collider
Method
5HD
ENCO
ENCO-acydic
0
0
1k
2k
Sk
1k
24
Sk
1k
2k
5k
Observational samples
Observational samples
Observational samples
full
jungle
random
5.0
SHD
SHD
2.5
0.0
2k
Observational samples
Observational samples
Observational samples
Figure 15: Results of ENCO for different graph structures under limited observational data sample
size. Note the different scale of the y-axis for the six graphs. The structure learning performance remains good for sparse graphs, and suffers for graphs with larger parent sets.
would expect the collider and full graph suffer the most from having less observational data, and this is indeed the case as shown by the results in Figure 15. The results of all other graphs are less affected, although interestingly, some become even better with less observational data. For the chain X1 X2 + ...XN, for instance, we observed that the learned conditional distributions picked up spurious correlations among variables, e.g., between X1 and X3 when modeling p( X1|X1. X²) which are, in the data limit, independent given X2. Since those correlations do not necessarily transfer to the interventional setting. it is casier to spot false positive edges, and we can obtain even better results than for the larger sample sizes. In conclusion, having sufficient observational data is crucial in ENCO for graphs with variables that have larger parent sets, while being less important for sparser graphs.
Limited interventional and observational data sample sizes Finally, we combine the smallest interventional and observational data sample sizes, and also include the results of the previously best baselines, SDI and DCDI, in Table 12. The results of ENCO show the combination the previous two effects: graphs consisting of variables with small parent sets can still be recovered well by ENCO, while errors increase for the collider and full graph. Similar trends are observed for SDI, while DCDI showed considerable decrease in performance for all graphs. In conclusion, ENCO still works well for graphs with smaller parent sets under a small observational and interventional data regime, and outperforms related baselines in this setting.
45

---

Published as a conference paper at ICLR 2022
Table 12: Repeating experiments of Table with very small sample sizes (20 samples per intervention. I observational samples). Despite the limited data, ENCO can recover graphs with small parent sets reasonably well, while the graphs collider and full suffer for all methods.
Graph type
biding
chain
collider
full
jungle
random
SDI DCDI
10.9 (=2.7)
(+1.5)
22.1 (#1.9)
211.0 (46.2)
10.4
22.7 (#7.4) 40.2 (18.4)
30.0 (14.2)
22.0 (11.5)
23.2 (11.3)
185.2 (14.5)
25.8 (12.7)
ENCO (Ours)
9.7 (43.6)
5.6 (11.7)
22.7 (#2.1)
132.6 (18.0)
8.1 (42.3)
18.4 (24.8)
ENCO-acyclic (Ours)
2.0 (+2.3)
2.7 (+2.
22.9 (#2.3)
88.4 (1666)
4.1 (+:2.0)
5.3 (+2.5)
bidiag
chain
collider
20
Method
SHD
20
SHD
20
SHD
DCDI
ENCO (Ours)
ENCO-acyclle (Ours)
4 6 810
15
20
25
4.6 810
15
20
25
6 810
15
25
Intervened variables
Intervened variables
Intervened variables
full
jungle
random
40
SHD
100
SHD
SHD
20
0
0
4 6 810 Intervened variables
15
20 25
15
20
25
4 810
15
Intervened variables
Intervened variables
Figure 16: Results of ENCO and DCDI for different graph structures under fewer interventions provided. Note the different scale of the y-axis for the six graphs. For four out of six graphs, ENCO outperforms DCDI even for few as 4 interventions, especially when enforcing acyclicity. The detailed numbers of the results are listed Table 13.
D.2 INTERVENTIONS ON FEWER VARIABLES
We have performed the experiments in Section 4.4 using fewer interventions for all six synthetic graph structures. The results are visualized in Figure 16, and presented in table-form in Table 13. From the experiments, we can see that ENCO with interventions on only 4 variables matches or outperforms DCDI with 10 interventions for out of the graph structures (bidiag, chain, full, jungle) when enforcing acyclicity. Especially SHD scores for the same number of interventions on variables, while DCDI incorrectly orientates many edges and predicts false positive edges between children. On the graph collider, we observed a high variance for settings with very few interventions. This is because when we intervene on the collider node itself, ENCO can deduce the orientation for all edges. Finally, we observe that enforcing acyclicity for ENCO reduces the error lot. This is because incorrectly orientated edges cause more false positives in this densely connected graph, which are removed with the cycles. We include a longer discussion of the limitations of ENCO on fewer interventions in Appendix B.4. Still, we conclude that, in practice, ENCO performs competitively DCDI, even when very few interventions are provided, and scales better to more interventions.
chain-like graphs such as jungle, ENCO achieves lower
the graph random,
D.3 ABLATION STUDY ON GRADIENT ESTIMATORS
To analyze the importance of the low-variance gradient estimators in ENCO, we repeat the experiments on synthetic graph structure from Section 4.2 where the gradient estimators of ENCO have been replaced by those from Bengio et al. (2020). The results are shown in Table 14. Overall, the scores very similar with minor differences for the graphs full and bidi ag. In comparisons to the learning curves in Figure 11, the curves with the gradient estimator Bengio et al. (2020) are more noisy, with recall and precision jumping up and down. While ENCO easily converged early to the correct graphs for all graph types, this model often required the full 30 iterations to reach the optimal recovery.
46

---

Published as a conference paper at ICLR 2022
Table 13: Detailed results of the experiments with fewer interventions. See Figure 16 for a visualiza- tion and discussion.
Graph type
bidiag
chain
collider
full
jungle
random
DCDI vars
25.8 (12.0)
(111.3)
(21.8)
143.8 (1.10.7)
38.5 (1:32)
26.3 (1.6.2
DCDI vars
24.3 (+222)
(+2.7)
142.2 (2.5
32.7 (+3.8)
23.1 (+7.2)
DCDI vars
23.5 (=1.4)
13.3 (#2.4)
122.0)
134.8 (+8.9)
28.8 (-7.4)
19.1 (#5.5)
DCDI vars
22.4(111.1)
(1.0
10.8
126.2 (14.2)
28.0
14.8 (16.3)
DCDI 15 vars
22.0 (11.9)
(42.1)
11.5 (43.7)
(17.1)
(-2.1)
12.2 (15.3)
DCDI 20 vars
20.8 (+:.4)
(+1.3)
(2209)
62.8
(+
10.2 (+3.0)
DCDI 25 vars
16.9 (=2.0)
(0.11)
(+3.6)
(+-08)))))
(7.11
7.7 (#3.2)
ENCO 4
(+5.6)
(+392)
(+6.6))
110.6
(0.61)
28.5 (+-33)
ENCO
(23.0)
(62.8)
(2.11))
97.5
23.4 (35.))
ENCO 8
213..8)
(1224)
(10.3))3)
86.3 (18.9)
19.4
(13..)
ENCO 10 vars
17.7
10..6
(.4
77.3
1
18.7
ENCO 15 vars
10.7 (+223)
(.1.5))))
3.8
52.8
9.3 (+2.3)
13.6 (+268)28)
ENCO 20 vars
(42.2)
(.1.5))))
(.19)
(2.4
(+1.1)
8,8 (.8.1
ENCO 25 vars
(+144)))
1.7 (-1.3)
(.16))
9.2 (#3.4)
(.7.))
(+19)
ENCO-acyclic 4 vars
(17.5))
10.6.)
(16.7))
(18.6))
19.8 (1160)
(
ENCO-acyclie 6 vars
15.8
(.8.44.)
10.1 (+6.1)
(++85)
15.6 (+7.))
6.3.2
ENCO-acyclic vars
(15.4)
5.4(1377)
(.6.)))3)
86.3
12.6 (14.6)
5.3
ENCO-acyclic 10 vars
(+5.8)
(-444.4)
(.69)
77.3 (+8.5)
9.3.9)
3.2 (+2.2
ENCO-acyclic 15 vars
(.3.7)
344228))
(.2.5)
(49.1
(.2.77
2.1.77
ENCO-acyclic 20 vars
0.3 (+:07)
(=1.22
(.4
(24.2)
(=1.5)
(
ENCO-acyclic 25 vars
0.0 (10.0)
(10.0)
1.6 (L1.))
5.3 (12.3)
(.1.1
0.2 (10.5)
Table 14: Extension of Table 11 with ablation study of using Bengio et al. (2020) gradients with ENCO.
Graph type
bidiag
chain
collider
full
jungle
random
ENCO (Ours)
Bengio et al. (2020) grads
0.1 (10.4)
(1000)
0.0 (10000
(+1.5)
0.0 (10.00
0.0 (+0.00
gradient estimator
0.0 (+-0.)
(+0.0)
0.0 (+0.0)
0.3 (-9.)
0.0 (+0.0)
0.0 (+00.0
The difference between the two gradient estimators becomes more apparent on large graphs. We repeated the experiments of Section 4.3 on the graphs with 100 nodes using the gradient estimator of Bengio et al. (2020). Within the 30 epochs, the model obtained an SHD of 15.4 on average over experiments, which is considerably higher than ENCO with the proposed gradient estimators (0.0). Still, this only half of the errors that SDI (Ke et al., 2019) with the same gradient estimator achieved. Hence, we can conclude that the proposed gradient estimators are beneficial for ENCO but not strictly necessary for small graphs. For large graphs, the low variance of the estimator becomes much more important.
D.4 DETERMINISTIC VARIABLES
In contrast to algorithms working on observational data, ENCO does not strictly require the faithful- ness assumption. Hence, we can apply ENCO to graphs with deterministic variables. Deterministic variables have a distribution that is defined by a one-to-one mapping of its parents' inputs to an output value. In other words, have the following distribution:
(((x)))(X.) = 1 [X, = ((p((X)))
(47)
where is arbitrary function. The difficulty of deterministic variables is that a variable X; can fully replaced by its parents pa(X,) ) in any conditional distribution. The only way we can identify deterministic variables is from interventional data, where intervention on X; breaks the dependency to its parents.
We have already tested ENCO on deterministic variables in the context of the real-world inspired graphs of Section 4.6. To have a more detailed analysis, we created synthetic graphs following the random graph setup with edge probability of 0.1 and average of two parents, and maximum of three parents. An example graph is shown in Figure 17. All variables except the leaf nodes have deterministic distributions, where the function / (X))) is randomly created by picking
47

---

Published as a conference paper at ICLR 2022
Table 15: Experiments on graphs with deterministic variables. The performance over 10 experiments is reported in terms of SHD with standard deviation in brackets. ENCO can recover most of the graphs with less than two errors.
Graph type
deterministic
SDI (Ke al., 2019)
20.6 (13.8)
ENCO (Ours)
1.4 (+1.3)
Figure 17: Example graph for the deterministic setup. We use the random setup with edge proba- bility 0.1 and limit number parents to All variables except the leaf nodes have deterministic distributions.
X1
X₃
X2
Figure 18: Example graph for cancelling paths X1 + X2 + X3 and X1 + X3. This can, for instance, occur when the conditional distribution of X₂ Dirac delta around X1: p(X)(X)) 6[X₂ = X1).
random output category for any pair of input values. We create 10 such graphs and the same hyperparameter setup as for the synthetic graphs, except that we increase the sparsity regularizer to Ayparse 0.02. We report the results in Table 15. line with the results on the real-world graphs, ENCO is able to recover most graphs with less than two errors. As a baseline, we apply SDI (Ke et al., 2019) and see significant higher eтror rate. The method predicts many false positives, including two-variable loops, but was also missing out some true positive edges. We conclude that ENCO also works well on deterministic graphs.
Cancellation of paths Besides deterministic nodes, a common example for faithfulness violation the cancellation of two paths. For instance, consider the causal graph with the three variables X1, X2, X3 shown in Figure 18, and the conditional distribution p( xx
8[X, X1). In this case, the two paths X, - X₂ + X3 and X + X3 cancel each other, i.e. X3 independent of X2 when conditioned on X,, and independent of X, when conditioned of the two graphs is necessary for describing the relations. Yet. ENCO can find the edge X₁ - X₃ by observing interventions
X₂. This implies that only one
X2, since in this case, X1 1 X2 and X1 L X2/ The remaining edges can be learned in the same manner. We also emperically verify this by running ENCO on the graph structure of Figure 18 with the three variables being binary. We set p X2 (1) 8[X2 Xi] for
48

---

Published as a conference paper at ICLR 2022
Table 16: Experiments on graph with continuous data from Brouillard al. (2020). The suffix "-G" denotes that the neural networks model Gaussian density, and "-DSF" a two-layer deep sigmoidal flow. ENCO outperforms all baselines in this scenario, verifying that ENCO also works on continuous data well.
Graph type
Linear
Nonlinear with additive noise
Nonlinear with non-additive noise
GIES
0.6 (#1.3)
9.1 (#))3)
(.6.1))1)
IGSP (best)
(=188)
5.3 (#3.0)
4.1 (+2.8)
DCDI-G
1.3 (+1.9)
5.2 (+7.5)
2.3 (+3.6)
DCDI-DSE
0.9 (+1.3)
4.2 (+5.6)
7.0 (+10.7)
GES + Orientating
2.5 (+))
12.5 (+9.2)
7.9 (+7.4)
ENCO-G (ours)
0.0 (100.)
(10.4)
0.1 (+0.3)
ENCO-DSF (ours)
0.3 (1007)
1.4 (1144
1.2 (11.5)
canceling the two paths, and the remaining distributions are randomly initialized. ENCO reconstructs the graph without errors, showing that it also works in practice.
D.5 CONTINUOUS DATA
We verify that ENCO works just as well with continuous data by performing the experiments on datasets from Brouillard et al. (2020) that contained interventions on all variables. In these datasets, the graphs consist of 10 variables with an average of one edge per variable, and deploy three different causal mechanisms: linear, nonlinear additive noise models, and nonlinear models with non-additive noise using neural networks. The datasets contain 909 observational samples and 909 samples per intervention. All results of GIES, IGSP, and DCDI have been taken from Brouillard et al. (2020) (Appendix C.7, Table 22-24). We follow the setup of Brouillard et al. (2020) and compare two different neural network setups. First, we use MLPs that model a Gaussian density by predicting a mean and variance variable (denoted by suffix G). The second setup uses normalizing flows, more specifically two-layer deep sigmoidal flow (Huang et al., 2018), which is flexible enough to model more complex distributions (denoted by suffix DSF). The rest of the experimental setup ENCO identical to the categorical case.
Results are shown in Table 16, and the observations are the same as with categorical data. ENCO outperforms all other methods in all settings, especially for the more complex distributions. The higher error rate for the DSF setup mostly due overfitting of the flow models. We conclude that ENCO works as accurately for both continuous and categorical data.
D.6 SKELETON LEARNING WITH OBSERVATIONAL BASELINE
To show the benefit of learning graph from observational and interventional data jointly, we compare ENCO to a simple observational baseline. This baseline first learns the skeleton the graph by applying greedy equivalence search (GES) (Chickering, 2002) the observational data. Then, for each interventional dataset, we apply GES as well and use those skeletons to orientate the edges of the original one. This can be done by checking for each undirected edge X Y whether X Y is in the skeleton of interventions on X or not. As a reference implementation of GES, we have used the one provided in the Causal Discovery Toolbox (Kalainathan al., 2020).
The results on continuous data are shown in Table 16. Since GES assumes linear mechanisms and gaussianity of the data, it is unsurprising that performs better on the linear Gaussian dataset than on the non-linear datasets. However, on all the three datasets, it constitutes the lowest performance compared to the other methods, including ENCO. This highlights the benefits of incorporating interventional data in the learning of the skeleton and graph structure. To gain further insights in comparison to the constraint-basee baseline, we repeat the experiments with smaller sample sizes. The original dataset has 909 samples for observational data and per intervention, and we sub-sample 500 and 100 of those respectively for simulating smaller dataset sizes. The results of those experiments can be found in Table 17. It is apparent that the results of GES on the linear dataset get considerably worse with fewer data samples being available, while ENCO-G able to reconstruct most graphs still without emrrs. Especially for the small dataset of 100 samples, noticed that the skeletons
49

---

Published as a conference paper at ICLR 2022
Table 17:
Experiments on graph with continuous data from Brouillard et al. (2020) with smaller sample sizes for both observational and interventional datasets (in brackets). ENCO shows to perform much better in smaller sample sizes than a skeleton+orientation method, underlining the benefit of learning the whole graph from observational and interventional jointl..
Graph type
Linear Gaussian
Nonlinear with additive noise
Nonlinear with non-additive noise
GES + Orientating (909)
2.5 (#:4.3)
12.5 (+992)
7.9 (#7.4)
ENCO-G (ours) (909)
0.0 (110.0
0.2 (1:0.)
0.1 (=0.3)
GES +
Orientating (500)
3.2 (4:4.3)
12.0 (#:922
9.2 (=7.6)
ENCO-G (ours) (500)
0.1 (+0.3)
0.5 (#0.7)
0.0 (=0.0)
GES Orientating (100)
6.1 (d:5.3
10.2 (16.7)
9.0 (16.5)
ENCO-G (ours) (100)
0.2 (=:0.4)
0.9 (=111)
1.3 (=0.8)
Table 18:
Comparing structure learning methods in terms of structural hamming distance (SHD) on common graph structures (lower is better), averaged over 25 graphs each. ENCO outperforms all baselines, and by enforcing acyclicity after training, can recover most graphs with minimal errors.
Graph type
biding
chain
collider
full
jungle
random
GIES
(+7.5)
(+7.3)
(+2.9)
216.5 (+1.22
33.1 (+-2.9)
(7:1.2)
IGSP
32.7 (=5.1)
(+2.3)
23.7 (#2.3)
35.9 (+5.2)
(++.0)
SDI
9.0.
(-2.0
16.1 (12.4)
(110.3
(.2.3)
10.8 (13,9)
DCDI
16.9 (+-20)
10.1 (+1.1)
10.9 (+3.6)
21.0 (+4.8)
(++..)
7.7 (+3.2)
GES + Orientating
14.8 (1:2.6)
(-0.)
20.8 (1:2.4)
282.8 (1:42)
14.7 (13.1)
60.1 (18899
ENCO (ours)
2.2 (=1.4)
1.7 (=13)
(+1.6)
9.2 (+3.4)
1.7 (+:13)
4.6 (+1.9)
ENCO-acyclic (ours)
0.0 (1000)
0.0 (+0.0)
1.6 (1116)
5.3 (12.3)
0.6 (41.1)
0.2 (1005)
Table 19:
Experiments with different data simulator, introducing independence among parents for each variable. Similar to the neural-based synthetic data, ENCO recovers most graphs with a minor error rate, outperforming other baselines.
Graph type
bidiag
chain
collider
full
jungle
random
GIES
30.7 (#3.1)
16.9
18.6 (#2.5)
(+-.0)
30.1 (#3.5)
110.0
IGSP
(14.2)
(12.8)
25.0 (114)
(13.5)
24.0
112.5 (14.9)
SDI
(42.7)
(+2.6)
2.8 (=16)
(+7.5)
14.8 (+3.2)
36.7 (+4.7)
DCDI
15.7 (1119)
(+1.3)
4.7 (=2.))
25.3 (1..2)
18.9 (13.7)
9.8 (#3.6)
ENCO (ours)
0.5 (1.0.6
(L0.6)
(-0.)
(11.)))))
1.2 (.11.)
(L1.8)
ENCO-scyclic (ours)
0.0 (++00)
0.0 (++.0)
0.8 (+0.8)
0.2 (=0.4)
0.3 (+006)
0.1 (+0.3)
found by GES on observational data already contained couple of mistakes. This shows that for small datasets, observational data alone might not be sufficient to find the correct skeleton while by jointly learning from observational and interventional data, we can yet find the graph up minor errors.
Further, we also apply GES on the categorical data with an additional hyperparameter search over the penalty discount. The results in Table 18 give similar conclusion as on the continuous data. While the baseline attains good scores for chains, it makes considerably more errors on all other graph structures than ENCO. This shows that ENCO is much more robust by jointly learning from observational and interventional data.
D.7 NON-NEURAL BASED DATA SIMULATORS
Using neural networks to generate the simulated data might give SDI, DCDI and ENCO an advantage in our comparisons since they rely on similar neural networks model the distribution. To verify that ENCO works for other simulated data similarly well, we run experiments on categorical data with other function forms for the causal mechanisms instead of neural networks. Since there is no straightforward way of defining 'linear' mechanisms for categorical data, we instead express conditional distribution product of independent, single conditionals:
Ixx p(X,
p(X:)pa(X:)) EE. Ixx p(X: ,;
(48)
50

---

Published as a conference paper at ICLR 2022
exp(ox,,
Σx; xxp R. (0, Hence, the effect of each variable in the parent set is independent of all others, similar to linear functions in the continuous case. The individual probability densities represent softmax distribution over variables sampled from Gaussian distribution.
with p(X,|X,)
We apply GIES, IGSP, DCDI, SDI and ENCO the same set of synthetic graph structures with these new causal mechanisms. Similar to the previous experiments, we provide 200 samples per intervention and 5k observational samples the algorithms, and repeat the experiments with 25 independently sampled graphs. The results Table 19 give the same conclusion as the experiments ceural-based causal mechanisms, namely that ENCO outperforms all baselines. Most methods experience a decrease in performance since the average causal effect of cach parent is lower than in the neural case where more complex interactions between parents can be modeled. Still, ENCO only shows minor decreases, having less than one mistake on average for every graph structure when applying the orientation heuristic.
51


Reviews:
[{"review_id": "oVivrZQ0sBz", "reviewer": "Reviewer_rJAC", "summary": "", "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": "I find the proposal very interesting, especially in light of strong empirical results. The paper is well written but in my opinion, many of the claims were not sufficiently justified.\n\n1. A starting question I have is why acyclicity is necessary to enforce if interventions are available on all variables. An approach such as ICP [Peters et al, 2016] applied to all variables separately should recover the underlying graph and does not require acyclicity. \n\n2. It is not clear that the authors ensure that the recovered graph is acyclic. It is mentioned before section 3.3 that \"... acyclicity can be enforced as a post-processing step\". What is this step? Where is it discussed?\n\n3. Interventions on all variables will obviously never be the case in practice, how do the authors think the proposed approach can be used in practice. What does the method converge to if interventions are not available on all variables? Can we characterize equivalence classes?\n\n4. Theorem 3.1 assumes access to conditional distributions if my reading is correct, is it plausible to expect that there will be no local optima with finite samples? With overparameterized models, this seems unlikely.\n\n5. After Theorem 3.1 the authors say \"we can guarantee to converge to the true causal graph if, for every edge Xi → Xj in the ground truth, the variables are not conditionally independent in the interventional data of Xi\". Isn't this a faithfulness condition, i.e. that a conditional independence implies a corresponding separation in the graph?\n\nMinor.\nScore-based methods are defined by the optimization of a score rather than by a search strategy in the space of DAGs which makes me uneasy about the new category of continuous-optimization methods.", "summary_of_the_paper": "The problem of score-based structure learning using observational data as well as a set of interventional datasets is increasingly relevant in domains where experiments are readily available. The authors propose a new parameterization of the DAG underlying the data and score showing, provided that experiments are available on all observed variables and the system be Markovian, that the true causal graph can be recovered. ", "summary_of_the_review": "Interesting proposal if some of the claims can be better justified and if we can bring the algorithm closer to practice. In its current formulation, both interventions on all variables and unlimited sample sizes to guarantee convergence are not plausible in my opinion.  ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "lOnSKXGo9GH", "reviewer": "Reviewer_qKQN", "summary": "", "questions": "", "limitations": "", "rating": 6, "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 is well-written. The claims seem sound and the experimental results are convincing. The assumptions are clearly described, and the limitations are also discussed. Also, the authors extend their method to detect hidden confounders, which is interesting and important.\n\nMy only question is about the distribution fitting stage. The authors build a neural network for each variable to model the conditional distributions of that variable, conditioning on any variable set. That is, given a variable $X_i$ and a conditioning set ${\\\\bf X}_I$, where the indexing set $I\\\\subset \\\\{1, 2, \\\\cdots, N\\\\}\\\\setminus\\\\{i\\\\}$, the authors define a mask ${\\\\bf M}\\\\in \\\\{0, 1\\\\}^{N}$ such that ${\\\\bf M}(j)=1$ if and only if $j=i$ or $j\\\\in I$, and model $P(X_i \\\\mid {\\\\bf X}_I)$ by the NN $f({\\\\bf M}\\\\circ {\\\\bf X}; {\\\\phi_i})$. I am wondering how this objective relates to the loss function given in Equation (1).\n", "summary_of_the_paper": "This paper proposes a new gradient-based method to learn causal DAGs from observational and interventional data. It models a probability for every possible directed edge between pairs of variables and searches for the graph which generalizes best from observational to interventional data.", "summary_of_the_review": "Overall, this paper is well-written. The claims are supported by rigorous theoretical analysis and extensive empirical studies. It seems that there is only a minor issue in Section 3.2. I think this is a good submission if the authors could clarify this point.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "f_kIvQ5diUT", "reviewer": "Reviewer_H7r5", "summary": "", "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": "Overall, the paper is well written and clearly structured. I appreciate Figure 1 which is a concise graphical overview of the method.\n\nThe method itself and how to train it is clearly explained in the sections 3.1-3.3. Regarding the convergence guarantees made in section 3.4, I am not an expert in this domain so I cannot confirm with high certainty that the proof is correct but the logical reasoning seems to be fine. I acknowledge that the authors have an example of how to check the conditions of Theorem 3.1 in appendix B.2.1. However, as far as I understand the paper, these steps can only be applied if the probabilistic model is known but this is not the case in an real world scenario. Especially condition 3 seems problematic as it sets the scale of the regularisation parameter $\\lambda_{sparse}$. For large graphs, it seems to be very complicated to estimate $\\lambda_{sparse}$ and, hence, there will be no guarantee that we recover the full causal graph.\n\nThe experiments are extensive and the authors proved that their method is clearly better than several baselines on different datasets. What seems strange to me is that they did only compare their method to SDI in section 4.6. They argue that it is in general the strongest baseline but in fact it is outperformed by DCDI on the sachs and child dataset according to Table 9. These results should be reported in the main paper as well, especially since they do do not need a lot of space. \nFor the paper to be truly outstanding, it would have been nice to see the method being applied to a real-word dataset. Although the true causal graph is not known in this setting, the plausibility of the result could be discussed in this case to evaluate the performance of the method.", "summary_of_the_paper": "The paper introduces a method, called Efficient Neural Causal Discovery (ENCO), to determine the structure of a causal graphical model of a set of random variables, which uses both observational and interventional data.\n\nThe training procedure of the method consists of two steps: a distribution of one random variable given the others is fitted through a neural network and by performing interventions the graph is fitted using the approximate distributions. By alternating the two steps the causal graph can be recovered. The authors state the conditions when this procedure converges and prove the respective convergence theorem. Furthermore, they show how latent confounders can be detected within their framework.\n\nThe authors apply their method to several synthetic datasets, among them real-world inspired data form the Bayesian Network Repository (BnLearn), to show that they can reconstruct the causal graph for various numbers of random variables as well as find latent confounders if they are present. Based on various performance metrics they outperform competing methods, such as Structural Discovery from Interventions (SDI).", "summary_of_the_review": "Overall, the paper is a significant contribution to the area of learning causal graphs. Besides several minor points, the article is well written and should be accepted.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "acE50LFI09_", "reviewer": "Reviewer_b5p5", "summary": "", "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": "\"Unfortunately, the solution space of DAGs grows super-exponentially with the number of variables. Current methods are thus typically limited to a few dozens of variables\"\nI am not sure if this reasoning is accurate since exhaustive search is not used.\n\n\"Yet, we show that under certain assumptions like being able to intervene on all variables, the proposed optimization is guaranteed to converge to the correct, acyclic graph\"\nThis is not an insignificant assumption so it should be front and center in the introduction instead of being mentioned in passing. \n\n\"Thirdly, under mild conditions including interventions on all variables, we show that ENCO is guaranteed to converge to the correct causal graph\"\nSimilarly, it is hard to justify this claim that this is a mild condition. I think it would be much more transparent to call the paper causal discovery from observational and interventional data - which describes the setting - rather than trying to present this as just another assumption.\n\n\"A DAG G and a joint distribution P are faithful to ea􏰍ch other\"\nDistribution is said to be faithful to a graph, not the other way around.\n\nImportant recent works that learn from obs+int data are missing from literature review under constraint-based methods:\n- Mooij et al. \"Joint causal inference from multiple contexts\"\n- Kocaoglu et al. \"Characterization and Learning of Causal Graphs with Latent Variables from Soft Interventions\" \n- Jaber et al. \"Causal Discovery from Soft Interventions with Unknown Targets: Characterization and Learning\"\n\n\"Since the search space of DAGs is super-exponential in the number of nodes, many score-based approaches rely on a greedy search\"\nAs far as I am aware of, there are consistency guarantees that these methods will output graphs in the true equivalence class. So it would be more accurate to mention this. \n\n\"We do not assume faithfulness\"\nThis is a great addition.\n\n\"perfect\" intervention is also called a \"hard\" intervention. \n\nIt would be nice to list all the assumptions in 3.1 under a numbered list. Specifically, authors assume that we have n interventional datasets for n variables where in each intervention a different node is intervened on (intervention size 1).\n\nThe discussion under \"Graph fitting\" on page 4 mentions that the method requires intervention on a node that is adjacent to the edge (although others might eventually be sufficient as well). However, if such an abundance of interventional data is available, one can simply learn each edge separately by contrasting two distributions and hypothesis testing without the need for a continuous optimization framework. Say X-Y and p_1 = p(y|x), p_2=p(y|do(x)) is available. Then we can check if p_1 and p_2 are the same or different to infer this edge. So essentially I am not sure if there is a need to leverage the power of continuous optimization or neural networks in this setting. More on this later.\n\nThe conditions of Theorem 3.1. seems to have been developed by essentially reverse-engineering the proof. In that sense, they are far from being intuitive about the data-generating conditions or system settings. Furthermore, due to this, it is not possible to assess whether these assumptions are realistic. This renders the theorem not really useful in my opinion and the authors can simply remove this unless the assumptions can be converted to intuitive ones descriptive of the setting irrespective of the proof. \n\nThe authors propose an interesting method to handle latent confounders. They have to assume that latent confounders cannot exist between adjacent variables. Even then though I am not sure if there is any guarantee because authors say \"which is maximized by latent confounders\" but I don't see a proof of this claim. It would help to clarify the exact claim here.\n\nIn the experiments, I am not sure why authors only compare with DCDI but not with the other mentioned methods.\n\nCould you comment on whether the improvements seen in Section 4.3 are due to the new edge-wise detection framework (not having to enforce acyclicity constraint) or due to the better-variance estimator described in Section 3.3? \n\nI really would like to see the performance of a very simple baseline: \n- Learn the skeleton somehow (e.g. w/ GES)\n- Orient each edge X-Y by learning the skeleton on the interventional distribution, say under an intervention on X. If X-Y still then X->Y else X<-Y. This is the basic idea that has been used in several constrained-based methods in the past when hard-interventional data is available. If only soft-interventional data is available one can check p_obs(y|x)==p_x(y|x) and orient X->Y if so and X<-Y otherwise. Of course, these require some sort of faithfulness but in settings with faithfulness, I think it is important to have these comparisons to clarify the contribution of the method. Perhaps your method is more sample-efficient and it would be great to demonstrate this empirically. \n\nFinally, I am not absolutely sure if your method can avoid faithfulness. Suppose X->Y is the true graph but we have X\\indep Y. Your (or any method really) will remove the edge between the two but that won't give the true graph. Any comments on this as well will be appreciated. \n\nAfter the rebuttal:   \n--------------------  \nI would like to thank the authors for the additional experiments they conducted and the further explanations about their assumptions. Accordingly, I will increase my score. However, due to the shortcomings mentioned in my original review and the theoretical assumptions still being non-intuitive, I will not be able to recommend a strong acceptance. In case the paper goes through, please implement all the discussed changes, and thank you again for keeping an open mind and enabling a productive rebuttal period!", "summary_of_the_paper": "The authors propose a differentiable causal discovery method from interventional data. Each edge is learned separately using a score-based approach that relies on Monte-Carlo sampling. The main difference is that - since edges are learned separately - acyclicity need not be enforced. ", "summary_of_the_review": "A new approach for causal discovery and some interesting ideas but the utility should be demonstrated by comparing with simple baselines since interventional data is assumed to be available. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "LEiNPkIuwx9", "reviewer": "Reviewer_Paee", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "The paper is in general well-written. The discussion on the latent confounder case and few interventions are appreciated. The strength of the work is its originality. It introduces a certain parameterization method for enforcing the acyclicity of the model without imposing acyclic constraints on the optimization problem.\n\nNevertheless, I do have some concerns about the work:\n1. The sparsity penalty term in (2).\nThe graph with the maximum likelihood may not be the one with the minimum loss value of (2) because of the penalty term. The graph fitting procedure is actually not picking up the graph with the maximum likelihood but selecting a graph by jointly considering likelihood and the penalty term. So why would this be reasonable?\n\n2. The convergence guarantee.\n* It requires intervention on all variables. It would be necessary to justify the significancy of the convergence result. According to the citation (Eberhardt et al., 2005) in the paper, roughly speaking, in the worst case, it requires log2(N) + 1 experiments to determine all causal relations, where N is the number of variables. \n\n* Moreover, in the case with intervention on fewer variables, it remains unclear that what would be the results in such case. Are they in the same Markov equivalent class as the ground-truth graph, or there can be arbitrary differences from the ground-truth graph?\n\n* It would also be necessary to justify condition 3 in Thm. 3.1. Especially, how realistic is the condition? When can it be satisfied? And what will happen if it is not satisfied? \nIt would be better to make it clear whether the improper sparsity parameter can only lead to slow convergence or it can even lead to wrong results.\n\n3. A minor suggestion. It would be better if the paper distinguishes causal discovery on observational from causal discovery on both observational and interventional data, and then introduce more about the later setting which is more related to the work.\nBecause when mentioning causal discovery, in common, it refers to only on the observational data, it can be misleading sometimes.\nMoreover, instead of only generally mentioning the works based on observational and interventional data, it would be helpful for having an overview and a better judgment and understanding of the paper if including the important theoretical results and the related methods/conclusions in such setting.", "summary_of_the_paper": "The paper works on causal discovery on both observational and interventional data.\nIt proposes a method for enforcing the acyclicity of DAGs which avoids imposing constraints on the optimization problem.\nIt proves that under some assumptions, when all the variables have interventional data, the resulting graph converges to the ground-truth one. The experimental results show the properties and performance of the proposed method, including the cases with latent confounders and interventions on fewer variables.", "summary_of_the_review": "The strength of the paper is the proposed way for enforcing acyclicity and the good properties of the proposed model. \nHowever, the convergence guarantee may require justification for its significancy and Condition 3. Moreover, the loss function in the graph fitting procedure may need a justification and discussion for the penalty term.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the problem of learning a graphical model given observational and experimental data. The main novelty is the use of interventions to avoid the acyclicity constraint that plagues existing methods. Although this idea is quite standard and well-known, the generality of the approach merits consideration. After the discussion, there was a consensus among the reviewers to accept this paper. Some valid concerns have been raised and we expect that the authors will take into account all of the suggestions raised by the reviewers.
Decision: Accept (Poster)