Efficient Neural Causal Discovery without Acyclicity Constraints

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 Efficient Neural Causal Discovery.

We make the following four contributions. Firstly, we propose ENCO, a causal structure learning method for observational and interventional data using continuous optimization. Different from recent methods, ENCO models the edge orientation as a separate parameter. Secondly, we derive unbiased, low-variance gradient estimators, which is crucial for scaling up the model to large numbers of variables. Thirdly, we show that ENCO is guaranteed to converge to the correct causal graph if interventions on all variables are available, despite not having any acyclicity constraints. Yet, we show in practice that the algorithm works on partial intervention sets 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.

Background and Related Work

2 Causal structure learning

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.

Efficient Neural Causal Discovery

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.

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

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).

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

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 Low-variance gradient estimators for edge parameters

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:

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:

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.

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:

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}$:

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}$:

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.

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

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.

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:

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.

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.

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.

2 Causal structure learning on common graph structures

We first experiment on synthetic graphs. We pick six common graph structures and sample 5,000 observational data points and 200 per intervention. The graphs 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.

3 Scalability

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 Interventions on fewer variables

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.

5 Detecting latent confounders

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.

6 Real-world inspired data

Finally, we evaluate ENCO on causal graphs from the Bayesian Network Repository (BnLearn) Scutari (2010). The repository contains graphs inspired by real-world applications that are used as benchmarks 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.

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.

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

Appendix A Gradient estimators

$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;

$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.

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

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:

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:

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:

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

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:

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}$:

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

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

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.

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:

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:

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.

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.

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).

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.

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.

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.

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:

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$.

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.

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.

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:

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:

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:

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:

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

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):

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$. ∎

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}$.

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. ∎

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:

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}$.

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.

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$.

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:

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}$.

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}$:

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:

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:

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

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 exampleThe 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:

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}$.

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:

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:

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.

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:

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:

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}}$:

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}}$.

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:

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.

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

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:

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:

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:

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.

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:

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:

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.

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}$.

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. ∎

For every pair of variables $X_{i},X_{j}$ where $X_{i}$ is an ancestor but not direct parents of $X_{j}$, the graph $\hat{G}$ that optimizes objective in Equation 2 does not include the edge $X_{j}\to X_{i}$ if the conditions in Theorem B.5 hold.

To show this statement, we need to consider different independence relations between $X_{i}$ and $X_{j}$. First, if $X_{i}$ and $X_{j}$ are independent in the observational dataset given any conditional set, the edge will be removed since any edge between two independent variables is removed for any $\lambda_{\text{sparse}}>0$. The same holds if $X_{i}$ and $X_{j}$ are independent for interventions on $X_{i}$ and $X_{j}$.

If they are dependent, we can follow a similar argument as in Theorem B.2. The causal mechanism $p(X_{j}|X_{i},...)$ transfers from observational to interventional data on $X_{i}$ since on interventions on $X_{i}$, the causal mechanism of $X_{j}$ is not changed. Further, when intervening on $X_{j}$, $X_{i}$ and $X_{j}$ become independent such that any mechanism $p(X_{i}|X_{j},...)$ cannot transfer except if $X_{i}$ and $X_{j}$ are independent under interventions on $X_{j}$. In this case, the edge will be again removed by the sparsity regularizer. This shows that for any setting, the orientation of the edge $X_{j}\to X_{i}$ cannot lead to a better estimate than $X_{i}\to X_{j}$, and in case of independence, the edge $X_{j}\to X_{i}$ will be removed as well. ∎

The graph $\hat{G}$ that optimizes objective in Equation 2 models the same parent sets for each variable under the following conditions:

For any variable $X_{i}$ with its true parent set $\text{pa}(X_{i})$, there does not exist a smaller parent set $\hat{\text{pa}}\subset\bm{X}\setminus X_{i},\text{descendants}(X_{i})$ which approximates the log-likelihood of $X_{i}$ up to $\lambda_{\text{sparse}}\cdot\left(|\text{pa}(X_{i})|-|\hat{\text{pa}}|\right)$ on average.

The regularization parameter $\lambda_{\text{sparse}}$ is greater than zero.

If the orientations for all edges are according to the ground truth graph in the global optimum following Theorem B.5 and B.6, the parent set for a variable $X_{i}$ is limited to those variables which are not descendants of $X_{i}$. From the ground truth graph, we know that, conditioned on the true parent set $\text{pa}(X_{i})$, $X_{i}$ is independent of all other non-descendant variables. Thus, the log-likelihood estimate of $X_{i}$, i.e. the left part of the objective in Equation 2, is optimized by $p(X_{i}|\text{pa}(X_{i}))$. To show that this is also the global optimum when combining with the regularizer, we need to consider those parent sets of $X_{i}$ which obtain a lower penalty, i.e. smaller parent sets. The difference between two parent sets $\text{pa}(X_{i})$ and $\hat{\text{pa}}$ in terms of the regularizer corresponds to $\lambda_{\text{sparse}}\cdot\left(|\text{pa}(X_{i})|-|\hat{\text{pa}}|\right)$. Thus, if there exists no parent set for which this difference is greater than the penalty for the worse log-likelihood estimate, the true parent set $\text{pa}(X_{i})$ constitutes the global optimum. ∎

B.3 Convergence conditions for latent confounder detection

In Section 3.5, we have discussed that ENCO can be extended to graph with latent confounders. For this, we have to record the gradients of $\gamma_{ij}$ for the interventional data on $X_{i}$ and all other interventional data separately. 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. The score to detect latent confounders is: