# Efficient Neural Causal Discovery
without Acyclicity Constraints

Phillip Lippe
University of Amsterdam
QUVA lab
<EMAIL_ADDRESS>
&Taco Cohen
Qualcomm AI Research
<EMAIL_ADDRESS>
&Efstratios Gavves
University of Amsterdam
QUVA lab
<EMAIL_ADDRESS>
Qualcomm AI Research is an initiative of Qualcomm Technologies, Inc.

###### 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
when interventions on all variables are available, 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.

##  1 Introduction

Uncovering and understanding causal mechanisms is an 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 is causal structure learning Pearl (2009); Peters et al. (2017),
which aims at learning a directed acyclic graph (DAG) in which edges represent
causal relations between variables. While observational data alone is in
general not sufficient to identify the DAG Yang et al. (2018); 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 to a 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).

A promising new direction for scaling up DAG discovery methods are continuous-
optimization 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-based Peters et al. (2017); Guo et al. (2020) methods,
continuous-optimization methods reinterpret the search over discrete graph
topologies as a continuous problem with neural networks as function
approximators, for which efficient solvers are 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 to apply a regularizer Zhu et al. (2020); Ke et
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 to difficulties in enforcing acyclicity.

In this work, we address both problems. By modeling the orientation of an edge
as a 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, if we are
able to intervene on all variables, the proposed optimization is guaranteed to
converge to the correct, acyclic 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 _E fficient 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 as well.
Fourthly, we extend ENCO to detecting latent confounders. In various
experimental settings, ENCO recovers graphs accurately, making less than one
error 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 $\mathbb{P}$ over
a set of random variables $\bm{X}=\left\\{X_{1},...,X_{N}\right\\}$ and a
directed acyclic graph (DAG) $G=(V,E)$. Each node $i\in V$ corresponds to the
random variable $X_{i}$ and an edge $(i,j)\in E$ represents a direct causal
relation from variable $X_{i}$ to $X_{j}$: $X_{i}\to X_{j}$. A joint
distribution $\mathbb{P}$ is faithful to a graph $G$ if all and only the
conditional independencies present in $\mathbb{P}$ are entailed by $G$ Pearl
(1988). Vice versa, a distribution $\mathbb{P}$ is Markov to a graph $G$ if
the joint distribution can be factorized as
$p(\bm{X})=\prod_{i=1}^{N}p_{i}(X_{i}|\text{pa}(X_{i}))$ where
$\text{pa}(X_{i})$ 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 a variable $X_{i}$ describes the local
replacement of its conditional distribution $p_{i}(X_{i}|\text{pa}(X_{i}))$ by
a new distribution $\tilde{p}(X_{i}|\text{pa}(X_{i}))$. $X_{i}$ is thereby
referred to as the intervention target. An intervention is “perfect” when the
new distribution is independent of the original parents, _i.e._
$\tilde{p}(X_{i}|\text{pa}(X_{i}))=\tilde{p}(X_{i})$.

###  2.2 Causal structure learning

Discovering the graph $G$ from samples of a joint distribution $\mathbb{P}$ is
called causal structure learning or causal discovery, a fundamental problem in
causality Pearl (2009); Peters et al. (2017). While often performed from
observational data, i.e. samples from $\mathbb{P}$ (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-based 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 et 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 et 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 is
to extend methods working on observations only to 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
Tsamardinos 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
log-likelihood, as well as regularizers for encouraging sparsity. Since the
search space of DAGs is super-exponential in the number of nodes, many methods
rely on a 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 a proposal graph 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.

Continuous-optimization methods are score-based methods that avoid the
combinatorial greedy search over DAGs by using gradient-based methods Zheng et
al. (2018); Ke et al. (2019); Lachapelle et al. (2020); Yu et al. (2019);
Zheng et al. (2020); Zhu et 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 acyclic graphs. One common
approach is to view the search as a 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)
propose to use a regularization term penalizing cyclic graphs while allowing
unconstrained optimization. However, the regularizer must be designed and
weighted such that the correct, acyclic causal graph is the global optimum of
the score function.

##  3 Efficient Neural Causal Discovery

Figure 1: Visualization of the two training stages of ENCO, distribution
fitting and graph fitting, on an example graph with 3 variables
($X_{1},X_{2},X_{3}$). The graph on the right further shows how the parameters
$\bm{\gamma}$ and $\bm{\theta}$ correspond to edge probabilities. We learn
those parameters by comparing multiple graph samples on how well they
generalize from observational to interventional data.

###  3.1 Scope and Assumptions

We consider the task of finding a directed acyclic graph $G=(V,E)$ 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 a different variable;
(3) The interventions are “perfect” and “stochastic”, meaning the intervention
does not set the variable necessarily to a single value. Thereby, we 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 to 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 a 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, we alternate between two learning stages,
that is distribution fitting and graph fitting, visually summarized in Figure
1.

Distribution fitting trains a neural network $f_{\phi_{i}}$ per variable
$X_{i}$ parameterized by $\phi_{i}$ to model its observational, conditional
data distribution $p(X_{i}|...)$. The input to the network are all other
variables, $\bm{X}_{-i}$. For simplicity, we want this neural network to model
the conditional of the variable $X_{i}$ 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_{j}$ to zero based on the
probability of its corresponding edge $X_{j}\to X_{i}$, and minimize

$\min_{\phi_{i}}\mathbb{E}_{\bm{X}}\mathbb{E}_{\bm{M}}\left[-\log
f_{\phi_{i}}(X_{i};\bm{M}_{-i}\odot\bm{X}_{-i})\right],$ (1)

where $M_{j}\sim\text{Ber}(p(X_{j}\to X_{i}))$. For categorical random
variables $X_{i}$, we apply a softmax output activation for $f_{\phi_{i}}$,
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: $\bm{\gamma}\in\mathbb{R}^{N\times N}$ represents the
existence of edges in a graph, and $\bm{\theta}\in\mathbb{R}^{N\times N}$ the
orientation of the edges. The likelihood of an edge is determined by
$p(X_{i}\to X_{j})=\sigma(\gamma_{ij})\cdot\sigma(\theta_{ij})$, with
$\sigma(...)$ being the sigmoid function and $\theta_{ij}=-\theta_{ji}$. The
probability of the two orientations always sum to one. The benefit of
separating the edge probabilities into two independent parameters
$\bm{\gamma}$ and $\bm{\theta}$ is that it gives us 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_{i}$ or $X_{j}$ for learning $\theta_{ij}$. While
other interventions eventually provide information on the edge direction,
e.g., intervening on a node $X_{k}$ which is a child of $X_{i}$ and a parent
of $X_{j}$, we do not know the relation of $X_{k}$ to $X_{i}$ and $X_{j}$ at
this stage, as we are in the process of learning the structure. Despite having
just one variable for the orientation, $\gamma_{ij}$ and $\gamma_{ji}$ are
learned as two separate parameters. One reason is that on interventional data,
an edge can improve the log-likelihood estimate in one direction, but not
necessarily the other, leading to conflicting gradients.

We optimize the graph parameters $\bm{\gamma}$ and $\bm{\theta}$ by minimizing

$\tilde{\mathcal{L}}=\mathbb{E}_{\hat{I}\sim
p_{I}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\bm{X})}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]+\lambda_{\text{sparse}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sigma(\gamma_{ij})\cdot\sigma(\theta_{ij})$
(2)

where $p_{I}(I)$ is the distribution over which variable to intervene on
(usually uniform), and $\tilde{p}_{\hat{I}}(\bm{X})$ the joint distribution of
all variables under the intervention $\hat{I}$. In other words, these two
distributions represent our interventional data distribution. With
$p_{\bm{\gamma},\bm{\theta}}(C)$, we denote the distribution over adjacency
matrices $C$ under $\bm{\gamma},\bm{\theta}$, where
$C_{ij}\sim\text{Ber}(\sigma(\gamma_{ij})\sigma(\theta_{ij}))$.
$\mathcal{L}_{C}(X_{i})$ is the negative log-likelihood estimate of variable
$X_{i}$ conditioned on the parents according to $C$:
$\mathcal{L}_{C}(X_{i})=-\log
f_{\phi_{i}}(X_{i};C_{\cdot,i}\odot\bm{X}_{-i})$. The second term of Equation
2 is an $\ell_{1}$-regularizer on the edge probabilities. It acts as a prior,
selecting the sparsest graph of those with similar likelihood estimates by
removing redundant edges.

Prediction. Alternating between the distribution and graph fitting stages
allows us 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 $\sigma(\gamma_{ij})$ and $\sigma(\theta_{ij})$ are
greater than 0.5. The orientation parameters prevent loops between any two
variables, since $\sigma(\theta_{ij})$ 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 $\bm{\gamma}$ and $\bm{\theta}$ based on Equation 2, we need to
determine their gradients through the expectation
$\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}$, where $C$ is a 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 $\gamma_{ij}$, we obtain the
following gradient:

$\frac{\partial}{\partial\gamma_{ij}}\tilde{\mathcal{L}}=\sigma^{\prime}(\gamma_{ij})\cdot\sigma(\theta_{ij})\cdot\mathbb{E}_{\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to
X_{j}}(X_{j})+\lambda_{\text{sparse}}\right]$ (3)

where $\mathbb{E}_{\bm{X},C_{-ij}}$ summarizes for brevity the three
expectations in Equation 2, excluding the edge $X_{i}\to X_{j}$ from $C$.
Further, $\mathcal{L}_{X_{i}\to X_{j}}(X_{j})$ denotes the negative log-
likelihood for $X_{j}$, if we include the edge $X_{i}\to X_{j}$ to the
adjacency matrix $C_{-ij}$, _i.e._ , $C_{ij}=1$, and
$\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})$ if we set $C_{ij}=0$. The gradient
in Equation 3 can be intuitively explained: if by the addition of the edge
$X_{i}\to X_{j}$, the log-likelihood estimate of $X_{j}$ is improved by more
than $\lambda_{\text{sparse}}$, we increase the corresponding edge parameter
$\gamma_{ij}$; otherwise, we decrease it.

We derive the gradients for the orientation parameters $\bm{\theta}$
similarly. As mentioned before, we only take the gradients for $\theta_{ij}$
when we perform an intervention on either $X_{i}$ or $X_{j}$. This leads us
to:

$\begin{split}\frac{\partial}{\partial\theta_{ij}}\tilde{\mathcal{L}}=\sigma^{\prime}(\theta_{ij})\Big{(}&p(I_{X_{i}})\cdot\sigma(\gamma_{ij})\cdot\mathbb{E}_{I_{X_{i}},\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})\right]-\\\\[-5.0pt]
&p(I_{X_{j}})\cdot\sigma(\gamma_{ji})\cdot\mathbb{E}_{I_{X_{j}},\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{j}\to
X_{i}}(X_{i})-\mathcal{L}_{X_{j}\not\to
X_{i}}(X_{i})\right]\Big{)}\end{split}$ (4)

The probability of taking an intervention on $X_{i}$ is represented by
$p(I_{X_{i}})$ (usually uniform across variables), and
$\mathbb{E}_{I_{X_{i}},\bm{X},C_{-ij}}$ the same expectation as before under
the intervention on $X_{i}$. When the oriented edge $X_{i}\to X_{j}$ improves
the log-likelihood of $X_{j}$ under intervention on $X_{i}$, then the first
part of the gradient increases $\theta_{ij}$. In contrast, when the true edge
is $X_{j}\to X_{i}$, the correlation between $X_{i}$ and $X_{j}$ learned from
observational data would yield a worse likelihood estimate of $X_{j}$ on
interventional data on $X_{i}$ than without the edge $X_{j}\to X_{i}$. This is
because $p(X_{j}|X_{i},...)$ does not stay invariant under intervening on
$X_{i}$. The same dynamic holds for interventions on $X_{j}$. Lastly, for
independent nodes, the expectation of the gradient is zero.

Based on 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 $p_{\bm{\gamma},\bm{\theta}}(C)$
($K$ usually between 20 and 100). We then evaluate the log likelihoods of all
variables for these graphs on the batch, and estimate $\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})$ and $\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})$ for all pairs of
variables $X_{i}$ and $X_{j}$ by simply averaging the results for the two
cases separately. Finally, the estimates are used to determine the gradients
for $\bm{\gamma}$ and $\bm{\theta}$.

Low variance. Previous methods Ke et al. (2019); Bengio et al. (2020) relied
on a different REINFORCE-like estimator proposed by Bengio et al. (2020).
Adjusting their estimator to our setting of the parameter $\gamma_{ij}$, for
instance, the gradient looks as follows:

$g_{ij}=\sigma(\theta_{ij})\cdot\mathbb{E}_{\bm{X}}\left[\frac{\mathbb{E}_{C}\left[(\sigma(\gamma_{ij})-C_{ij})\cdot\mathcal{L}_{C}(X_{j})\right]}{\mathbb{E}_{C}\left[\mathcal{L}_{C}(X_{j})\right]}+\lambda_{\text{sparse}}\right]$
(5)

Figure 2: ENCO estimates gradients of significantly lower variance compared to
Bengio et al. (2020).

where $g_{ij}$ represents the gradient of $\gamma_{ij}$. Performing Monte-
Carlo sampling for estimating the gradient leads to a biased estimate which
becomes asymptotically unbiased with increasing number of samples Bengio et
al. (2020). The division by the expectation of $\mathcal{L}_{C}(X_{j})$ is
done for variance reduction Mahmood et al. (2014). Equation 5, however, is
still sensitive to the proportion of sampled $C_{ij}$ being one or zero. A
major benefit of our gradient formulation in Equation 3, instead, is that it
removes this noise by considering the difference of the two independent Monte-
Carlo estimates $\mathcal{L}_{X_{i}\to X_{j}}(X_{j})$ and
$\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})$. Hence, we can use a smaller sample
size than previous methods and attain 10 times lower standard deviation, as
visualized in Figure 2.

###  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 a 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 $p(\bm{X})$ is
Markovian with respect to the true graph $\mathcal{G}$. In other words, the
parent set $\text{pa}(X_{i})$ reflects the inputs to the causal generation
mechanism of $X_{i}$. We assume that there exist no latent confounders in
$\mathcal{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:

######  Theorem 3.1.

Given a causal graph $\mathcal{G}$ with variables $X_{1},...,X_{N}$ and
conditional observational distributions $p(X_{i}|...)$, the proposed method
ENCO will converge to the true, causal graph $\mathcal{G}$, if the following
conditions hold for all edges $X_{i}\to X_{j}$ in $\mathcal{G}$:

  1. 1.

For all possible sets of parents of $X_{j}$ excluding $X_{i}$, by adding
$X_{i}$ the log-likelihood estimate of $X_{j}$ is improved or unchanged under
the intervention on $X_{i}$:

$\forall\widehat{\text{pa}}(X_{j})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{i}},\mathbf{X}}\left[\log
p(X_{j}|\widehat{\text{pa}}(X_{j}),X_{i})-\log
p(X_{j}|\widehat{\text{pa}}(X_{j}))\right]\geq 0$ (6)

  2. 2.

For at least one set of nodes $\widehat{\text{pa}}(X_{j})$, for which the
probability to be sampled as parents of $X_{j}$ is greater than 0, Equation 6
must be strictly greater than zero.

  3. 3.

The effect of $X_{i}$ on $X_{j}$ cannot be described by other variables up to
$\lambda_{\text{sparse}}$:

$\min_{\hat{\text{pa}}\subseteq\text{gpa}_{i}(X_{j})}\mathbb{E}_{\hat{I}\sim
p_{I_{-j}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{j}|\hat{\text{pa}},X_{i})-\log
p(X_{j}|\hat{\text{pa}})\big{]}>\lambda_{\text{sparse}}$ (7)

where $\text{gpa}_{i}(X_{j})$ is the set of nodes excluding $X_{i}$ which,
according to the ground truth graph, could have an edge to $X_{j}$ without
introducing a cycle, and $p_{I_{-j}}(I)$ refers to the distribution over
interventions $p_{I}(I)$ excluding the intervention on variable $X_{j}$.

Further, for all other pairs $X_{i},X_{j}$ for which $X_{j}$ is a descendant
of $X_{i}$, condition 1 and 2 must hold.

Condition 1 and 2 ensure that the orientations can be learned from
interventions. Intuitively, ancestors and descendants in the graph have to be
dependent when intervening on the ancestors. This aligns with the technical
interpretation in Theorem 3.1 that the likelihood estimate of the child
variable must improve when intervening and conditioning on its ancestor
variables. Condition 3 states intuitively 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 a gradient-based algorithm.
We emphasize that this condition only gives an upper bound for
$\lambda_{\text{sparse}}$ when sufficiently large datasets are available. In
practice, the graph can thus be recovered with a sufficiently small sparsity
regularizer and dependencies among variables under interventions. We provide
more details for various settings and further intuition in Appendix B.

Interventions on fewer variables. It is straightforward to 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 $\theta_{ij}$ of
only those edges that connect to an intervened variable, either $X_{i}$ or
$X_{j}$, as before. For all other orientation parameters, we extend the
gradient estimator to 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\in S_{N}$, with $S_{N}$ being the set of permutations, that
maximizes the pairwise orientation probabilities:
$\arg\max_{O}\prod_{i=1}^{N}\prod_{j=i+1}^{N}\sigma(\theta_{O_{i},O_{j}})$.
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 is 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, in the context of ENCO latent confounders for two
variables $X_{i},X_{j}$ cause a unique pattern of learned parameters. When
intervening on $X_{i}$ or $X_{j}$, having an edge between the two variables is
disadvantageous, as in the intervened graph $X_{i}$ and $X_{j}$ are
(conditionally) independent. For interventions on all other variables,
however, an edge can be beneficial as $X_{i}$ and $X_{j}$ 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 a child of any other
observed variable. For all other edges besides between $X_{i}$ and $X_{j}$, we
can still rely on the guarantees in Section 3.4 since Equation 7 already
includes the possibility of additional edges in such cases. After convergence,
we score every pair of variables on how likely they share a latent confounder
using a function $\text{lc}(\cdot)$ that is maximized in the scenario
mentioned above. For this, we define
$\gamma_{ij}=\gamma^{(I)}_{ij}+\gamma^{(O)}_{ij}$ where $\gamma^{(I)}_{ij}$ is
only updated with gradients from Equation 3 under interventions on $X_{i}$,
and $\gamma^{(O)}_{ij}$ on all others. With this separation, we define the
following score function which is maximized by latent confounders:

$\text{lc}(X_{i},X_{j})=\sigma\left(\gamma^{(O)}_{ij}\right)\cdot\sigma\left(\gamma^{(O)}_{ji}\right)\cdot\left(1-\sigma\left(\gamma^{(I)}_{ij}\right)\right)\cdot\left(1-\sigma\left(\gamma^{(I)}_{ji}\right)\right)$
(8)

To converge to the mentioned values, especially of $\gamma^{(O)}_{ij}$, we
need a similar condition as in Equation 7: the improvement on the log-
likelihood estimate gained by the edge $X_{i}\to X_{j}$ and conditioned on all
other parents of $X_{j}$ needs to be larger than $\lambda_{\text{sparse}}$ on
interventional data excluding $X_{i}$ and $X_{j}$. If this is not the case,
the sparsity regularizer will instead remove the edge between $X_{i}$ and
$X_{j}$ preventing any false predictions among observed variables. For all
other pairs of variables, at least one of the terms in Equation 8 converges to
zero. Thus, we can detect latent confounders by checking whether the score
function $\text{lc}(X_{i},X_{j})$ is greater than a threshold hyperparameter
$\tau\in(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://github.com/phlippe/ENCO.

###  4.1 Experimental setup

Graphs and datasets. Given a ground-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 et 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 a common observational baseline, we apply GES Chickering
(2002) on the observational data to obtain a graph skeleton, and orient each
edge by learning the skeleton on the corresponding interventional
distribution. We perform a 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

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$ ($\pm 7.5$) |  $17.5$ ($\pm 7.3$) |  $24.0$ ($\pm 2.9$) |  $216.5$ ($\pm 15.2$) |  $33.1$ ($\pm 2.9$) |  $57.5$ ($\pm 14.2$)
IGSP |  $32.7$ ($\pm 5.1$) |  $14.6$ ($\pm 2.3$) |  $23.7$ ($\pm 2.3$) |  $253.8$ ($\pm 12.6$) |  $35.9$ ($\pm 5.2$) |  $65.4$ ($\pm 8.0$)
GES + Orientating |  $14.8$ ($\pm 2.6$) |  $0.5$ ($\pm 0.7$) |  $20.8$ ($\pm 2.4$) |  $282.8$ ($\pm 4.2$) |  $14.7$ ($\pm 3.1$) |  $60.1$ ($\pm 8.9$)
SDI |  $9.0$ ($\pm 2.6$) |  $3.9$ ($\pm 2.0$) |  $16.1$ ($\pm 2.4$) |  $153.9$ ($\pm 10.3$) |  $6.9$ ($\pm 2.3$) |  $10.8$ ($\pm 3.9$)
DCDI |  $16.9$ ($\pm 2.0$) |  $10.1$ ($\pm 1.1$) |  $10.9$ ($\pm 3.6$) |  $21.0$ ($\pm 4.8$) |  $17.9$ ($\pm 4.1$) |  $7.7$ ($\pm 3.2$)
ENCO (ours) |  $2.2$ ($\pm 1.4$) |  $1.7$ ($\pm 1.3$) |  $\bm{1.6}$ ($\pm 1.6$) |  $9.2$ ($\pm 3.4$) |  $1.7$ ($\pm 1.3$) |  $4.6$ ($\pm 1.9$)
ENCO-acyclic (ours) |  $\bm{0.0}$ ($\pm 0.0$) |  $\bm{0.0}$ ($\pm 0.0$) |  $\bm{1.6}$ ($\pm 1.6$) |  $\bm{5.3}$ ($\pm 2.3$) |  $\bm{0.6}$ ($\pm 1.1$) |  $\bm{0.2}$ ($\pm 0.5$)

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 chain and full represent the minimally- and maximally-connected DAGs.
The graph bidiag is a chain with 2-hop connections, and jungle is a tree-like
graph. In the collider graph, one node has all other nodes as parents.
Finally, 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 1 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 on 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 to recover four out
of six graphs with less than one error on average. We further include
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.

###  4.3 Scalability

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.

Next, we test ENCO on graphs with large sets of variables. We create random
graphs ranging from $N=100$ to $N=1$,$000$ nodes with larger sample sizes.
Every node has on average 8 edges and a 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.

We compare ENCO to the two best performing baselines from Table 1, SDI and
DCDI. All methods 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 errors except for the $1$,$000$-node graph, for which it misses only one
out of 1 million edges in 2 out 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.

###  4.4 Interventions on fewer 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.

We perform experiments on 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 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.4 for detailed discussion). We
conclude that ENCO works competitively with partial interventions too.

###  4.5 Detecting latent confounders

Table 2: Results of ENCO on detecting latent confounders. The missed confounders do not affect other edge predictions.  Metrics | ENCO
---|---
SHD |  $0.0$ ($\pm 0.0$)
Confounder recall |  $96.8\%$ ($\pm 9.5\%$)
Confounder precision |  $100.0\%$ ($\pm 0.0\%$)

To test the detection of latent confounders, we create a set of 25 random
graphs with 5 additional latent confounders. The dataset is generated in the
same way as before, except that we remove the latent variable from the input
data and increase the observational and interventional sample size (see
Appendix C.3 for ablation studies). After training, we predict the existence
of a latent confounder on any pair of variables $X_{i}$ and $X_{j}$ if
$\text{lc}(X_{i},X_{j})$ is greater than $\tau$. We choose $\tau=0.4$ but
verify in Appendix C.3 that the method is not sensitive to the specific value
of $\tau$. As shown in Table 2, ENCO 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 are recovered perfectly.

###  4.6 Real-world inspired data

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 | asia | sachs | child | alarm | diabetes | pigs
---|---|---|---|---|---|---|---
| (5 nodes) | (8 nodes) | (11 nodes) | (20 nodes) | (37 nodes) | (413 nodes) | (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) | $\mathbf{0.0}$ | $\mathbf{0.0}$ | $\mathbf{0.0}$ | $\mathbf{0.0}$ | $\mathbf{1.0}$ | $\mathbf{2.0}$ | $\mathbf{0.0}$

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 in literature. In
comparison to the synthetic graphs, the real-world graphs are sparser with a
maximum of 6 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 on 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 a
maximum of 2 parents per node, and hence is easier to learn. The most
challenging graph is diabetes Andreassen et 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 a wide variety of settings, including
deterministic variables.

##  5 Conclusion

We propose ENCO, an efficient causal structure learning method leveraging
observational and interventional data. Compared to previous work, ENCO models
the edge orientations as 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 a
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-art on an
extensive array of settings on 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 on 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 are missing
transitivity: if $X_{1}\succ X_{2}$ and $X_{2}\succ X_{3}$, then $X_{1}\succ
X_{3}$ must hold. A potential direction is incorporating transitive relations
for improving convergence speed and results on fewer interventions.

## 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 a desired way. For instance, by
knowing that the environmental pollution in a city has an impact on 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 are ranging
across many scientific disciplines, including computational biology Friedman
et al. (2000); Sachs et al. (2005); Opgen-Rhein & Strimmer (2007),
epidemiology Robins et 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 is the field of
genomics. Recent advances have enabled to perform gene knockdown experiments
in a 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 a first step towards this goal, and
our work can foster future work for creating algorithms scaling beyond
$10$,$000$ variables.

Since the possible applications are fairly wide-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 as 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 to use this method in a
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, we have published the source code of the proposed
method ENCO at https://github.com/phlippe/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 a 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 on a computer
with a 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 on smaller GPUs as 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
discussion which helped to improve this work further. Finally, we thank
SURFsara for the support in using the Lisa Compute Cluster.

## References

  * Andreassen et al. (1991) 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 Informatics_ , volume 44 of _Lecture Notes in Medical Informatics_ , pp. 239–249, Germany, 1991. Springer.  Conference date: 24-06-1991 Through 27-06-1991.
  * Beinlich et al. (1989) 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-93437-7.
  * Bengio et al. (2020) 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.
  * Brouillard et al. (2020) Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin.  Differentiable Causal Discovery from Interventional Data.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin (eds.), _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Chickering (2002) David Maxwell Chickering.  Optimal structure identification with greedy search.  _Journal of machine learning research_ , 3(Nov):507–554, 2002.
  * Christina et al. (2018) Heinze-Deml Christina, Meinshausen Nicolai, and Peters Jonas.  Invariant Causal Prediction for Nonlinear Models.  _Journal of Causal Inference_ , 6(2):1–35, 2018\.
  * Cover & Thomas (2005) Thomas M. Cover and Joy A. Thomas.  _Elements of Information Theory_.  John Wiley and Sons, Ltd, 2005.  ISBN 9780471748823.
  * Dixit et al. (2016) A. Dixit, O. Parnas, B. Li, J. Chen, C. P. Fulco, L. Jerby-Arnon, 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_ , 167(7):1853–1866, 2016.
  * Eberhardt (2008) Frederick Eberhardt.  Almost Optimal Intervention Sets for Causal Discovery.  In _Proceedings of the Twenty-Fourth Conference on Uncertainty in Artificial Intelligence_ , UAI’08, pp. 161–168, Arlington, Virginia, USA, 2008. AUAI Press.  ISBN 0974903949.
  * Eberhardt et al. (2005) Frederick Eberhardt, Clark Glymour, and Richard Scheines.  On the Number of Experiments Sufficient and in the Worst Case Necessary to Identify All Causal Relations among N Variables.  In _Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence_ , UAI’05, pp. 178–184, Arlington, Virginia, USA, 2005\. AUAI Press.  ISBN 0974903914.
  * Friedman et al. (2000) Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe’er.  Using Bayesian networks to analyze expression data.  _Journal of computational biology_ , 7(3-4):601–620, 2000.
  * Glymour et al. (2019) 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.
  * Goudet et al. (2017) Olivier Goudet, Diviyan Kalainathan, Philippe Caillou, Isabelle Guyon, David Lopez-Paz, and Michèle Sebag.  Causal generative neural networks.  _arXiv preprint arXiv:1711.08936_ , 2017.
  * Guo et al. (2020) Ruocheng Guo, Lu Cheng, Jundong Li, P. Richard Hahn, and Huan Liu.  A Survey of Learning Causality with Data: Problems and Methods.  _ACM Comput. Surv._ , 53(4), 2020.  ISSN 0360-0300.
  * Hauser & Bühlmann (2012) Alain Hauser and Peter Bühlmann.  Characterization and Greedy Learning of Interventional Markov Equivalence Classes of Directed Acyclic Graphs.  _Journal of Machine Learning Research_ , 13(1):2409–2464, 2012.  ISSN 1532-4435.
  * Hicks et al. (1980) John Hicks et al.  _Causality in economics_.  Australian National University Press, 1980.
  * Hornik et al. (1989) Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  _Neural Networks_ , 2(5):359–366, 1989.  ISSN 0893-6080.
  * Huang et al. (2018) Chin-Wei 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ässan, Stockholm, Sweden, July 10-15, 2018_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 2083–2092. PMLR, 2018\.
  * Hyttinen et al. (2014) Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo.  Constraint-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\.
  * Ivanov et al. (2019) 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.
  * Jaber et al. (2020) 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 Processing Systems_ , volume 33, pp. 9551–9561. Curran Associates, Inc., 2020.
  * Kalainathan et al. (2020) Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta.  Causal Discovery Toolbox: Uncovering causal relationships in Python.  _J. Mach. Learn. Res._ , 21:37–1, 2020.
  * Ke et al. (2019) 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.01075_ , 2019.
  * Kingma & Ba (2015) 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.
  * Kocaoglu et al. (2019) Murat Kocaoglu, Amin Jaber, Karthikeyan Shanmugam, and Elias Bareinboim.  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. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Korb & Nicholson (2010) Kevin B. Korb and Ann E. Nicholson.  _Bayesian artificial intelligence, second edition_.  CRC Press, Australia, 2010.  ISBN 9781439815915.
  * Lachapelle et al. (2020) 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.
  * Lauritzen & Spiegelhalter (1988) S. Lauritzen and D. Spiegelhalter.  Local computations with probabilities on graphical structures and their application to expert systems.  _Journal of the royal statistical society series b-methodological_ , 50:415–448, 1988.
  * Li et al. (2020) Yang Li, Shoaib Akbar, and Junier Oliva.  ACFlow: Flow Models for Arbitrary Conditional Likelihoods.  In _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 5831–5841. PMLR, 2020\.
  * Luce (1959) R Duncan Luce.  _Individual choice behavior_.  John Wiley, Oxford, England, 1959.
  * Macosko et al. (2015) E. Z. Macosko, A. Basu, R. Satija, J. Nemesh, K. Shekhar, M. Goldman, I. Tirosh, A. R. Bialas, N. Kamitaki, E. M. Martersteck, J. 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):1202–1214, 2015.
  * Mahmood et al. (2014) 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 Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada_ , pp. 3014–3022, 2014.
  * Meek (1997) Christopher Meek.  _Graphical Models: Selecting causal and statistical models_.  PhD thesis, Carnegie Mellon University, 1997.
  * Monti et al. (2019) Ricardo Pio Monti, Kun Zhang, and Aapo Hyvärinen.  Causal Discovery with General Non-Linear Relationships using Non-Linear ICA.  In Amir Globerson and Ricardo Silva (eds.), _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence, UAI 2019, Tel Aviv, Israel, July 22-25, 2019_ , volume 115 of _Proceedings of Machine Learning Research_ , pp. 186–195. AUAI Press, 2019.
  * Mooij et al. (2020) Joris M. Mooij, Sara Magliacane, and Tom Claassen.  Joint Causal Inference from Multiple Contexts.  _Journal of Machine Learning Research_ , 21(99):1–108, 2020.
  * Ng et al. (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.
  * Opgen-Rhein & Strimmer (2007) Rainer Opgen-Rhein and Korbinian Strimmer.  From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data.  _BMC Systems Biology_ , 1(1):37, 2007.  ISSN 1752-0509.
  * Paszke et al. (2019) 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 DeVito, 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’Alché-Buc, 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–8035, 2019.
  * Pearl (1988) Judea Pearl.  _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988.  ISBN 0934613737.
  * Pearl (2009) Judea Pearl.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, USA, 2nd edition, 2009.  ISBN 052189560X.
  * Peters & Bühlmann (2015) Jonas Peters and Peters Bühlmann.  Structural Intervention Distance (SID) for Evaluating Causal Graphs.  _Neural Computation_ , 27(3):771–799, 2015.
  * Peters et al. (2016) Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen.  Causal inference by using invariant prediction: identification and confidence intervals.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , pp. 947–1012, 2016.
  * Peters et al. (2017) Jonas Peters, Dominik Janzing, and Bernhard Schlkopf.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  The MIT Press, 2017.  ISBN 0262037319.
  * Plackett (1975) R. L. Plackett.  The Analysis of Permutations.  _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , 24(2):193–202, 1975.  ISSN 00359254, 14679876.
  * Rezende & Mohamed (2015) 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.
  * Robins et al. (2000) J.M. Robins, M.A. Hernan, and B Brumback.  Marginal Structural Models and Causal Inference in Epidemiology.  _Epidemiology_ , 11(5):550–560, 2000.
  * Sachs et al. (2005) K. Sachs, O. Perez, D. Pe’er, D. Lauffenburger, and G. Nolan.  Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data.  _Science_ , 308:523–529, 2005.
  * Schölkopf et al. (2012) 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_. icml.cc / Omnipress, 2012.
  * Schölkopf et al. (2021) Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio.  Toward causal representation learning.  _Proceedings of the IEEE_ , 109(5):612–634, 2021\.
  * Scutari (2010) Marco Scutari.  Learning Bayesian Networks with the bnlearn R Package.  _Journal of Statistical Software_ , 35(3):1–22, 2010.
  * Spiegelhalter & Cowell (1992) D. J. Spiegelhalter and R. G. Cowell.  Learning in probabilistic expert systems.  _Bayesian Statistics_ , 4:447–466, 1992.
  * Spirtes et al. (2000) 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-19440-2.
  * Sun et al. (2007) Xiaohai Sun, Dominik Janzing, Bernhard Schölkopf, and Kenji Fukumizu.  A kernel-based causal learning algorithm.  In Zoubin Ghahramani (ed.), _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.
  * Tsamardinos et al. (2006) 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.
  * Vandenbroucke et al. (2016) Jan P Vandenbroucke, Alex Broadbent, and Neil Pearce.  Causality and causal inference in epidemiology: the need for a pluralistic approach.  _International journal of epidemiology_ , 45(6):1776–1786, 2016.
  * Varando (2020) Gherardo Varando.  Learning DAGs without imposing acyclicity.  _arXiv preprint arXiv:2006.03005_ , 2020.
  * Wang et al. (2017) Yuhao Wang, Liam Solus, Karren D. Yang, and Caroline Uhler.  Permutation-based 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.
  * Williams (1992) Ronald J. Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  In _Machine Learning_ , pp. 229–256, 1992.
  * Yang et al. (2018) 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\.
  * Yu et al. (2019) 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 International 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\.
  * Zheng et al. (2018) 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, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pp. 9492–9503, 2018.
  * Zheng et al. (2020) 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, Italy]_ , volume 108 of _Proceedings of Machine Learning Research_ , pp. 3414–3425. PMLR, 2020\.
  * Zhu et al. (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.

Supplementary material
Efficient Neural Causal Discovery without Acyclicity Constraints

[sections]  Table of Contents [sections]l1

##  Appendix 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 we parameterize the graph by edge existence
parameters $\bm{\gamma}$ and orientation parameters $\bm{\theta}$. Our
objective is to optimize $\bm{\gamma}$ and $\bm{\theta}$ such that we maximize
the probability of interventional data, _i.e._ , data generated from the true
graphs under (arbitrary) interventions on single variables. Thereby, the
likelihood estimates have been trained on observational data only.
Additionally, we want to ensure that the graph is as sparse as possible to
prevent unnecessary connections. Thus, an $\ell_{1}$ regularizer is added on
top of the edge probabilities. The full objective can be written as follows:

$\tilde{\mathcal{L}}=\mathbb{E}_{\hat{I}\sim
p_{I}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\bm{X})}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]+\lambda_{\text{sparse}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sigma(\gamma_{ij})\sigma(\theta_{ij})$
(9)

where:

  * •

$N$ is the number of variables in the causal graph ($X_{1},...,X_{N}$);

  * •

$p_{I}(I)$ is the distribution over interventions that are performed. This
distribution can be set as a hyperparameter to weight certain interventions
higher than others. In our experiments, we assume it to be uniform across
interventions on variables;

  * •

$\tilde{p}_{\hat{I}}(\bm{X})$ is the joint distribution of all variables under
the intervention $\hat{I}$;

  * •

$p_{\bm{\gamma},\bm{\theta}}(C)$ is the distribution over adjacency matrices
$C$, which we model as a product of independent edge probabilities:
$p_{\bm{\gamma},\bm{\theta}}(C)=\prod_{i=1}^{N}\prod_{j=1,j\neq
i}^{N}\sigma(\gamma_{ij})\cdot\sigma(\theta_{ij})$;

  * •

$\mathcal{L}_{C}(X_{i})$ is the negative log-likelihood estimate of variable
$X_{i}$ under sampled adjacency matrix $C$: $\mathcal{L}_{C}(X_{i})=-\log
f_{\phi_{i}}(X_{i};C_{\cdot,i}\odot\bm{X}_{-i})$;

  * •

$\lambda_{\text{sparse}}$ is a hyperparameter representing the regularization
weight.

Based on this objective, we derive the gradient estimators for optimizing both
edge existence and orientation parameters.

###  A.1 Low-variance gradient estimator for edge parameters

In order to optimize the edge parameters via SGD, we need to determine the
gradient $\frac{\partial}{\partial\gamma_{ij}}\tilde{\mathcal{L}}$. Since
$\mathcal{L}$ consists of a 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,l$ as indices for the
parameter $\gamma_{kl}$ for which we determine the gradient, _i.e._ ,
$\frac{\partial}{\partial\gamma_{kl}}\tilde{\mathcal{L}}$, and $i,j$ as
indices for sums.

As a first step, we determine the gradients for the regularization term. Those
can be found by taking the derivative of the sigmoid:

$\frac{\partial}{\partial\gamma_{kl}}\lambda_{\text{sparse}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sigma(\gamma_{ij})\sigma(\theta_{ij})=\sigma(\gamma_{kl})\cdot(1-\sigma(\gamma_{kl}))\cdot\sigma(\theta_{kl})\lambda_{\text{sparse}}$
(10)

Thus, it is straight-forward to calculate for any edge parameter. In the
following, we use $\sigma^{\prime}(...)$ to abbreviate the derivate of the
sigmoid:
$\sigma^{\prime}(\gamma_{kl})=\sigma(\gamma_{kl})(1-\sigma(\gamma_{kl}))$.

For the log-likelihood term, we start by reorganizing the expectations to
simplify the gradient expression. The derivate term
$\frac{\partial}{\partial\gamma_{kl}}$ can be moved inside the two
expectations over interventional data since those are independent of the graph
parameters. Thus, we can write:

$\frac{\partial}{\partial\gamma_{kl}}\tilde{\mathcal{L}}^{\prime}=\mathbb{E}_{\hat{I}\sim
p_{I}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\bm{X})}\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]$
(11)

For readability, we denote $\tilde{\mathcal{L}}^{\prime}$ to be the objective
in Equation 9 without the regularizer.

Next, we take a closer look at the derivate of the expectation over adjacency
matrices. Note that we have defined the adjacency matrix distribution as
$p_{\bm{\gamma},\bm{\theta}}(C)=\prod_{i=1}^{N}\prod_{j=1,j\neq
i}^{N}\sigma(\gamma_{ij})\sigma(\theta_{ij})$, with $C_{ij}=1$ representing
the edge $X_{i}\to X_{j}$. Since a parameter $\gamma_{ij}$ only influences the
likelihood of the edge $X_{i}\to X_{j}$ and no other edges, we can reduce the
expectation to a single binary variable over which we need to differentiate
the expectation:

$\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]=\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C_{-kl})}\left[\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\gamma,\theta}(C_{kl})}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]\right]$
(12)

where $p_{\gamma,\theta}(C_{kl})=\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})$.
The first expectation over $p_{\bm{\gamma},\bm{\theta}}(C_{-kl})$ is
independent of $\gamma_{kl}$ as we have defined the adjacency matrix
distribution to be a product of independent edge probabilities.

The log-likelihood estimate of a variable, $\mathcal{L}_{C}(X_{i})$, depends
on the adjacency matrix column $C_{\cdot,i}$ which represents the input
connections to the node $X_{i}$. All other edges have no influence on the log-
likelihood estimate of $X_{i}$. Hence, the parameter $\gamma_{kl}$ only
influences $\mathcal{L}_{C}(X_{l})$, and thus we can reduce the sum inside the
expectation to:

$\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\gamma,\theta}(C_{kl})}\left[\sum_{i=1}^{N}\mathcal{L}_{C}(X_{i})\right]=\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\gamma,\theta}(C_{kl})}\left[\mathcal{L}_{C}(X_{l})\right]$
(13)

Figure 5: Visualizing the gradient calculation for the incoming edges of
$X_{2}$ in an example graph with three variables. The intervention is being
performed on $X_{1}$, and the data is used to calculate the log-likelihood
estimates under the three randomly sampled graphs:
$\mathcal{L}_{C^{1}}(X_{2}),\mathcal{L}_{C^{2}}(X_{2})$ and
$\mathcal{L}_{C^{3}}(X_{2})$. Those terms are assigned to the Monte-Carlo
estimators for $\mathcal{L}_{X_{i}\to X_{2}}(X_{2})$ and
$\mathcal{L}_{X_{i}\not\to X_{2}}(X_{2})$, and finally used to determine the
gradients for $\bm{\gamma}$ and $\bm{\theta}$. The same process is performed
for $X_{3}$ as well.

The REINFORCE trick is a simple method to move the derivative of a discrete
distribution inside the expectation. Applied to our situation, we obtain:

$\frac{\partial}{\partial\gamma_{kl}}\mathbb{E}_{p_{\gamma,\theta}(C_{kl})}\left[\mathcal{L}_{C}(X_{l})\right]=\mathbb{E}_{p_{\gamma,\theta}(C_{kl})}\left[\mathcal{L}_{C}(X_{l})\frac{\partial\log
p_{\gamma,\theta}(C_{kl})}{\partial\gamma_{kl}}\right]$ (14)

This leaves us with two cases in the expectation: $C_{kl}=0$ and $C_{kl}=1$.
In other words, we need to distinguish between samples of $C$ where we have
the edge $X_{k}\to X_{l}$, and where we do not have such an edge
($X_{k}\not\to X_{l}$). Thus, we can also write the expectation as a weighted
sum of those two cases:

$\begin{split}\mathbb{E}_{p_{\gamma,\theta}(C)}\left[\mathcal{L}_{C}(X_{l})\frac{\partial\log
p_{\gamma,\theta}(C_{kl})}{\partial\gamma_{kl}}\right]&=\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})\cdot\mathcal{L}_{X_{k}\to
X_{l}}(X_{l})\cdot\frac{\partial\log\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})}{\partial\gamma_{kl}}+\\\
&\hskip
14.22636pt(1-\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl}))\cdot\mathcal{L}_{X_{l}\not\to
X_{k}}(X_{k})\cdot\frac{\partial\log\left(1-\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})\right)}{\partial\gamma_{kl}}\end{split}$
(15)

We use $\mathcal{L}_{X_{k}\to X_{l}}(X_{l})$ to denote the (expected) negative
log likelihood for $X_{l}$ under adjacency matrices where we have an edge from
$X_{k}$ to $X_{l}$:

$\displaystyle\mathcal{L}_{X_{k}\to X_{l}}(X_{l})$
$\displaystyle=\mathbb{E}_{p_{\gamma,\theta}(C_{-kl}),C_{kl}=1}\left[\mathcal{L}_{C}(X_{l})\right]$
(16) $\displaystyle\mathcal{L}_{X_{k}\not\to X_{l}}(X_{l})$
$\displaystyle=\mathbb{E}_{p_{\gamma,\theta}(C_{-kl}),C_{kl}=0}\left[\mathcal{L}_{C}(X_{l})\right]$
(17)

The final step is to solve the two derivative terms in Equation 14. This is
done as follows:

$\displaystyle\frac{\partial\log\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})}{\partial\gamma_{kl}}$
$\displaystyle=\frac{\partial\log\sigma(\gamma_{kl})+\log\sigma(\theta_{kl})}{\partial\gamma_{kl}}=1-\sigma(\gamma_{kl})$
(18)
$\displaystyle\frac{\partial\log\left(1-\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})\right)}{\partial\gamma_{kl}}$
$\displaystyle=-\frac{\sigma(\gamma_{kl})\cdot(1-\sigma(\gamma_{kl}))\cdot\sigma(\theta_{kl})}{1-\sigma(\gamma_{kl})\cdot\sigma(\theta_{kl})}$
(19)

Putting these results back in the original equation and adding the sparsity
regularizer, we get:

$\frac{\partial}{\partial\gamma_{ij}}\tilde{\mathcal{L}}=\sigma(\gamma_{ij})\cdot(1-\sigma(\gamma_{ij}))\cdot\sigma(\theta_{ij})\cdot\mathbb{E}_{\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to
X_{j}}(X_{j})+\lambda_{\text{sparse}}\right]$ (20)

To align the result with the gradient in Section 3.3, we switch the index
notation from $k,l$ to $i,j$ again. The expectation
$\mathbb{E}_{\bm{X},C_{-ij}}$ is a short form for the expectations
$\mathbb{E}_{\hat{I}\sim
p_{I}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\bm{X})}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C_{-ij})}$.
From this expression, we can see that the gradients of $\gamma_{ij}$ are
proportional to the difference of the expected negative log-likelihood of
$X_{j}$ with having an edge between $X_{i}\to X_{j}$, and the cases where
$X_{i}\not\to X_{j}$. The sparsity regularizer thereby biases the difference
towards the no-edge case. The value of $\gamma_{ij}$ and $\theta_{ij}$ 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 $\mathcal{L}_{X_{i}\to X_{j}}(X_{j})$
and $\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})$ are not required to be the same.
For efficiency, we instead sample $K$ adjacency matrices from
$p_{\bm{\gamma},\bm{\theta}}(C)$, evaluate the likelihood of a batch $\bm{X}$
under all these graphs. Afterwards, we assign the evaluated samples to one of
the two cases, depending on $C_{ij}$ being zero or one. This way, we can reuse
the same graph samples for all edge parameters $\bm{\gamma}$. We visualize the
gradient calculation in Figure 5. In the cases where we perform an
intervention on $X_{i}$, we do not optimize $\gamma_{ij}$ 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 log-likelihood estimates.

####  A.1.1 Comparison to previous gradient estimators

Figure 6: Plotting the gradient estimate variance of ENCO for three variables
of $\gamma$ compared to previous REINFORCE-like approach by Bengio et al.
(2020) on an example graph with $K=100$. The gradients have been scaled to
match in terms of averages. We can see a clear reduction in variance with the
gradient estimator of ENCO allowing us to use lower sample sizes.

As discussed in Section 3, previous work on similar structure learning methods
Bengio et al. (2020); Ke et al. (2019) relied on a different estimator. In
terms of derivation, the main difference is the continuation from Equation 14
on. In our proposed method, we write the expectation as the sum of two terms
that can independently be approximated via Monte-Carlo sampling. In
comparison, Bengio et al. (2020) proposed to 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) to
which we have experimentally compared our method.

Figure 2 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 it 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 random graphs from our experiments on synthetic graphs (see Section
4.2), and evaluated $64,000$ sampled adjacency matrices in terms of log-
likelihood estimates. These $64,000$ samples are grouped into sets of $K$
samples which we have used to estimate the gradients of $\bm{\gamma}$. 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).

To verify that the improvement of ENCO is not just because of the gradient
estimators, we 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 $\bm{\theta}$, we can
mostly follow the same approach as for the edge existence parameters
$\bm{\gamma}$. However, we have to keep in mind the constraint
$\theta_{kl}=-\theta_{lk}$ which ensures that the orientation probability sums
to one: $\sigma(\theta_{kl})+\sigma(\theta_{lk})=1$.

To determine the gradient of the likelihood term, we can separate the two
gradients of $\theta_{kl}$ and $\theta_{lk}$. This is because $\theta_{kl}$
only influences the expectation over $\mathcal{L}_{C}(X_{l})$, while
$\theta_{lk}$ concerns $\mathcal{L}_{C}(X_{k})$. We can follow Equation 11 to
Equation 20 of Section A.1 by swapping $\theta_{kl}$ and $\gamma_{kl}$. For
the derivative through the expectation, we obtain the following gradient:

$\frac{\partial}{\partial\theta_{kl}}\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C)}\left[\mathcal{L}_{C}(X_{l})\right]=\sigma^{\prime}(\theta_{kl})\cdot\sigma(\gamma_{kl})\cdot\mathbb{E}_{p_{\bm{\gamma},\bm{\theta}}(C_{-kl})}\left[\mathcal{L}_{X_{k}\to
X_{l}}(X_{l})-\mathcal{L}_{X_{k}\not\to X_{l}}(X_{l})\right]$ (21)

Since we have the condition that $\theta_{kl}=-\theta_{lk}$, the full gradient
for $\theta_{kl}$ would therefore consist of the gradient above minus the
gradient of Equation 21 with respect to $\theta_{ji}$. However, as 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 $\theta_{kl}$ if we intervene on node $X_{k}$, which gives us the following
gradient expression:

$\begin{split}\frac{\partial}{\partial\theta_{ij}}\tilde{\mathcal{L}}=\sigma^{\prime}(\theta_{ij})\Big{(}&p(I_{X_{i}})\cdot\sigma(\gamma_{ij})\cdot\mathbb{E}_{I_{X_{i}},\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})\right]-\\\
&p(I_{X_{j}})\cdot\sigma(\gamma_{ji})\cdot\mathbb{E}_{I_{X_{j}},\bm{X},C_{-ij}}\left[\mathcal{L}_{X_{j}\to
X_{i}}(X_{i})-\mathcal{L}_{X_{j}\not\to
X_{i}}(X_{i})\right]\Big{)}\end{split}$ (22)

To align the equation with the one in Section 3.3, we swap the indices $k,l$
with $i,j$ again. The first line represents cases where we have an
intervention on the variable $X_{i}$, while we have it over interventions on
the variable $X_{j}$ in the second line. The two terms are weighted based on
the edge existence likelihood $\sigma(\gamma_{ij})$ and $\sigma(\gamma_{ji})$
respectively, and the likelihood of performing an intervention on $X_{i}$ or
$X_{j}$. In our experiments, we use a uniform probability across interventions
on variables, but emphasize that this is not strictly required. Moreover, one
could design heuristics that selects the 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 case.

We clarify that we do not consider the gradients of $\theta_{ij}$ 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 $\theta_{ij}$ if the edge
existence for the opposite direction would be greater than for the direction
from $X_{i}$ to $X_{j}$, _i.e._ $\gamma_{ij}<\gamma_{ji}$, and decrease
$\theta_{ij}$ if we have $\gamma_{ij}>\gamma_{ji}$. However, the orientation
should only model the causal direction of an edge. Hence, we do not gain any
value from a causal perspective when adding the regularizer to the gradient.
Secondly, the regularizer would require us to take additional assumptions for
guaranteeing the discovery of the true graph upon convergence. In experiments
with using a regularizer in the $\theta$-gradient, we did not observe any
difference to the experiments without the regularizer.

We note that the orientation parameters are considered to be pairwise
independent. In other words, $\theta_{ij}$ and $\theta_{kl}$ are considered
independent parameters if $i\neq k,l$ and $j\neq k,l$. 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 experiments. 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(I)$ is set to a uniform distribution for
all our experiments. However, the distribution can also be replaced by a
heuristic which selects interventions to increase computational efficiency. To
keep the convergence guarantees, $p(I)$ would have to guarantee a non-zero
probability to pick any variable. In experiments, we experienced that the Adam
optimizer Kingma & Ba (2015) speeds up the convergence of the parameters
$\bm{\gamma}$ and $\bm{\theta}$ while not interfering with the convergence
guarantees in practice.

Data: $N$ variables $\\{X_{1},...,X_{N}\\}$; observational dataset
$D_{\text{obs}}$; interventional datasets $D_{\text{int}}(\hat{I})$ for
sparse, perfect interventions on all variables; distribution over
interventions $p(I)$

Result: A graph structure corresponding to the cuasal relations among
variables

Initialize $\bm{\gamma}=\bm{0},\bm{\theta}=\bm{0}$

for  _number of epochs_ do

/* Distribution fitting */

for  _$F$ iterations_ do

$\bm{X}\sim D_{\text{obs}}$

for  _$i\leftarrow 1$ to $N$_ do

$M_{l}\sim\text{Ber}(\sigma(\gamma_{li})\cdot\sigma(\theta_{li}))$

$L=-\log f_{\phi_{i}}(X_{i};\bm{M}_{-i}\odot\bm{X}_{-i})$

$\phi_{i}\leftarrow\text{Adam}(\phi_{i},\nabla_{\phi_{i}}L)$

end for

end for

/* Graph fitting */

for  _$G$ iterations_ do

/* Sample an intervention */

$\hat{I}\sim p(I)$ with intervention target $X_{t}$

$\bm{X}\sim D_{\text{int}}(\hat{I})$

/* Evaluate multiple graph samples for gradient estimator */

for  _$k=1$ to $K$_ do

$C^{(k)}\sim\text{Ber}(\sigma(\bm{\gamma})\cdot\sigma(\bm{\theta}))$

for  _$i=1$ to $N$_ do

$\mathcal{L}_{C^{(k)}}(X_{i})\leftarrow-\log
f_{\phi_{i}}(X_{i};C^{(k)}_{:,i}\odot\bm{X}_{-i})$

end for

end for

/* Update parameters */

for  _$i=1$ to $N$_ do

for  _$j=1$ to $N$ where $j\neq i$ and $j\neq t$_ do

/* Considering edge $X_{i}\to X_{j}$ */

$\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})\leftarrow\frac{\sum_{k=1}^{K}C^{(k)}_{ij}\cdot\mathcal{L}_{C^{(k)}}(X_{j})}{\sum_{k=1}^{K}C^{(k)}_{ij}}$

$\mathcal{L}_{X_{i}\not\to
X_{j}}(X_{j})\leftarrow\frac{\sum_{k=1}^{K}(1-C^{(k)}_{ij})\cdot\mathcal{L}_{C^{(k)}}(X_{j})}{\sum_{k=1}^{K}(1-C^{(k)}_{ij})}$

$\gamma_{ij}\leftarrow\gamma_{ij}-\sigma^{\prime}(\gamma_{ij})\cdot\sigma(\theta_{ij})\cdot\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to
X_{j}}(X_{j})+\lambda_{\text{sparse}}\right]$

if  _$i==t$_ then

$\theta_{ij}\leftarrow\theta_{ij}-\sigma^{\prime}(\theta_{ij})\cdot\sigma(\gamma_{ij})\cdot\left[\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})\right]$

$\theta_{ji}\leftarrow-\theta_{ij}$

end if

end for

end for

end for

end for

return $G=(V,E)$ with $V=\bm{X}$ and $E=\\{(X_{i},X_{j})\mid i,j\in[1,N]\text{
and }i\neq j\text{ and }\sigma(\gamma_{ij})>0.5\text{ and
}\sigma(\theta_{ij})>0.5\\}$

Algorithm 1 Learning algorithm of ENCO

##  Appendix 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, it returns the correct causal graph. This is a
stronger statement than showing that the global optimum corresponds to the
true graph, since a gradient-based algorithm can get stuck in a 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
are needed for the guarantee, and then give a sketch of the proof. The proof
will first assume that we work in the data limit, _i.e._ have given sufficient
data, such that we 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

##### Assumption 1

We are given a dataset of observational data from the joint distribution
$p(\bm{X})$. Additionally, we have $N$ interventional datasets for $N$
variables where in each intervention a different node is intervened on (the
intervention size for each dataset is therefore 1).

##### Assumption 2

A common assumption in causal structure learning is that the data distribution
over all variables $p(\bm{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 $\text{pa}(X_{i})$ reflects the
inputs to the causal generation mechanism of $X_{i}$. This allows us to also
handle deterministic variables.

##### Assumption 3

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 4

ENCO relies on neural networks to determine the conditional data distributions
$p(X_{i}|...)$. Hence, for providing a guarantee, we 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_{i}|...)$. 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 $\epsilon$, based on the
universal approximation theorem Hornik et al. (1989). For the limited data
setting, see Appendix B.2.3.

##### 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 $\theta_{ij}$ will
converge to $+\infty$, i.e. $\sigma(\theta_{ij})\to 1$, if $X_{i}$ is an
ancestor of $X_{j}$. Similarly, if $X_{i}$ is a descendant of $X_{j}$, the
parameter $\theta_{ij}$ will converge to $-\infty$, i.e.
$\sigma(\theta_{ij})\to 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, $\gamma_{ij}$ will converge to 1.
Note that we need to take additional assumptions/conditions with respect to
$\lambda_{\text{sparse}}$ here.

Step 3

Once the parameters $\gamma_{ij}$ and $\theta_{ij}$ have converged for the
edges in the ground truth graph, we 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 $\bm{\theta}$ and
$\bm{\gamma}$ 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_{i}\to X_{j}$ in the true causal graph. The orientation
parameter $\theta_{ij}$ converges to $\sigma(\theta_{ij})=1$ if the following
two conditions are fulfilled:

  1. (1)

for all possible sets of parents of $X_{j}$ excluding $X_{i}$, adding $X_{i}$
improves the log-likelihood estimate of $X_{j}$ under the intervention on
$X_{i}$, or leaves it unchanged:

$\forall\widehat{\text{pa}}(X_{j})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{i}},\bm{X}}\left[\log
p(X_{j}|\widehat{\text{pa}}(X_{j}),X_{i})-\log
p(X_{j}|\widehat{\text{pa}}(X_{j}))\right]\geq 0$ (23)

  1. (2)

there exists a set of nodes $\widehat{\text{pa}}(X_{j})$, for which the
probability to be sampled as parents of $X_{j}$ is greater than 0, and the
following condition holds:

$\exists\widehat{\text{pa}}(X_{j})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{i}},\bm{X}}\left[\log
p(X_{j}|\widehat{\text{pa}}(X_{j}),X_{i})-\log
p(X_{j}|\widehat{\text{pa}}(X_{j}))\right]>0$ (24)

###### Proof.

Based on the conditions in Equations 23 and 24, we need to show that the
gradient of $\theta_{ij}$ is negative in expectation, independent of other
values of $\bm{\gamma}$ and $\bm{\theta}$. For readability, we define the
following function:

$T(X_{k},X_{l})=\mathbb{E}_{I_{X_{k}},\bm{X},C_{-kl}}\left[\mathcal{L}_{X_{k}\to
X_{l}}(X_{l})-\mathcal{L}_{X_{k}\not\to X_{l}}(X_{l})\right]$ (25)

Hence, the gradient of $\theta_{ij}$ can be written as:

$\frac{\partial}{\partial\theta_{ij}}\tilde{\mathcal{L}}=\sigma^{\prime}(\theta_{ij})\cdot\Big{(}p(I_{X_{i}})\cdot\sigma(\gamma_{ij})\cdot
T(X_{i},X_{j})-p(I_{X_{j}})\cdot\sigma(\gamma_{ji})\cdot
T(X_{j},X_{i})\Big{)}$ (26)

Looking at the gradient of $\theta_{ij}$ in Equation 26, the conditions
correspond to $T(X_{i},X_{j})$ being smaller or equals to zero. Note that the
sign is flipped because in $T(X_{i},X_{j})$, we have negative log-likelihoods
represented by $\mathcal{L}_{X_{k}\to X_{l}}(X_{l})$, while in Equations 23
and 24, we have log-likelihoods. Further, the other factors of
$\sigma^{\prime}(\theta_{ij})$, $\sigma(\gamma_{ij})$ and $p(I_{X})$ are all
limited in the range of $(0,1)$ meaning that the sign of the gradient is
solely determined by $T(X_{i},X_{j})$ and $T(X_{j},X_{i})$. If
$T(X_{i},X_{j})-T(X_{j},X_{i})$ is smaller than zero, then the gradient of
$\theta_{ij}$ is negative, _i.e._ increasing $\theta_{ij}$.

First, we look at when $T(X_{i},X_{j})<0$. The condition in Equation 23
ensures that conditioning $X_{j}$ on a true parent $X_{i}$ when intervening on
$X_{i}$ 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.2.4). The second
condition, Equation 24, guarantees that there is at least one parent set for
which $T(X_{i},X_{j})$ is negative. Therefore, in expectation over all
possible adjacency matrices $p_{\bm{\gamma},\bm{\theta}}(C)$, $T(X_{i},X_{j})$
is smaller than zero if the two conditions hold.

To guarantee that the whole gradient of $\theta_{ij}$ is negative, we also
need to show that for interventions on $X_{j}$, $T(X_{j},X_{i})$ can only be
positive. When intervening on $X_{j}$, $X_{i}$ and $X_{j}$ become independent
as the edge $X_{i}\to X_{j}$ is removed in the intervened graph. A
distribution $p(X_{i}|X_{j},...)$ relying on correlations between $X_{i}$ and
$X_{j}$ from observational data cannot achieve a better estimate than the same
distribution when removing $X_{j}$. This is because the cross entropy is
minimized when the sampled distribution, in this case $p(X_{i})$, is equal to
the log-likelihood estimator Cover & Thomas (2005):

$-\sum_{x_{i},x_{j}}p(X_{i}=x_{i})\log
p(X_{i}=x_{i})\leq\sum_{x_{i},x_{j}}p(X_{i}=x_{i})\log
q(X_{i}=x_{i}|X_{j}=x_{j})$ (27)

The only situation where $X_{i}$ and $X_{j}$ can become conditionally
dependent under interventions on $X_{j}$ is if $X_{i}$ and $X_{j}$ share a
collider $X_{k}$, and $X_{i}$ is being conditioned on the collider $X_{k}$ and
$X_{j}$. However, this requires that $\theta_{ki}$ has negative gradients,
i.e. $\theta_{ki}$ increasing, when intervening on $X_{k}$. This cannot be the
case since under interventions on $X_{k}$, $X_{i}$ and $X_{k}$ become
conditionally independent, and the correlations learned from observational
data cannot be transferred to the interventional setting. If $X_{k}$ and
$X_{i}$ again share a collider, we can apply this argumentation recursively
until a node $X_{n}$ does not share a collider with $X_{i}$. The recursion
will always come to an end as we have a finite set of nodes, and the causal
graph is assumed to be acyclic.

Thus, if the conditions in Equations 23 and 24 hold for an edge $X_{i}\to
X_{j}$ in the causal graph, we can guarantee that with sufficient time and
data, the corresponding orientation parameter $\theta_{ij}$ will converge to
$\sigma(\theta_{ij})=1$. ∎

######  Theorem B.2.

Consider a pair of variables $X_{i},X_{j}$ for which $X_{i}$ is an ancestor of
$X_{j}$ without direct edge in the true causal graph. Then, the orientation
parameter of the edge $X_{i}\to X_{j}$ converges to $\sigma(\theta_{ij})=1$ if
the same conditions as in Theorem B.1 hold for the pair of $X_{i},X_{j}$.

###### Proof.

To show this theorem, we need to consider two cases for a pair of variables
$X_{i}$ and $X_{j}$: $X_{i}$ and $X_{j}$ are conditionally independent under a
sampled adjacency matrix $C$, or $X_{i}$ and $X_{j}$ are not independent. Both
cases need to be considered for an intervention on $X_{i}$ with the log-
likelihood estimate of $X_{j}$, and an intervention on $X_{j}$ with the log-
likelihood estimate of $X_{i}$.

First, we discuss interventions on $X_{i}$. If under the sampled adjacency
matrix $C$, $X_{j}$ is conditionally independent of $X_{i}$, the difference in
the log-likelihood estimates $T(X_{i},X_{j})$ is zero in expectation. The
variables can be independent if, for example, the parents of $X_{j}$ are all
parents of the true causal graph. If $X_{j}$ is not conditionally independent
of $X_{i}$, the conditions in Equations 23 and 24 from Theorem B.1 ensure that
$X_{i}$ has, in expectation, a positive effect on the log-likelihood estimate
of $X_{j}$. Thus, under interventions on $X_{i}$, the gradient of
$\theta_{ij}$ is smaller or equals to zero, _i.e._ increases $\theta_{ij}$.

Next, we consider interventions on $X_{j}$. If under the sampled adjacency
matrix $X_{i}$ is conditionally independent of $X_{j}$, the difference in the
log-likelihood estimates $T(X_{j},X_{i})$ is zero. The variables can be
independent if $X_{i}$ is conditioned on variables that d-separate $X_{i}$ and
$X_{j}$ in the true causal graph. For instance, having the children of $X_{i}$
as parents of $X_{i}$ creates this scenario. However, for this scenario to
take place, one or more orientation parameters of parent-child or ancestor-
descendant pairs must be incorrectly converged. In case of a parent-child pair
$X_{i},X_{k}$, Theorem B.1 shows that $\sigma(\theta_{ik})$ will converge to
one removing any possibility of a reversed edge to be sampled. In case of an
ancestor-descendant pair $X_{i},X_{l}$, we can apply a recursive argument: as
$X_{l}$ d-separates $X_{i}$ and $X_{j}$, $X_{l}$ must come before $X_{j}$ in
the causal order. If for the gradient $\theta_{il}$, we have a similar
scenario with $X_{i}$ being conditionally independent of $X_{j}$, the same
argument applies. This can be recursively applied until no more variables
except direct children of $X_{i}$ can d-separate $X_{i}$ and $X_{j}$. In that
case, $\sigma(\theta_{ik})$ will converge to one, which leads to all other
orientation parameters to converge to one as well. If $X_{i}$ is not
conditionally independent of $X_{j}$, we can rely back on the argumentation of
Theorem B.1 when we have an edge $X_{i}\to X_{j}$: as in the intervened causal
graph, $X_{i}$ and $X_{j}$ are independent, any correlation learned from
observational data can only lead to a worse log-likelihood estimate. In cases
of colliders, we can rely on the recursive argument from before. Thus, under
interventions on $X_{j}$, the gradient of $\theta_{ij}$ must be smaller or
equals to zero in expectation, _i.e._ increases $\theta_{ij}$.

Therefore, we can conclude that $\sigma(\theta_{ij})$ converges to one for any
ancestor-descendant pairs $X_{i},X_{j}$ under the conditions in Theorem B.1. ∎

######  Theorem B.3.

Consider an edge $X_{i}\to X_{j}$ in the true causal graph. The parameter
$\gamma_{ij}$ converges to $\sigma(\gamma_{ij})=1$ if the following condition
holds:

$\min_{\hat{\text{pa}}\subseteq\text{gpa}_{i}(X_{j})}\mathbb{E}_{\hat{I}\sim
p_{I_{-j}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\bm{X})}\big{[}\log
p(X_{j}|\hat{\text{pa}},X_{i})-\log
p(X_{j}|\hat{\text{pa}})\big{]}>\lambda_{\text{sparse}}$ (28)

where $\text{gpa}_{i}(X_{j})$ is the set of nodes excluding $X_{i}$ which,
according to the ground truth graph, could have an edge to $X_{j}$ without
introducing a cycle, and $p_{I_{-j}}(I)$ refers to the distribution over
interventions $p_{I}(I)$ excluding the intervention on variable $X_{j}$.

###### Proof.

To show this convergence, we assume that the orientation parameters have
converged corresponding to Theorem B.1 and B.2. The parameter $\gamma_{ij}$
converges to $\sigma(\gamma_{ij})=1$ if its gradient,
$\frac{\partial}{\partial\gamma_{ij}}\mathcal{\tilde{L}}$, is negative
independent of other values of $\bm{\gamma}$ and orientation parameters
$\bm{\theta}$ that are not included in Theorem B.1 and B.2. The gradient of
$\gamma_{ij}$ includes an expectation over adjacency matrices
$p_{\bm{\gamma},\bm{\theta}}(C)$. Based on the converged $\bm{\theta}$-values,
we only need to consider sets of nodes as parents for $X_{j}$ that contain
parents, ancestors, or (conditionally) independent nodes according to the
ground truth graph. This sets of parents is represented by
$\text{gpa}_{i}(X_{j})$. Among those remaining parent sets, we need to ensure
that for any such set, the gradient is negative. The condition in Equation 28
corresponds to the inequality
$\frac{\partial}{\partial\gamma_{ij}}\mathcal{\tilde{L}}<0$ since the term on
the left represents the log-likelihood difference $\mathcal{L}_{X_{i}\to
X_{j}}(X_{j})-\mathcal{L}_{X_{i}\not\to X_{j}}(X_{j})$ in the gradients of
$\gamma_{ij}$ in Equation 20 with a flipped sign. For readability and better
interpretation, $\lambda_{\text{sparse}}$ has been moved on the right site of
the inequality. This is possible as $\lambda_{\text{sparse}}$ is independent
of the two expectations in Equation 28. If the inequality holds for all parent
sets $\hat{\text{pa}})$, the gradient of $\gamma_{ij}$ can be guaranteed to be
negative in expectation, independent of the other values of $\bm{\gamma}$.
Since the distribution over parent sets $\hat{\text{pa}})$ depends on other
values of $\bm{\gamma}$, the condition in Equation 28 ensures that even for
the parent set with the lowest log-likelihood difference, it is still larger
than $\lambda_{\text{sparse}}$. If this condition holds, then the gradient of
$\frac{\partial}{\partial\gamma_{ij}}\mathcal{\tilde{L}}$ will be smaller than
zero independent of other values of $\bm{\gamma}$. ∎

The condition in Equation 28 introduces a dependency between convergence
guarantees and the regularizer parameter $\lambda_{\text{sparse}}$. The lower
we set the regularization weight $\lambda_{\text{sparse}}$, the more edges we
can guarantee to recover. If the regularization weight is set too high, we can
eventually obtain 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 a limited
data regime. Nonetheless, if sufficient computational resources and data is
provided, any value of $\lambda_{\text{sparse}}>0$ can be used.

######  Theorem B.4.

Assume for all edges $X_{i}\to X_{j}$ in the true causal graph,
$\sigma(\theta_{ij})$ and $\sigma(\gamma_{ij})$ have converged to one. Then,
the likelihood of all other edges, _i.e._
$\sigma(\theta_{lk})\cdot\sigma(\theta_{lk})$, will converge to zero under the
condition that $\lambda_{\text{sparse}}>0$.

###### Proof.

If all edges in the ground truth graph have converged, all other pairs of
variables $X_{l},X_{k}$ are (conditionally) independent in the graph. This
statement follows from the Markov property of the graph and excludes ancestor-
descendant pairs $X_{i},X_{j}$. The possibility of having edges from
descendants to ancestors has been removed by the fact that the orientation
parameters $\theta_{ij}$ have converged according to Theorem B.2. Thus, for
those cases, we already have the guarantee that
$\sigma(\theta_{ij})\cdot\sigma(\theta_{ij})$ converges to zero.

For a conditionally independent pair $X_{l},X_{k}$, the difference of the log-
likelihood estimate in the gradient of $\gamma_{lk}$, _i.e._
$\mathcal{L}_{X_{l}\to X_{k}}(X_{k})-\mathcal{L}_{X_{l}\not\to X_{k}}(X_{k})$,
is zero in expectation since independent nodes do not share any information.
Thus, the gradient remaining is:

$\frac{\partial}{\partial\gamma_{lk}}\tilde{\mathcal{L}}=\sigma^{\prime}(\gamma_{lk})\cdot\sigma(\theta_{lk})\cdot\lambda_{\text{sparse}}$
(29)

Since the gradient is positive independent of the values of $\gamma_{lk}$ and
$\theta_{lk}$, $\gamma_{lk}$ will decrease until it converges to
$\sigma(\gamma_{lk})=0$.

Hence, if $\gamma_{lk}$ decreases for all pairs of (conditionally) independent
variables $X_{l},X_{k}$ in the ground truth graph, and $\sigma(\theta_{lk})$
converged to zero for children and descendants, the product
$\sigma(\gamma_{lk})\cdot\sigma(\theta_{lk})$ will converge to zero for all
edges not existing in the ground truth graph. ∎

For graphs that fulfill all conditions in the Theorems B.1 and B.4, ENCO is
guaranteed to converge given sufficient data and time. The conditions in 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 $\bm{\gamma}$ and
$\bm{\theta}$.

##### Summary

We can summarize the conditions discussed above as follows. Given a causal
graph $\mathcal{G}$ with variables $X_{1},...,X_{N}$ and sparse interventions
on all variables, the proposed method ENCO will converge to the true, causal
graph $\mathcal{G}$, if the following three conditions hold for all edges
$X_{i}\to X_{j}$ in the true causal graph $\mathcal{G}$:

  1. 1.

For all possible sets of parents of $X_{j}$ excluding $X_{i}$, adding $X_{i}$
improves the log-likelihood estimate of $X_{j}$ under the intervention on
$X_{i}$, or leaves it unchanged:

$\forall\widehat{\text{pa}}(X_{j})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{i}},\mathbf{X}}\left[\log
p(X_{j}|\widehat{\text{pa}}(X_{j}),X_{i})-\log
p(X_{j}|\widehat{\text{pa}}(X_{j}))\right]\geq 0$ (30)

  2. 2.

There exists a set of nodes $\widehat{\text{pa}}(X_{j})$, for which the
probability to be sampled as parents of $X_{j}$ is greater than 0, and the
following condition holds:

$\exists\widehat{\text{pa}}(X_{j})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{i}},\mathbf{X}}\left[\log
p(X_{j}|\widehat{\text{pa}}(X_{j}),X_{i})-\log
p(X_{j}|\widehat{\text{pa}}(X_{j}))\right]>0$ (31)

  3. 3.

The effect of $X_{i}$ on $X_{j}$ cannot be described by other variables up to
$\lambda_{\text{sparse}}$:

$\min_{\hat{\text{pa}}\subseteq\text{gpa}_{i}(X_{j})}\mathbb{E}_{\hat{I}\sim
p_{I_{-j}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{j}|\hat{\text{pa}},X_{i})-\log
p(X_{j}|\hat{\text{pa}})\big{]}>\lambda_{\text{sparse}}$ (32)

where $\text{gpa}_{i}(X_{j})$ is the set of nodes excluding $X_{i}$ which,
according to the ground truth graph, could have an edge to $X_{j}$ without
introducing a cycle.

Further, for all other pairs $X_{i},X_{j}$ for which $X_{j}$ is a descendant
of $X_{i}$, 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 on the precise
calculations, we provide a Jupyter Notebook that contains all calculations in
this example111The calculations can be found in the notebook called
convergence_guarantees_ENCO.ipynb, see
https://github.com/phlippe/ENCO/blob/main/convergence_guarantees_ENCO.ipynb ..

Suppose we have a graph with 3 binary variables, $X_{1}$, $X_{2}$, $X_{3}$,
with the causal graph being $X_{1}\to X_{2}\to X_{3}$, i.e., a small chain.
For simplicity, let us assume that the true, conditional distributions are the
following:

$\displaystyle p(X_{1})$ $\displaystyle=$ $\displaystyle\text{Bern}(0.7)$ (33)
$\displaystyle p(X_{2}|X_{1})$ $\displaystyle=$
$\displaystyle\begin{cases}X_{1}&\text{with prob. }0.6\\\ X_{1}\oplus
1&\text{with prob. }0.4\end{cases}$ (34) $\displaystyle p(X_{3}|X_{2})$
$\displaystyle=$ $\displaystyle\begin{cases}X_{2}&\text{with prob. }0.2\\\
X_{2}\oplus 1&\text{with prob. }0.8\end{cases}$ (35)

In other words, $X_{2}$ is equals to the value of $X_{1}$ with a probability
of $0.6$, and the opposite binary value otherwise. Similarly, $X_{3}$ is
equals to the value of $X_{2}$ with a 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
$X_{1}=1,X_{2}=1,X_{3}=0$. Further, we assume that all interventions replace
the respective conditional distribution by a uniform distribution, i.e.,
$p_{I_{X_{i}}}(X_{i})=\text{Bern}(0.5)$. Next, we will check the conditions
for the edges in $\mathcal{G}$, i.e., $X_{1}\to X_{2}$ and $X_{2}\to X_{3}$,
and the remaining ancestor-descendant pair $X_{1},X_{3}$.

$X_{1}$$X_{2}$$X_{3}$ 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.

Edge $X_{1}\to X_{2}$:

  * •

Condition 1: the possible parent sets that exclude $X_{1}$ and $X_{2}$ are
$X_{3}$ and the empty set. For the empty set, we get:

$\mathbb{E}_{I_{X_{1}},\mathbf{X}}\left[\log p(X_{2}|X_{1})-\log
p(X_{2})\right]\approx 0.023\geq 0$

For conditioning on $X_{3}$, we obtain:

$\mathbb{E}_{I_{X_{1}},\mathbf{X}}\left[\log p(X_{2}|X_{1},X_{3})-\log
p(X_{2}|X_{3})\right]\approx 0.015\geq 0$

Since both values are greater than zero, condition 1 is fulfilled for
$X_{1}\to X_{2}$.

  * •

Condition 2: is already fulfilled by the equations in condition 1 since all
parent sets have a difference greater than zero.

  * •

Condition 3: the set $\text{gpa}_{1}(X_{2})$ is the empty set since $X_{3}$ is
a descendant of $X_{2}$, and no other nodes exist in the graph. Thus, the
parent set minimizing the expression on the left can only be the empty set,
and we can calculate it as follows:

$\mathbb{E}_{\hat{I}\sim
p_{I_{-2}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{2}|X_{1})-\log p(X_{2})\big{]}\approx\underbrace{1/2\cdot
0.023}_{I_{X_{1}}}+\underbrace{1/2\cdot
0.017}_{I_{X_{3}}}=0.020>\lambda_{\text{sparse}}$

with assuming $p_{I}(I)$ being the uniform distribution, and excluding
$I_{X_{2}}$ since we do not update $\gamma_{12}$ in this case. Hence, as long
as $\lambda_{\text{sparse}}$ is smaller than $0.02$, the condition is
fulfilled.

Edge $X_{2}\to X_{3}$:

  * •

Condition 1: the possible parent sets that exclude $X_{2}$ and $X_{3}$ are
$X_{1}$ and the empty set. For the empty set, we get:

$\mathbb{E}_{I_{X_{2}},\mathbf{X}}\left[\log p(X_{3}|X_{2})-\log
p(X_{3})\right]\approx 0.194\geq 0$

For conditioning on $X_{1}$, we obtain:

$\mathbb{E}_{I_{X_{2}},\mathbf{X}}\left[\log p(X_{3}|X_{2},X_{1})-\log
p(X_{3}|X_{1})\right]\approx 0.200\geq 0$

Since both values are greater than zero, condition 1 is fulfilled for
$X_{2}\to X_{3}$.

  * •

Condition 2: is already fulfilled by the equations in condition 1 since all
parent sets have a difference greater than zero.

  * •

Condition 3: the set $\text{gpa}_{2}(X_{3})$ contains the variable $X_{1}$
since we can introduce an edge $X_{1}\to X_{3}$ without introducing acyclicity
in the true, causal graph. Thus, we need to compare two parent sets for
finding the minimum of the left-side term: $X_{1}$ and the empty set. First,
we consider the empty set:

$\mathbb{E}_{\hat{I}\sim
p_{I_{-3}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{3}|X_{2})-\log p(X_{3})\big{]}\approx\underbrace{1/2\cdot
0.194}_{I_{X_{1}}}+\underbrace{1/2\cdot
0.194}_{I_{X_{2}}}=0.194>\lambda_{\text{sparse}}$

Again, we exclude $I_{X_{3}}$ since we do not update $\gamma_{23}$ in this
case. The second case considers $X_{1}$ as additional parent set
$\hat{\text{pa}}$:

$\begin{split}\mathbb{E}_{\hat{I}\sim
p_{I_{-3}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{3}|X_{2},X_{1})-\log p(X_{3}|X_{1})\big{]}&\approx\underbrace{1/2\cdot
0.186}_{I_{X_{1}}}+\underbrace{1/2\cdot 0.200}_{I_{X_{2}}}\\\
&=0.193>\lambda_{\text{sparse}}\\\ \end{split}$

The minimum of both values is $0.193$. Hence, the edge $X_{2}\to X_{3}$ can be
recovered if $0.193>\lambda_{\text{sparse}}$.

Ancestor-descendant pair $X_{1},X_{3}$:

  * •

Condition 1: the possible parent sets that exclude $X_{1}$ and $X_{3}$ are
$X_{2}$ and the empty set. For the empty set, we get:

$\mathbb{E}_{I_{X_{1}},\mathbf{X}}\left[\log p(X_{3}|X_{1})-\log
p(X_{3})\right]\approx 0.008\geq 0$

For conditioning on $X_{2}$, we obtain:

$\mathbb{E}_{I_{X_{1}},\mathbf{X}}\left[\log p(X_{3}|X_{1},X_{2})-\log
p(X_{3}|X_{2})\right]=0\geq 0$

The difference is zero because $X_{3}$ is independent of $X_{1}$ when
conditioned on $X_{2}$: $p(X_{3}|X_{1},X_{2})=p(X_{3}|X_{2})$ Since both
values are greater or equals to zero, condition 1 is fulfilled for the pair
$X_{1},X_{3}$.

  * •

Condition 2: from condition 1, we can see that the parent set of
$\widehat{\text{pa}}(X_{3})$ being the empty set is the only option that
fulfills the condition being greater than zero. Since we start the
optimization process with an initialization that assigns a non-zero
probability to all possible parent sets, it follows that
$\widehat{\text{pa}}(X_{3})$ being the empty set has a probability greater
than zero throughout the optimization process. Hence, condition 2 is fulfilled
as well.

Summary: in conclusion, for the discussed example, we can guarantee that ENCO
converges to the correct causal graph if $\lambda_{\text{sparse}}<0.02$. To
experimentally verify this results, we applied ENCO on this graph with two
hyperparameter settings for the sparsity regularizer:
$\lambda_{\text{sparse}}=0.019$ and $\lambda_{\text{sparse}}=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
$\lambda_{\text{sparse}}=0.019$, ENCO was able to recover the graph without
errors while for $\lambda_{\text{sparse}}=0.021$, the edge $X_{1}\to X_{2}$
was, as expected, missed. This verifies the theoretical result above with
respect to $\lambda_{\text{sparse}}$. Note that if the condition is not
fulfilled by selecting a 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 1 and 2

As mentioned in the text, condition 1 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 a 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 a linear Gaussian
setting when the ancestor has a 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 1 and 2, it 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_{i}\to X_{j}$, there exist one or more parent(s) of $X_{j}$ sharing a
common confounder $X_{k}$ with $X_{i}$, where $X_{k}$ is not a direct parent
of $X_{j}$. An example of this structure is the following: $X_{1}\to
X_{2},X_{3}$; $X_{2},X_{3}\to X_{4}$ where the orientations of the edges
$X_{2}\to X_{4}$ and $X_{3}\to X_{4}$ 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 $X_{1},X_{2},X_{3}$ with
$X_{1}\to X_{2}$, $X_{3}\to X_{2}$, and $X_{2}$ having no parents besides
$X_{1}$ and $X_{3}$. If $X_{1}$ and $X_{3}$ do not share a confounder, then,
from do-calculus, we know that
$p(X_{2}|\text{do}(X_{1}=x_{1}))=p(X_{2}|X_{1}=x_{1})$ and
$p(X_{2}|\text{do}(X_{3}=x_{3}))=p(X_{2}|X_{3}=x_{3})$. Furthermore, since the
conditional entropy of a variable can only be smaller or equals to the
marginal, i.e. $H(X)\geq H(X|Y)$, estimating $X_{2}$ under interventions on
$X_{1}$ can only be improved by conditioning on $X_{1}$, and similarly for
$X_{3}$. Thus, condition 1 is 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 $X_{1}$ and $X_{3}$ share a common
confounder. Then, from do-calculus, we can state that there can exist a
parameterization of the conditional distributions for which
$p(X_{2}|\text{do}(X_{1}=x_{1}))\neq p(X_{2}|X_{1}=x_{1})$. Under this
setting, we cannot guarantee that condition 1 is always fulfilled. However,
whether this parent-confounder structure above actually leads to a local
optimum or 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 to a graph
with the correct orientation.

Besides, if the confounder $X_{j}$ of $X_{1}$ and $X_{3}$ is a parent of
$X_{2}$, then the local optimum would disappear with learning that edge since
$p(X_{2}|\text{do}(X_{1}=x_{1}),X_{j})=p(X_{2}|X_{1}=x_{1},X_{j})$. 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, we often have a limited data set. The
proof presented for the data limit is straightforward to extend to this
setting with the following modification:

  * •

The conditional distributions $p(X|...)$ are replaced by the conditional
distributions that follow from the given, observational data.

  * •

The expectations over the interventional data $\tilde{p}_{\hat{I}}(\bm{X})$ is
replaced by the joint distribution over samples given for the intervention
$\hat{I}$.

  * •

Theorem B.1 and B.2 for the edge $X_{i}\to X_{j}$ are extended as follows:

    1. (3)

For all possible sets of parents of $X_{i}$ excluding $X_{j}$, adding $X_{j}$
does not improve the log-likelihood estimate of $X_{i}$ under the intervention
on $X_{j}$, or leaves it unchanged:

$\forall\widehat{\text{pa}}(X_{i})\subseteq
X_{-i,j}:\mathbb{E}_{I_{X_{j}},\bm{X}}\left[\log
p(X_{i}|\widehat{\text{pa}}(X_{i}),X_{j})-\log
p(X_{i}|\widehat{\text{pa}}(X_{i}))\right]\leq 0$ (36)

This condition is the inverse statement of Equation 23, in the sense that we
consider interventions on the child/descendant $X_{j}$. 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 an overestimation of edges. Thus, we add the following condition:

    1. (1)

For all pairs of variables $X_{i},X_{j}$ for which there exists no direct
causal relation in the true causal graph, and $X_{j}$ not being the ancestor
of $X_{i}$, the following condition has to hold:

$\min_{\hat{\text{pa}}\subseteq(\text{gpa}_{i}(X_{j})\setminus\text{pa}(X_{j}))}\mathbb{E}_{\hat{I}\sim
p_{I_{-j}}(I)}\mathbb{E}_{\tilde{p}_{\hat{I}}(\mathbf{X})}\big{[}\log
p(X_{j}|\text{pa}(X_{j}),\hat{\text{pa}},X_{i})-\log
p(X_{j}|\text{pa}(X_{j}),\hat{\text{pa}})\big{]}<\lambda_{\text{sparse}}$ (37)

where $\text{gpa}_{i}(X_{j})$ is the set of nodes excluding $X_{i}$ which,
according to the ground truth graph, could have an edge to $X_{j}$ without
introducing a 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 to 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 guarantees

Most common causal graphs fulfill the conditions mentioned above, as long as a
small enough value for $\lambda_{\text{sparse}}$ is chosen. Still, there are
situations where we cannot guarantee that ENCO convergences to the correct
causal graph independent of the chosen value of $\lambda_{\text{sparse}}$.
Here, we want to discuss two scenarios visualized in Figure 8 under which the
guarantees fail. Still, we want to 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
$\bm{\theta}$ and $\bm{\gamma}$ in all situations.

$X_{3}$$X_{1}$$X_{2}$ (a) Fork example

$X_{1}$$X_{2}$$X_{3}$ (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 first example we discuss is based on a fork structure where we have three
binary variables, $\\{X_{1},X_{2},X_{3}\\}$, and the edges $X_{1}\to X_{3}$
and $X_{2}\to X_{3}$ (see Figure 8(a)). The parameterization we look at is a
(noisy) XOR-gate for $X_{3}$ with its two input variables $X_{1},X_{2}$ being
independent of each other and uniformly distributed. The conditional
probability distribution $p(X_{3}|X_{1},X_{2})$ can be summarized in the
following probability function:

$p(X_{3}=1|X_{1},X_{2})=\begin{cases}\epsilon&\text{if }X_{1}=0,X_{2}=0\\\
1-\epsilon&\text{if }X_{1}=1,X_{2}=0\\\ \epsilon&\text{if }X_{1}=0,X_{2}=1\\\
1-\epsilon&\text{if }X_{1}=1,X_{2}=1\\\ \end{cases}$ (38)

In other words, if $X_{1}\neq X_{2}$, $X_{3}$ is equals $1$ with a likelihood
of $1-\epsilon$. If $X_{1}=X_{2}$, $X_{3}$ is equals $1$ with a likelihood of
$\epsilon$. The issue that this probability table creates is the following.
Knowing only one out of the two variables does not improve the log likelihood
estimate for the output. This is because $X_{1}$ and $X_{2}$ are independent
of each other, and $p(X_{3}|X_{1})=p(X_{3})$ 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
$\lambda_{\text{sparse}}$ 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, we will sample adjacency matrices with
both input variables being a parent of $X_{3}$ with a non-zero probability.
Hence, the log-likelihood difference for $X_{1}$ and $X_{2}$ to $X_{3}$ 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 $\epsilon=\\{0.1,0.2,0.3,0.4,0.45\\}$ using a sparsity
regularizer of $\lambda_{\text{sparse}}=1$e-$4$, and in all cases, ENCO
converged to the correct, acyclic graph. Note that values close to 0.5 for
$\epsilon$ 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.1, more specifically Equation 23. The graph we consider
is a fully connected graph with three variables $X_{1},X_{2},X_{3}$ (see
Figure 8(b)). The scenario can be described as follows: if knowing $X_{2}$
informs the log-likelihood estimate of $X_{3}$ more about $X_{1}$ than about
$X_{2}$ itself, an intervention on $X_{2}$ and a sampled graph with the edge
$X_{2}\to X_{3}$ could lead to a worse likelihood estimate of $X_{3}$ than
without the edge. For this scenario to happen, $p(X_{2}|X_{1})$ must be close
to deterministic. Additionally, $p(X_{3}|X_{1},X_{2})$ must be much less
reliant on $X_{2}$ than on $X_{1}$, such as in the following probability
density:

$p(X_{3}=1|X_{1},X_{2})=\begin{cases}\epsilon_{1}&\text{if }X_{1}=0,X_{2}=0\\\
1-\epsilon_{1}&\text{if }X_{1}=1,X_{2}=0\\\ \epsilon_{2}&\text{if
}X_{1}=0,X_{2}=1\\\ 1-\epsilon_{2}&\text{if }X_{1}=1,X_{2}=1\\\ \end{cases}$
(39)

The two variables $\epsilon_{1},\epsilon_{2}$ represent small constants close
to zero. In this case, the graph can violate the condition in Equation 23
since intervening on $X_{2}$ breaks the dependency between $X_{1}$ and
$X_{2}$. The conditional distribution $p(X_{3}|X_{2})$ learned from
observational data relies on the dependency between $X_{1}$ and $X_{2}$ which
can make it to a worse estimate than $p(X_{3})$. Note that if the edge
$X_{1}\to X_{3}$ is learned by ENCO though, this will not constitute a problem
anymore since with conditioning on $X_{1}$, i.e. $p(X_{3}|X_{2},X_{1})$, the
edge $X_{2}\to X_{3}$ will gain a gradient towards the correct graph. Thus,
when $\bm{\gamma}$ and $\bm{\theta}$ are not initialized with the worst-case
values, the graph with both $X_{1}$ and $X_{2}$ as parents of $X_{3}$ 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 $\bm{\gamma}$
and $\bm{\theta}$ 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 to 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 a similar structure to those as 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.

######  Theorem B.5.

For every pair of variables $X_{i},X_{j}$ where $X_{i}$ is a parent of
$X_{j}$, the graph $\hat{G}$ that optimizes objective in Equation 2 models the
orientation $X_{i}\to X_{j}$, if there exists an edge between $X_{i}$ and
$X_{j}$ in $\hat{G}$, under the following conditions:

  * •

$X_{i}$ and $X_{j}$ are not independent under observational data.

  * •

Under interventions on $X_{i}$, $X_{i}$ and $X_{j}$ are not independent given
the true parent set of $X_{j}$.

###### Proof.

If $X_{i}$ and $X_{j}$ 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_{i}|...)=p(X_{i}|X_{j},...)$, thus making the objective invariant to
the orientation of the edge, and removing any edge between $X_{i}$ and $X_{j}$
for sparsity.

If $X_{i}$ and $X_{j}$ are dependent, we can prove the statement by showing
that modeling the orientation $X_{j}\to X_{i}$ will strictly lead to a worse
estimate under the intervention on $X_{j}$ since the orientation parameters
are optimized by comparing the interventions of the two adjacent variables.
Under interventions on $X_{i}$, the causal mechanism
$p(X_{j}|\text{pa}(X_{j}))$, with $\text{pa}(X_{j})$ being the parent set of
the ground truth graph including $X_{i}$, remains invariant under
interventions on $X_{i}$, and is strictly better than
$p(X_{j}|\text{pa}(X_{j})\setminus X_{i})$ for estimating $X_{j}$ due to the
direct causal relation. Under interventions on $X_{j}$, the causal mechanism
$p(X_{i}|X_{j},...)$ leads to a strictly worse estimate as discussed in
Theorem B.1, since the dependency between $X_{i}$ and $X_{j}$ does not exist
in the interventional regime. Hence, the inverse orientation of the edge
$X_{i}\to X_{j}$, i.e. $X_{j}\to X_{i}$ cannot be part of the global optimum.
∎

######  Theorem B.6.

For every pair of variables $X_{i},X_{j}$ where $X_{i}$ is an ancestor but not