# Scalable Differentiable Causal Discovery in the Presence of Latent
Confounders with Skeleton Posterior (Extended Version)

Pingchuan Ma<EMAIL_ADDRESS>HKUSTHong Kong SAR ,  Rui Ding
<EMAIL_ADDRESS>Microsoft ResearchBeijing, China ,  Qiang Fu
<EMAIL_ADDRESS>Microsoft ResearchBeijing, China ,  Jiaru Zhang
<EMAIL_ADDRESS>Shanghai Jiao Tong UniversityShanghai, China ,  Shuai
Wang<EMAIL_ADDRESS>HKUSTHong Kong SAR ,  Shi Han<EMAIL_ADDRESS>Microsoft ResearchBeijing, China and  Dongmei Zhang<EMAIL_ADDRESS>Microsoft ResearchBeijing, China

(2024)

###### Abstract.

Differentiable causal discovery has made significant advancements in the
learning of directed acyclic graphs. However, its application to real-world
datasets remains restricted due to the ubiquity of latent confounders and the
requirement to learn maximal ancestral graphs (MAGs). To date, existing
differentiable MAG learning algorithms have been limited to small datasets and
failed to scale to larger ones (e.g., with more than 50 variables).

The key insight in this paper is that the causal skeleton, which is the
undirected version of the causal graph, has potential for improving accuracy
and reducing the search space of the optimization procedure, thereby enhancing
the performance of differentiable causal discovery. Therefore, we seek to
address a two-fold challenge to harness the potential of the causal skeleton
for differentiable causal discovery in the presence of latent confounders: (1)
scalable and accurate estimation of skeleton and (2) universal integration of
skeleton estimation with differentiable causal discovery.

To this end, we propose SPOT (Skeleton Posterior-guided OpTimization), a two-
phase framework that harnesses skeleton posterior for differentiable causal
discovery in the presence of latent confounders. On the contrary to a “point-
estimation”, SPOT seeks to estimate the posterior distribution of skeletons
given the dataset. It first formulates the posterior inference as an instance
of amortized inference problem and concretizes it with a supervised causal
learning (SCL)-enabled solution to estimate the skeleton posterior. To
incorporate the skeleton posterior with differentiable causal discovery, SPOT
then features a skeleton posterior-guided stochastic optimization procedure to
guide the optimization of MAGs.

Extensive experiments on various datasets show that SPOT substantially
outperforms SOTA methods for MAG learning. SPOT also demonstrates its
effectiveness in the accuracy of skeleton posterior estimation in comparison
with non-parametric bootstrap-based, or more recently, variational inference-
based methods. Finally, we observe that the adoption of skeleton posterior
exhibits strong promise in various causal discovery tasks.

††journalyear: 2024††copyright: acmlicensed††conference: Proceedings of the
30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining; August
25–29, 2024; Barcelona, Spain††booktitle: Proceedings of the 30th ACM SIGKDD
Conference on Knowledge Discovery and Data Mining (KDD ’24), August 25–29,
2024, Barcelona, Spain††doi: 10.1145/3637528.3672031††isbn:
979-8-4007-0490-1/24/08††ccs: Mathematics of computing Bayesian networks††ccs:
Mathematics of computing Causal networks††ccs: Computing methodologies Causal
reasoning and diagnostics

##  1\. Introduction

Causal discovery in the presence of latent confounders is a long-standing
problem (Spirtes et al., 2000; Zhang, 2008). Under this setting, the causal
relations are typically represented by a maximal ancestral graph (MAG)
(Spirtes et al., 2000), a special class of acyclic directed mixed graphs
(ADMGs). MAG learning has historically been conducted by either constraint-
based methods, such as FCI (Spirtes et al., 2000; Zhang, 2008), RFCI (Colombo
et al., 2012) and ICD (Rohekar et al., 2021), or by score-based methods, such
as M3HC (Tsirlis et al., 2018), AGIP (Chen et al., 2021) and GPS (Claassen and
Bucur, 2022). In recent years, differentiable causal discovery has emerged as
a promising approach to enhance the accuracy and efficiency of existing
methods (Zheng et al., 2018; Vowels et al., 2022). By recasting the
combinatorial constraint of graph structure into a differentiable form,
continuous optimization techniques can be applied in an “out-of-the-box”
manner.

The goal of differentiable methods is to identify the ancestral ADMG, and
apply maximal ancestral projection, a standard process, to generate the
corresponding MAG. Despite the encouraging results achieved thus far, they
struggle with large-scale causal graphs, particularly those containing more
than 20 variables. One widely-used method, ABIC (Bhattacharya et al., 2021),
acknowledges its limitations with datasets of only 10–15 variables (i.e., a
causal graph with 10-15 nodes), which restricts its broader applicability.
Likewise, N-ADMG (Ashman et al., 2022), another differentiable method for ADMG
learning, suffers from similar scalability issues. Upon closer examination, we
find that these methods are both inefficient and inaccurate when working with
large-scale causal graphs. For example, ABIC may take days to converge and
sometimes produces cyclic graphs when dealing with 50 nodes. This is likely
due to the inherent challenges of learning from causally insufficient data, as
ABIC must manage a more complicated objective function, additional
optimization variables, and complex constraints. Consequently, the search
space of large datasets poses significant challenges for standard optimization
techniques, such as L-BFGS (Fletcher, 1987) or gradient descent.

On the other hand, it has been demonstrated that adjacency information (i.e.,
the skeleton) can be utilized as a pre/post-processing step to facilitate
differentiable methods for learning DAGs (directed acyclic graphs). For
instance, CDRL (Zhu et al., 2019) uses the CAM (Bühlmann et al., 2014) to rule
out spurious edges while ML4S (Ma et al., 2022a) investigates using skeleton
in a pre-processing step to preclude superfluous variables for NOTEARS (Zheng
et al., 2018). These approaches show the potential for boosting differentiable
DAG learning by leveraging skeleton information.

Key Observation. These promising results on DAG learning motivate us to
investigate whether MAG learning can benefit from a high-quality estimation of
the skeleton. To this end, we implement two ABIC variants, one using the
ground-truth skeleton as prior knowledge and the other using the skeleton
learnt by FCI, to validate our hypotheses. In both variants, we black out the
edges that are not in the skeleton and then run ABIC to only optimize the
remaining edges. In this way, the optimization procedure is only applied on a
subset of variables rather than the entire adjacency matrix. Then, we evaluate
them on ten synthetic datasets with 50-100 variables and 1000 samples (see
Table 1). In particular, we find that ABIC is impressively accurate when the
true skeleton is known (3rd row in Table 1). This result, to a considerable
extent, validates our hypothesis that the skeleton can be used to reduce the
search space of the optimization procedure and therefore boost the performance
of differentiable MAG learning.

On the contrary, we also find that when the ground-truth skeleton is unknown
and the used skeleton is learned by FCI (4th row in Table 1), the improvement
becomes modest because many erroneous/missing edges induced by FCI propagate
to the subsequent optimization procedure. This indicates the practical hurdle
of using the skeleton to facilitate differentiable causal discovery.

Indeed, due to finite samples, learning a high-quality skeleton is
challenging. Simply putting a “point estimation” of the skeleton as a hard
constraint on the optimization procedure can result in considerable missing
edges (i.e., false negatives) and finally impair performance.

Table 1. Boosting differentiable causal discovery with skeleton information.
(avg. on ten datasets with 50-100 variables.)

Method | Skeleton F1 | Arrowhead F1 | Tail F1
---|---|---|---
ABIC | 0.84 | 0.76 | 0.67
w/ True Skeleton | 1.00 +19% | 0.96 +26% | 0.96 +33%
w/ FCI | 0.87 +4% | 0.84 +11% | 0.71 +6%
SPOT | 0.91 +8% | 0.86 +13% | 0.78 +16%

The above preliminary results shed light on a possible solution that
synergistically combines differentiable causal discovery with skeleton
information, while alleviating error propagations. Recently, much research has
promoted the value of posterior distribution of DAGs (Cundy et al., 2021;
Lorch et al., 2021, 2022). Enlightened by these works, we advocate estimating
a posterior distribution (i.e., $p(S\mid\mathcal{D})$ where $S$ is the
skeleton and $\mathcal{D}$ is the dataset) over all skeletons to replace the
conventional “point estimation.” (i.e., $\arg\max_{S}p(S\mid\mathcal{D})$ )
The posterior effectively quantifies epistemic uncertainty and the degree of
confidence to any skeletons w.r.t. the given dataset. Nonetheless, unlike
DAGs, skeletons over causally insufficient variables do not have an explicit
form of likelihood function (i.e., $p(\mathcal{D},S)$).

To address these challenges, we propose SPOT (Skeleton Posterior-guided
OpTimization) as a two-phase framework for facilitating differentiable causal
discovery in the presence of latent confounders. SPOT first performs amortized
variational inference to estimate skeleton posterior, SPOT then employs a
novel optimization procedure for boosting subsequent differentiable causal
discovery. Specifically, SPOT leverages a recent advancement to amortize the
variational inference into the joint distribution of data and corresponding
skeleton $p(\mathcal{D},S)$ and alleviates the need of an explicit likelihood
function. Then, it concretizes the amortized inference with a supervised model
to estimate the skeleton posterior. To effectively facilitate the optimization
procedure, SPOT employs the skeleton posterior to stochastically update
variables in each optimization step, instead of deterministically updating
optimization variables with gradients.

As shown in the last row of Table 1, SPOT improves the performance of
differentiable causal discovery methods by a notable margin. In Sec. 4, we
conduct extensive experiments on various large-scale datasets and show that
SPOT delivers 8% improvement on skeleton F1 score and 13% and 16% improvement
on arrowhead and tail F1 scores, respectively, which are representative
metrics for evaluating the accuracy of the learned ADMGs. We also demonstrate
that the skeleton posterior estimated by SPOT is highly accurate compared to
other variational inference-based and non-parametric bootstrap-based
solutions. Finally, we explore the versatile applications of SPOT in MAG
learning for non-linear causally insufficient datasets and also in DAG
learning methods. Our empirical results indicate the strong potential of SPOT
in such scenarios.

In summary, we make the following contributions:

  1. (1)

Conceptually, we advocate a novel focus of using the skeleton posterior to
facilitate differentiable causal discovery in the presence of latent
confounders.

  2. (2)

Technically, we formulate the problem of skeleton posterior estimation under
amortized inference framework and propose a supervised learning-based solution
to estimate the skeleton posterior from observational data.

  3. (3)

On the basis of the skeleton posterior, we propose SPOT, a novel stochastic
optimization procedure, to facilitate differentiable causal discovery in the
presence of latent confounders, which incorporates the skeleton posterior in a
stochastic manner.

  4. (4)

Empirically, SPOT demonstrates superior performance on nearly all evaluation
metrics and various datasets, substantially improving upon its counterpart. We
also explore the extension of SPOT to other causal discovery settings,
including non-linear ADMG and DAG. Our results show that SPOT consistently
improves state-of-the-art methods in these settings. Our code will be made
publicly available at (Microsoft, 2024).

##  2\. Preliminary

As in many previous works (Tsirlis et al., 2018; Chen et al., 2021; Claassen
and Bucur, 2022; Bhattacharya et al., 2021), we focus on discovering MAG from
a linear Gaussian structural causal model (SCM) in the presence of latent
confounders and assume the absence of selection bias (no undirected edges in
ADMGs). In this section, we introduce preliminary knowledge of linear Gaussian
SCM with latent confounders, differentiable causal discovery and amortized
inference. In Appendix, we provide definitions to terminologies that are not
explicitly defined in the main text.

###  2.1. Linear Gaussian SCM with Latent Confounders

We first start with the definition of a linear Gaussian SCM without latent
confounders. Consider a linear SCM with $d$ observable variables parameterized
by a coefficient matrix $\delta\in\mathbb{R}^{d\times d}$. The SCM can be
written as

(1) $\small
V_{i}\leftarrow\sum_{V_{j}\in\textbf{PA}_{i}}\delta_{ji}V_{j}+\epsilon_{i}$

where $\textbf{PA}_{i}$ are the parents of $V_{i}$ and $\epsilon_{i}$ is a
noise term that is mutually independent of all other noise terms.

The noise term $\epsilon_{i}$ is mutually independent of others if and only if
$V_{i}$ has no latent confounder. Otherwise, $\epsilon_{i}$ is correlated with
$\epsilon_{j}$ if $V_{j}$ shares a latent confounder with $V_{i}$. When the
noise terms are Gaussian, the correlation can be expressed by a covariance
matrix $\beta=\mathbb{E}[\epsilon\epsilon^{T}]$ and the joint distribution
marginalized over observable variables $\bm{V}_{\text{O}}$ forms a zero-mean
multivariate Gaussian distribution with covariance matrix as
$\Sigma=(I-\delta)^{-T}\beta(I-\delta)^{-1}$. The induced graph ${G}$ is an
ADMG and contains two types of edges, including directed edges ($\rightarrow$)
and bidirected edges ($\leftrightarrow$), which implies two adjacency
matrices, namely $D$ and $B$. In the special case where there is no latent
confounder, the ADMG is equivalent to a DAG and the adjacency matrix of
bidirected edges $B$ is all zeros. $V_{i}\to V_{j}$ exists in ${G}$ and
$D_{ij}=1$ if and only if $\delta_{ij}\neq 0$. $V_{i}\leftrightarrow V_{j}$
exists in ${G}$ and $B_{ij}=1$ if and only if $\beta_{ij}\neq 0$. It is
commonly assumed that $\delta_{ij}=\delta_{ji}=\beta_{ij}=\beta_{ji}=0$ if and
only if $D_{ij}=D_{ji}=B_{ij}=B_{ji}=0$.

###### Remark 0.

In the context of ADMGs, the noise term $\epsilon$ does not exclusively
represent the exogenous variables. In particular, $\epsilon_{i}$ in Eqn. 1
includes both the exogenous variable and the confounding effect from $V_{u}$
(the latent confounder). Similarly, the noise term on another variable
$V_{j}$, $\epsilon_{j}$, also includes its own exogenous variable and the
confounding effect from $V_{u}$. Therefore, $\epsilon_{i}$ and $\epsilon_{j}$
are jointly influenced by the confounding effect from $V_{u}$, which makes
them dependent. We provide a more detailed explanation below.

Consider the causal graph $V_{u_{i}}\rightarrow V_{i}\leftarrow
V_{u}\rightarrow V_{j}\leftarrow V_{u_{j}}$. Without loss of generality on the
linear setting, we have $V_{i}=V_{u_{i}}+aV_{u}$ and $V_{j}=V_{u_{j}}+bV_{u}$,
where $a,b$ are non-zero coefficients. Given that only $V_{i}$ and $V_{j}$ are
observed, the resulting ADMG should be $V_{i}\leftrightarrow V_{j}$, based on
Sec. 2.2 of (Zhang, 2008). In Eqn. 1, $V_{i}=\sum_{V_{k}\in
Pa(V_{i})}\delta_{ki}V_{k}+\epsilon_{i}$ where $Pa(V_{i})$ are the parents of
$V_{i}$ in the ADMG. Since there is a bidirectional edge $V_{i}\leftrightarrow
V_{j}$, neither $V_{i}$ is the parent of $V_{j}$ nor $V_{j}$ the parent of
$V_{i}$. Hence, $Pa(V_{i})=Pa(V_{j})=\emptyset$. Thus, we have
$V_{i}=\epsilon_{i}$ and $V_{j}=\epsilon_{j}$. However, since $V_{i}$ and
$V_{j}$ are dependent, it follows that $\epsilon_{i}$ and $\epsilon_{j}$ must
also be dependent.

According to the causal Markov and faithfulness assumption, ancestral ADMG,
MAG and skeleton can be defined as follows.

###### Definition 0 (Ancestral Acyclic Directed Mixed Graph).

${G}$ is an ancestral ADMG if it is a mixed graph with directed edges
($\rightarrow$) and bidirected edges ($\leftrightarrow$) and contains no
directed cycles or almost directed cycles (Zhang, 2008).

###### Definition 0 (Maximal Ancestral Graph).

An ancestral ADMG $G$ is a Maximal Ancestral Graph if for each pair of non-
adjacent nodes, there exists a set of nodes that make them m-separated (Zhang,
2008).

###### Definition 0 (Skeleton).

An undirected graph $S$ is a skeleton of an MAG $G$ if $S$ is obtained from
$G$ by replacing all directed edges with undirected edges.

Therefore, two nodes $V_{i},V_{j}$ are adjacent in the skeleton if and only if
$\forall\bm{Z}\subseteq\bm{V}_{\text{O}}\setminus\\{V_{i},V_{j}\\},V_{i}\not\perp
V_{j}\mid\bm{Z}$, where $\perp$ and $\not\perp$ denote conditional
independence and dependence, respectively. Hence, it is clear that given a
dataset whether two nodes are adjacent is not influenced by the adjacencies of
other nodes.

Parameter Estimation of Linear Gaussian SCM. For DAGs, one easily can estimate
the parameters of the SCM by solving a simple least squares regression
problem. However, for ADMGs, the estimation of the parameters is more
challenging due to the presence of latent confounders. Drton et al. (2009)
proposed an iterative algorithm to estimate the parameters of the SCM called
Residual Iterative Conditional Fitting (RICF). RICF generally works for bow-
free ADMGs and ancestral ADMGs, which are our focus in this paper, are special
cases of bow-free ADMGs (Bhattacharya et al., 2021). This algorithm
iteratively fits the SCM to the data and updates the covariance matrix. The
algorithm is guaranteed to converge to a local minimum when the corresponding
ADMGs are ancestral.

###  2.2. Differentiable Causal Discovery

Score-based methods aim to maximize the score (e.g., log-likelihood) of the
graph on the given data, which can be written in the following form:

(2) $\small\arg\max_{{G}}f(G)\text{ s.t. }{G}\text{\leavevmode\nobreak\ is
acyclic}$

where $f(\cdot)$ is the score function. Given the acyclicity constraint, the
optimization process is combinatorial. Recently, differentiable methods
reformulate this combinatorial constraint into a constraint
$h_{\text{DAG}}(W)=\text{tr}(e^{W\circ W})-d$ such that

(3) $\small h_{\text{DAG}}(W)=0\iff{G}\text{\leavevmode\nobreak\ is acyclic}$

where $d=|\bm{V}_{\text{O}}|$ and $W$ is a weighted adjacency matrix of ${G}$.
For linear SCM, an element in $W$ represents the linear coefficient. As a
continuous optimization with equality constraints, the augmented Lagrangian
method (ALM) is commonly used to convert the constrained optimization problem
into several unconstrained subproblems and use standard optimizers to solve
them separately (Zheng et al., 2018).

Despite the success of differentiable methods for DAG learning, the algebraic
characterization in Eqn. 3 cannot be directly applied to ADMGs (and MAGs). As
aforementioned, ADMGs requires two adjacency matrices, $D$ and $B$, to
represent directed edges and bidirected edges, respectively. To extend the
algebraic characterization to ADMGs, ABIC (Bhattacharya et al., 2021) modified
it as

(4) $h_{\text{ADMG}}(D,B)=\text{tr}(e^{D})-d+\text{sum}(e^{D}\circ B)$

where $\circ$ denotes the Hadamard product, i.e., element-wise multiplication,
and $\text{sum}(\cdot)$ is the sum of all elements in a matrix. It has been
proved that $h_{\text{ADMG}}(D,B)=0\iff{G}\text{\leavevmode\nobreak\ is
ancestral ADMG}$. In a nutshell, $\text{tr}(e^{D})-d$ implies the standard
directed acyclicity constraint and the last term $\text{sum}(e^{D}\circ B)$ is
a term to ensure that the bidirected edges does not introduce “almost directed
cycles”. In the following, we use $h$ to denote $h_{\text{ADMG}}$ for
simplicity.

Input: Dataset $\mathcal{D}$

Output: Adjacency Matrix $D,B$

1 Initialize $\delta^{(1,1)},\beta^{(1,1)}$;

2 Define $h(\delta,\beta)$ according to Eqn. 4;

3 foreach  _$t_{1}=1,\cdots,T_{\text{ALM}}$_ do

4 foreach  _$t_{2}=1,\cdots,T_{\text{RICF}}$_ do

5 Update pseudovariables by $\delta^{(t_{1},t_{2})},\beta^{(t_{1},t_{2})}$;

6 Constitute $f(\delta,\beta)$;

7 $\delta^{(t_{1},t_{2}+1)},\beta^{(t_{1},t_{2}+1)}\leftarrow\arg\min
f(\delta,\beta)$;

8

9 end foreach

10 Update $\alpha^{(t_{1})},\rho^{(t_{1})},\lambda^{(t_{1})}$;

11

12 end foreach

13$D\leftarrow\bm{1}(|\delta^{(T_{\text{ALM}},T_{\text{RICF}})}>\omega|)$;

14 $B\leftarrow\bm{1}(|\beta^{(T_{\text{ALM}},T_{\text{RICF}})}>\omega|)$;

15 return $D,B$;

Algorithm 1 ABIC (Bhattacharya et al., 2021)

Now, we introduce ABIC (Bhattacharya et al., 2021) in Alg. 1. To ease
presentation, we omit technical details (e.g., calculations of
$\alpha,\rho,\lambda$, pseudovariables, and stopping criteria) and refer
readers to the original paper for the full version. The ABIC algorithm aims to
maximize the likelihood of ADMG over the observational data $\mathcal{D}$,
ensuring $h_{\text{ADMG}}(\delta,\beta)=0$. It employs two nested loops: the
outer loop applies the augmented Lagrangian method (ALM) to convert the
equality-constrained problem into $T_{\text{ALM}}$ steps of unconstrained ones
(lines 3–10 in Alg. 1). After ALM, $\delta$ and $\beta$ represent the ADMG’s
linear coefficients and covariance, respectively. The adjacency matrix is
derived by applying a threshold to remove negligible coefficients (lines
22–23), yielding the ancestral ADMG. If a MAG is required, the maximal
ancestral projection, a standard post-processing procedure, is performed.

Inside the inner loop, Alg. 1 iteratively adjust the objective function (lines
5–6) and solve the unconstrained optimization problem accordingly (line 7). In
particular, the objective function $f(\delta,\beta)$ (line 6) consists of the
least square loss to fit the ADMG to the data, the constraint
$h(\delta,\beta)$, and a regularization term to enforce sparsity. The form of
$f(\delta,\beta)$ is updated according to the current pseudovariables (line
5). Then, in the optimization phase (line 7), optimization algorithms are used
in an “out-of-the-box” manner to find a minimum of $f(\delta,\beta)$. However,
due to the complexity and non-convex nature of $f(\delta,\beta)$, it is
challenging for standard optimization techniques to find a plausible solution,
especially when the ADMG is large, as we will show in Sec. 4.

Figure 1. SPOT workflow.

###  2.3. Amortized Inference for Causal Discovery

In the conventional variational inference framework, the posterior
distribution is approximated by a parametric distribution
$q_{\phi}(G|\mathcal{D})$ with parameters $\phi$. Here, $G$ is the causal
graph and $\mathcal{D}$ is the observed data in the context of causal
discovery. The parameters are learned by minimizing the KL divergence between
the true posterior $p(G|\mathcal{D})$ and the parametric distribution
$q_{\phi}(G|\mathcal{D})$. Typically, variational inference requires the
likelihood function $p(\mathcal{D},G)$ to compute the KL divergence. Recently,
amortized inference has been proposed to alleviate the need of an explicit
likelihood function by introducing a simulator that can yield samples from the
joint distribution $p(G,\mathcal{D})$ (Ambrogioni et al., 2019). Lorch et al.
(2022) follows this line of work and introduces AVICI, which amortizes the
inference process by introducing a sampler that can generate samples from the
simulator and capture the domain-specific inductive biases that would have
been difficult to characterize (e.g., gene regulatory networks or non-linear
functions of random Fourier features).

Recently, SCL (supervised causal learning) has emerged as a promising paradigm
for causal discovery (Ke et al., 2022; Ma et al., 2022a; Dai et al., 2023).
Compared to standard machine learning models, SCL models features many
advantages tailored for causal discovery (e.g., invariant to the permutation
of variables). It is worth noting that the concept of “permutation invariance”
is different from the one in machine learning setting. One dataset (a $M\times
N$ matrix) is a sample to the SCL model. It implies that the model is
permutation in- and equivariant with respect to the observation and variable
dimensions of the provided dataset, respectively (Ke et al., 2022). SCL models
can be trained on samples from a known simulator. In essence, SCL can be
viewed as an instantiation of amortized inference when the simulator is in
line with the underlying causal mechanism of the observational data.

##  3\. SPOT

As aforementioned, skeleton information can be leveraged to enhance the
optimization procedure of differentiable causal discovery in the presence of
latent confounders. However, a “point estimation” of the skeleton is prone-to-
error. In this section, we propose SPOT (Skeleton Posterior-guided
OpTimization) to leverage a probabilistic skeleton posterior to guide the
optimization of ADMGs.

As shown in Fig. 1, SPOT consists of four steps: ➀ a static simulator first
generates an initial set of data/skeleton pairs to obtain an initial amortized
inference model (blue one in Fig. 1). ➁ while the initial amortized inference
model already provides a good estimation of the skeleton posterior, SPOT can
optionally use a dynamic simulator to generate more in-domain training
instances with respect to the input dataset. ➂ with samples from the dynamic
simulator, SPOT adapts the initial amortized inference model to obtain an
updated model (yellow one in Fig. 1) for skeleton posterior inference. ➃ SPOT
uses the inferred skeleton posterior to enhance the optimization procedure of
differentiable causal discovery.

Conceptual Complexity. While our solution provides an additional layer of
complexity compared to standard differentiable causal discovery, we argue that
the level of complexity is generally manageable and comparable to other causal
discovery algorithms. First, the two stages in our pipeline are independent
and do not involve joint training. Second, constraint-based causal discovery
algorithms, such as FCI, also involve two similar steps: learning the skeleton
in the first step and orienting the edges in the second step.

###  3.1. Skeleton Variational Objective

Let $\mathcal{D}=\\{\bm{x}^{1},\cdots,\bm{x}^{n}\\}\sim p(V)$ be the
observational dataset, where $\bm{x}^{i}$ is sampled from the joint
distribution $p(V)$ and $n$ is sample size. We aim to approximate the
posterior over skeletons $p(S\mid\mathcal{D})$ with a variational distribution
$q(S;a)$ where $S$ is the (symmetric) adjacency matrix of the skeleton and $a$
is the variational parameters. Thus, the skeleton estimation can be decomposed
into a set of edge estimations (i.e., the probability of adjacency)
independently. Then, the variational family of $q(S;a)$ is defined as

(5) $\small q(S;a)=\prod_{i<j}q(S_{ij};a_{ij})\text{ with
}S_{ij}\sim\text{Bernoulli}(a_{ij}).$

In this regard, we aim to find an inference model $f_{\phi}(\mathcal{D})$ that
predicts $a$. This procedure can be attained by minimizing the expected
forward KL divergence from $p(S\mid\mathcal{D})$:

(6)
$\small\min_{\phi}\mathbb{E}_{p(\mathcal{D})}\left[\text{KL}\left(p(S\mid\mathcal{D})\mid\mid
q(S;f_{\phi}(\mathcal{D}))\right)\right]$

Following the principle of amortized inference (Ambrogioni et al., 2019), we
amortize the inference and rewrite the objective as

(7)
$\displaystyle\mathbb{E}_{p(\mathcal{D})}\text{KL}\left(p(S\mid\mathcal{D})\mid\mid
q(S;f_{\phi}(\mathcal{D}))\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{p(\mathcal{D})}\mathbb{E}_{p(S\mid\mathcal{D})}\left[\log
p(S\mid\mathcal{D})-\log q(S;f_{\phi}(\mathcal{D}))\right]$ $\displaystyle=$
$\displaystyle-\mathbb{E}_{p(S)}\mathbb{E}_{p(\mathcal{D}\mid S)}\left[\log
q(S;f_{\phi}(\mathcal{D}))\right]+\text{ const.}$ $\displaystyle=$
$\displaystyle-\mathbb{E}_{p(\mathcal{D},S)}\left[\log
q(S;f_{\phi}(\mathcal{D}))\right]+\text{ const.}$

Since the constant does not depend on $\phi$, we can merely minimize
$\mathcal{L}(\phi)\coloneqq-\mathbb{E}_{p(\mathcal{D},S)}\left[\log
q(S;f_{\phi}(\mathcal{D}))\right]$ to obtain $\phi$. In other words, the
problem is recast to train a predictive model
$f_{\phi}(\cdot):\mathcal{D}\mapsto S$ over the distribution
$p(\mathcal{D},S)$.

Estimate $p(\mathcal{D},S)$ via Static Simulator. Naturally, we can estimate
$p(\mathcal{D},S)$ by using a simulator (e.g., Erdős-Rényi random graph model)
then generate the corresponding dataset with pre-defined functional forms. The
simulator serves a direct sampler of $p(\mathcal{D},S)$ if the test data
$\mathcal{D}$ (i.e., the input dataset on which we need to conduct causal
discovery) is in the same distribution of the static simulator. In other
words, the pair of $\mathcal{D}$ and the skeleton $S$ is known to be drawn
from $p$ a priori.

Estimate $p(\mathcal{D},S)$ via Dynamic Simulator. To further feed the model
with more in-domain training instances, we propose to use nonparametric
bootstrap method (Friedman et al., 1999) to estimate a MAG $\hat{G}$ with a
random subset of observational data. Then, we fit the parameters of the
underlying SCM (i.e., $\delta,\beta$; see Sec. 2.1) on the given test
observational data using the RICF algorithm (Drton et al., 2009). Finally, we
regenerate the new observational data $\hat{\mathcal{D}}$ from the fitted SCM.
By repeating the above procedure, we obtain a set of data/MAG pairs
$\\{(\hat{\mathcal{D}}_{1},\hat{G}_{1}),\cdots\\}$ from which we can derive
the samples of $p(\mathcal{D},S)$ as
$\\{(\hat{\mathcal{D}}_{1},\hat{G}_{1}+\hat{G}_{1}^{T}),\cdots\\}$.111$G$ and
$G^{T}$ are the adjacency matrix and its transpose,
$\hat{G}_{1}+\hat{G}_{1}^{T}$ is the skeleton adjacency matrix. We note that,
while the way latent variables affect the observed variables is implicit, in
the linear Gaussian setting, we can characterize the marginalized distribution
over observed variables as a zero-mean multivariate Gaussian distribution with
a covariance matrix defined as $\Sigma=(I-\delta)^{-T}\beta(I-\delta)^{-1}$,
where $\delta$ and $\beta$ are parameters fitted by RICF. In this way, we
presume this dynamic (dataset-dependent) simulator would provide more relevant
training instances with respect to the input dataset $\mathcal{D}$.

###### Remark 0.

In summary, the static simulator features a robust model that generalizes well
to the test data, even with potential domain shift, as validated in (Lorch et
al., 2022; Ke et al., 2022). On the other hand, the dynamic simulator provides
more in-domain training instances, which may lead to a better performance at a
fairly lightweight cost on the runtime, as will be shown shortly in Sec. 3.2.

###  3.2. Skeleton Posterior Inference

Figure 2. Skeleton Inference Procedure.

In light of Sec. 3.1, the predictive model $f_{\phi}(\cdot):\mathcal{D}\mapsto
S$ can be instantiated by any supervised learning model. Enlightened by the
prosperous progress in SCL, we propose to use SCL model to enable the
amortized inference and estimate the skeleton posterior $p(S\mid\mathcal{D})$.
Below, we elaborate the design considerations of instantiating the SCL model
in SPOT.

➀ Model Architecture. Existing SCL methods either use an end-to-end model
(e.g., Transformer) to directly predict the adjacency matrix (Ke et al., 2022)
or use a simple classifier (e.g., xgboost) to predict the local structure
(e.g., the adjacency of two variables or the v-structure) (Ma et al., 2022a;
Dai et al., 2023). While the design of SPOT is agnostic to the model
architecture, we anticipate to use a simple classifier, as it is more
effective and efficient for estimating $p_{ij}$ in practice. Hence, in our
implementation, we adopt a ML4S-like cascade model (Ma et al., 2022a) to
predict the adjacency and constitute the skeleton posterior. We formulate the
problem of efficient skeleton posterior inference as predicting a series of
adjacencies and construct the skeleton posterior. Specifically, we aim to
estimate

(8) $\small
f_{\phi}(\mathcal{D})_{ij}=\mathbb{E}_{p(\mathcal{D},S)}\left[p(V_{i}-V_{j}\in
S\mid\mathcal{D})\right],$

We focus on the probability $p(V_{i}-V_{j}\mid\mathcal{D})$ for the adjacency
of $V_{i},V_{j}$ given $\mathcal{D}$, which helps constitute the skeleton
posterior (we use $p_{ij}$ as an abbreviation of
$p(V_{i}-V_{j}\mid\mathcal{D})$). Besides, we also reproduce Ke et al.
(2022)’s Transformer-based model, and the results show that the cascade model
is more effective and efficient.

➁ Objective Refinement. Instead of standard causal discovery objectives (e.g.,
the higher F1 score), we prefer higher recall over precision to prevent missed
edges. This is because the spurious edges will be offloaded to and eventually
removed by differentiable causal discovery algorithms. Specifically, even
though the SCL model may predict a spurious edge and corresponding parameters
$\delta_{ij},\beta_{ij}$ tend to be non-zero, the differentiable causal
discovery algorithm will gradually move the coefficients to zero to minimize
the objective function. In this regard, the spurious edges are eventually
removed. However, the missed edges are hard to recover due to the sparse
nature of the ADMG. To this end, we can either set a conservative threshold
(subject to the specific design of the SCL model) or apply label smoothing
trick to the training instances to prevent the model from predicting zero
probability for the true edges.

➂ Efficient Training. When using the dynamic simulator, the model has to be
trained from scratch with the training instances sampled from non-parametric
bootstrap. However, this would impose a heavy sampling and training overheads
in the runtime. To alleviate this issue, we adopt a widely-used domain
adaptation paradigm to reduce the training overheads by leveraging the static
model. Specifically, given a dataset, we only sample a few training instances
from the dynamic simulator. Then, we use the output and intermediate values of
the static model to augment the original features of the samples from the
dynamic simulator and train an adapted dynamic model efficiently. In this way,
the dynamic model can rapidly adapt to the new dataset with few-shot training
instances and yields a better performance with mild cost.

###  3.3. Posterior-guided Optimization

Non-adjacency between $V_{i},V_{j}$ implies
$\delta_{ij}=\delta_{ji}=\beta_{ij}=\beta_{ji}=0$. When ground-truth skeleton
is available, we enforce extra equality constraints to optimize. However, the
true skeleton is unattainable in practice, necessitating estimating a skeleton
from observational data. We advocate using the skeleton posterior for
optimization as it encodes epistemic uncertainty.

It is unclear how to incorporate a skeleton posterior in the optimization
procedure. The estimated skeleton posterior from Sec. 3.2 is a continuous
adjacency matrix $p_{ij}\in[0,1]$. Thus, we cannot derive meaningful
constraints for coefficients/covariances. Formally,
$p(V_{i}-V_{j}\mid\mathcal{D})=p_{ij}$ implies the probability of the union of
(disjoint) events $P(\delta_{ij}\neq 0\cup\delta_{ji}\neq 0\cup\beta_{ij}\neq
0\mid\mathcal{D})$ equals $p_{ij}$. Higher posterior probability indicates
higher likelihood of non-zero coefficients. In this regard, $p_{ij}$ cannot
provide any additional information regarding the value of non-zero
coefficients. For example, $\delta_{ij}=0.1$ and $\delta_{ij}=0.9$ both comply
with arbitrary positive posterior $p_{ij}>0$ equally well. To account for the
probabilistic nature of the skeleton posterior, we propose the posterior-
guided optimizer for structure learning as a replacement of the standard
optimizer in differentiable causal discovery (e.g., the gradient descent or
the L-BFGS).

Input: Skeleton Posterior $p$, Optimization Variables $\delta,\beta$,
Objective Function $f$, ALM Step $t$, Temperature Constant $c$

Output: Optimized Variables $\delta,\beta$

1 $\delta^{*},\beta^{*}\leftarrow\arg\min f(\delta,\beta)$;

2 forall  _$\delta_{ij}$ and $\beta_{ij}$_ do

3 if  _$\delta_{ij}\times\frac{\partial f(\delta,\beta)}{\partial\delta_{ij}}
>0$_ then

4 $\delta_{ij}\leftarrow\delta_{ij}^{*}$

5 else

6 Update $\delta_{ij}\leftarrow\delta_{ij}^{*}$ with probability
$(p_{ij}+c)^{\frac{1}{t+1}}$;

7

8 end if

9

10 end forall

11return $\delta,\beta$;

Algorithm 2 Skeleton Posterior-guided Optimizer

Alg. 2 presents the skeleton posterior-guided optimizer. In accordance with
many standard differentiable causal discovery methods, we first optimize the
objective function $f(\delta,\beta)$ with the standard optimizer (e.g.,
gradient descent) to obtain the optimal coefficients and covariances
$\delta^{*},\beta^{*}$ (line 1). Then, instead of directly taking
$\delta^{*},\beta^{*}$ as the final result, we introduce a stochastic update
scheme to update $\delta,\beta$ based on the skeleton posterior $p$ (lines
2–8) as follows:

(9) $\small
P(\delta_{ij}\leftarrow\delta_{ij}^{*})=\begin{cases}1&\delta_{ij}\times\frac{\partial
f(\delta,\beta)}{\partial\delta_{ij}}>0\\\
(p_{ij}+c)^{\frac{1}{t+1}}&\text{otherwise}\end{cases}$

If $\delta_{ij}\times\frac{\partial f(\delta,\beta)}{\partial\delta_{ij}}>0$,
the update is accepted unconditionally (line 3–4); otherwise, it is accepted
with a probability of $(p_{ij}+c)^{\frac{1}{t+1}}$ (line 12–15; Fig. 3). If
rejected, we keep $\delta_{ij}$ unchanged (line 13). $\beta$ is also optimized
in a similar way.

Conceptually, the overall procedure can be viewed as a modified simulated
annealing process (Kirkpatrick et al., 1983). $\delta_{ij}\times\frac{\partial
f(\delta,\beta)}{\partial\delta_{ij}}>0$ implies that the coefficient is
“moving to zero” (i.e., a weakened edge) in this update. Given that real-world
ADMGs are often sparse, “moving-to-zero” update is encouraged and therefore
accepted unconditionally. This design consideration is thus an analogy of the
unconditional update in simulated annealing.

Figure 3. The probability of accepting an update for different $p_{ij}$ and
$t$ with $c=0.1$.

In contrast, if the coefficient is “moving to non-zero” (i.e., a strengthened
edge), the update is accepted with a probability computed by $p_{ij}$ and $t$
jointly. Recall that $p_{ij}$ is the probability of the edge existence. The
update is thus more encouraged with higher $p_{ij}$. $\frac{1}{t}$ is a
analogy of the temperature parameter in simulated annealing. As shown in Fig.
3, an update is more likely to be rejected with lower $t$ and higher $t$
yields $(p_{ij}+c)^{\frac{1}{t}}\approx 1$. With the increase of $t$, the
probability of accepting an update is gradually similar to what would be in
the vanilla optimization procedure even with a low $p_{ij}$. In this way, the
optimization procedure is gradually stabilized. Furthermore, it is worth
noting that the scheme is robust to the accuracy of the skeleton posterior.
For spurious edges, even if $p_{ij}$ is high, the optimization procedure of
differentiable causal discovery is still able to move the coefficients to zero
given the unconditional update for “edge-weakening” updates.

Table 2. End-to-end comparison (avg. on ten datasets with nodes
$\in[50,100]$).

Method | Skeleton | Arrowhead | Tail | #Failed
---|---|---|---|---
F1 | TPR | FDR | F1 | TPR | FDR | F1 | TPR | FDR | Datasets
FCI | 0.84 $\pm$ 0.06 | 0.79 $\pm$ 0.10 | 0.10 $\pm$ 0.05 | 0.57 $\pm$ 0.05 | 0.78 $\pm$ 0.15 | 0.54 $\pm$ 0.04 | 0.53 $\pm$ 0.18 | 0.46 $\pm$ 0.20 | 0.31 $\pm$ 0.15 | 0/10
RFCI | 0.85 $\pm$ 0.07 | 0.77 $\pm$ 0.09 | 0.05 $\pm$ 0.03 | 0.62 $\pm$ 0.09 | 0.76 $\pm$ 0.12 | 0.45 $\pm$ 0.12 | 0.51 $\pm$ 0.15 | 0.42 $\pm$ 0.17 | 0.27 $\pm$ 0.14 | 0/10
ICD | 0.82 $\pm$ 0.09 | 0.81 $\pm$ 0.07 | 0.17 $\pm$ 0.11 | 0.56 $\pm$ 0.10 | 0.82 $\pm$ 0.06 | 0.57 $\pm$ 0.12 | 0.48 $\pm$ 0.15 | 0.41 $\pm$ 0.15 | 0.39 $\pm$ 0.18 | 1/10
M3HC | 0.73 $\pm$ 0.09 | 0.59 $\pm$ 0.12 |  0.03$\pm$ 0.03 | 0.56 $\pm$ 0.15 | 0.46 $\pm$ 0.16 | 0.27 $\pm$ 0.11 | 0.32 $\pm$ 0.11 | 0.24 $\pm$ 0.10 | 0.48 $\pm$ 0.13 | 0/10
GPS | 0.73 $\pm$ 0.10 | 0.81 $\pm$ 0.11 | 0.33 $\pm$ 0.11 | 0.62 $\pm$ 0.09 | 0.85 $\pm$ 0.11 | 0.50 $\pm$ 0.10 | 0.43 $\pm$ 0.12 | 0.62 $\pm$ 0.19 | 0.66 $\pm$ 0.10 | 0/10
ABIC | 0.84 $\pm$ 0.04 | 0.89 $\pm$ 0.07 | 0.19 $\pm$ 0.07 | 0.76 $\pm$ 0.07 | 0.85 $\pm$ 0.15 | 0.29 $\pm$ 0.11 | 0.67 $\pm$ 0.07 | 0.80 $\pm$ 0.15 | 0.40 $\pm$ 0.07 | 0/10
SPOT |  0.91 $\pm$ 0.03 |  0.94 $\pm$ 0.04 | 0.11 $\pm$ 0.04 |  0.86 $\pm$ 0.05 |  0.89 $\pm$ 0.10 |  0.16 $\pm$ 0.06 |  0.78 $\pm$ 0.09 |  0.86 $\pm$ 0.13 |  0.27 $\pm$ 0.13 | 0/10

Convergence. In the context of differentiable causal discovery, the concept of
convergence is three-fold. First, in accordance with the ABIC (Bhattacharya et
al., 2021), asymptotically, convergence to the global optimum of its objective
function implies that the corresponding ADMG is within the Markov equivalence
class of the true ADMG. Second, given the non-convexity of the objective
function, it is not guaranteed that the optimization procedure converges to
the global optimum; instead, it may result in a local minimum or a saddle
point. Finally, we discuss the impact of the stochastic update scheme on
convergence to the local minima in the following sense: Since the probability
of an update reaches 1 (see Fig. 3), the optimization procedure becomes
equivalent to the vanilla optimizer when $t$ is sufficiently large. In this
regard, the stochastic update scheme does not negatively affect the
convergence to the local minima. In fact, our empirical results show that our
stochastic update scheme yields much better results than ABIC and also
surpasses other SOTA methods.

Probabilistic Update vs. Regularization. A natural question arises: why not
use regularization (e.g., a penalty term) to incorporate the skeleton
posterior? We argue that the penalty term may face several hurdles in
practice. First, the objective function is often unstable due to the acyclic
term (Eqn. 4). Its value can range from $1\times 10^{-5}$ to $1\times 10^{5}$
over different stages of the optimization. Second, as aforementioned,
$\delta,\beta$ are not directly comparable with posterior probability. This
necessitates considerable efforts for implementation, thus we did not adopt
it.

##  4\. Evaluation

In evaluation, we aim to evaluate the performance of SPOT from three
perspectives: ➊ end-to-end causal discovery performance; ➋ standalone skeleton
posterior inference performance; ➌ extension to other differentiable causal
structure learning algorithms. We also compare SPOT with a wide range of
baselines.

Data Generation. We follow the similar setting from previous works (Bernstein
et al., 2020; Bhattacharya et al., 2021) to generate synthetic datasets. We
generates ADMGs from Erdős-Rényi (ER) model with various nodes whose average
indegree is $\in[1,1.5]$. Each ADMG contains $5-15\%$ bidirected edges. Then,
it is parameterized by the same rules of (Bhattacharya et al., 2021): if
$V_{i}\to V_{j}$, $\delta_{ij}$ is uniformly sampled from $\pm[0.5,2.0]$; if
$V_{i}\leftrightarrow V_{j}$, $\beta_{ij}=\beta_{ji}$ is uniformly sampled
from $\pm[0.4,0.7]$; and $\beta_{ii}$ is uniformly sampled from $\pm[0.7,1.2]$
and add $\text{sum}(|\beta_{i,-i}|)$ to ensure positive definiteness. For each
ADMG, we generate 1000 samples.

Baselines. We compare our method with a wide range of baselines, including
constraint-based methods (FCI (Spirtes et al., 2000), RFCI (Colombo et al.,
2012) and ICD (Rohekar et al., 2021)), score-based methods (M3HC (Tsirlis et
al., 2018) and GPS (Claassen and Bucur, 2022)) and also the differentiable
causal discovery methods (ABIC (Bhattacharya et al., 2021)). The comparison is
made on synthetic datasets with 50–100 nodes. We observe that ABIC
occasionally outputs cyclic graphs with the default threshold, which is
unwarranted for maximal ancestral projection. In such cases, we use the
minimal threshold that produces an acyclic graph. We also tentatively tried
AGIP (Chen et al., 2021) algorithm and excluded it from comparison, because
its preprocessing step is considerably slow (i.e., taking several days) for
large datasets. Additionally, we omit algorithms lacking open-source
implementation (e.g., GreedySPo (Bernstein et al., 2020)). Since we focus on
linear Gaussian data, we omit the algorithms that do not comply with this
assumption (e.g., (Salehkaleybar et al., 2020; Wang and Drton, 2023)) in the
main comparison. In Sec. 4.3, we explore the effectiveness of SPOT on other
settings.

Running Time. Many algorithms can take a long time to run on large datasets
and the convergence condition is not always satisfied. To make the comparison
fair, we set a timeout of 24 CPU hours for each algorithm. For algorithms that
take more than 24 CPU hours, we halt the algorithm and use the best-so-far
graph as their output.

Metrics. Under linear Gaussian SCM (Spirtes et al., 2000), we can only up to
identify a Markov equivalence class. To make all baselines comparable, for
those algorithms that output an MAG, we convert it to a PAG (Partial Ancestral
Graph) and compare it with the true PAG. Following (Bhattacharya et al.,
2021), we report the F1 score, TPR (True Positive Rate), and FDR (False
Discovery Rate) on the skeleton, arrowheads, and tails.

###  4.1. End-to-end Comparison

We report the results in Table 2 on ten datasets with nodes sampled between
$[50,100]$. We observe that the large datasets pose a considerable challenge
to all methods while SPOT accurately identifies the underlying causal
structure. ICD fails on one dataset; M3HC and GPS frequently terminate with
convergence warnings on these datasets. We also observe that constraint-based
methods are usually more powerful at identifying the correct skeleton while
all methods have notable difficulties to accurately identify the arrowheads
and tails. In contrast, SPOT consistently manifests its superiority across
nearly all criteria and yields more precise and stable estimation of
skeletons, arrowheads, and tails. Over the ten datasets, it significantly
improves the performance of its counterpart, ABIC (with p-value of 0.008). We
also investigate the criterion (i.e., Skeleton’s FDR) on which SPOT is sub-
optimal while M3HC is the best. Having said that, we observe that M3HC also
suffers from the lowest TPR and F1 scores on skeleton. Hence, we conclude that
M3HC may be too strict on confirming an edge while SPOT is more robust and
attains the best F1 scores.

It is worth noting that SPOT’s strong performance can be attributed to its
ability to effectively balance precision and recall in estimating the
underlying causal relationships. Additionally, the more consistent performance
of SPOT could potentially make it more suitable for various real-world
applications where accuracy and robustness are crucial factors.

Figure 4. Comparison of ABIC and SPOT on different node sizes.

Comparison by Node Size. In addition to the aggregated results, we also report
the edge-wise F1 scores under different node sizes in Fig. 4 where the data
point of each node size is the average of three datasets and the error bar is
the standard deviation. First, we observe that SPOT consistently outperforms
all baselines across all node sizes except for 20-node graphs. For 20-node
graphs, potentially due to randomness, though ABIC slightly outperforms SPOT,
the difference is not statistically significant ($p=0.51$). In other settings,
especially for large graphs, SPOT outperforms ABIC (the second-best method) by
a large margin (e.g., $p=0.001$ on 90-node graphs). Second, the performance of
ABIC significantly drops on 100-node graphs where it fails to learn a
meaningful causal structure. In contrast, SPOT, by incorporating the skeleton
posterior in the optimization, manifests stable and strong performance across
all node sizes with a small standard deviation.

###  4.2. Effectiveness of Skeleton Posterior Inference

Through this experiment, our goal is to determine if SPOT is a reliable
skeleton posterior estimator. We assess the posterior quality by calculating
the KL-divergence between the estimated posterior and ground-truth skeleton.
Following earlier research (Lorch et al., 2022), we report AUPRC and AUROC for
edge predictions to evaluate the quality of the estimated posterior.

Using the same datasets from Sec. 4.1, we compare SPOT to three baselines:
AVICI (Lorch et al., 2022), FCI* (Nonparametric Bootstrap FCI), and RFCI*
(Nonparametric Bootstrap RFCI). We tried to include more baselines (e.g., DiBS
(Lorch et al., 2021) and N-ADMG (Ashman et al., 2022)) but they either
frequently crash or encounter out-of-memory issues when processing large
graphs used in the experiment. These baselines are not designed for skeleton
posterior inference; however, we sum the probabilities of all edge types,
calibrate to a maximum of 1, and recast their output as a skeleton posterior
for a meaningful comparison in Table 3. In addition, we also conduct an
ablation study by replacing the ML4S-like model with an end-to-end
Transformer-based model (denoted as SPOT (w/ Ke et al. (2022))) to justify our
design consideration in Sec. 3.2.

Table 3. Evaluation on skeleton posterior inference (avg. on ten datasets).
AUROC and AUPRC: higher is better; KL: lower is better.

Metric | AVICI | FCI* | RFCI* | SPOT
---|---|---|---|---
AUROC | 0.96 | 0.97 | 0.95 | 0.99
AUPRC | 0.83 | 0.92 | 0.91 | 0.97
KL | 0.05 | 0.06 | 0.10 | 0.03

Comparison with Baselines. Overall, we observe that SPOT consistently
outperforms the baselines across all criteria in Table 3. This finding is
interpreted as encouraging and reasonable. First, SPOT employs a machine
learning model to estimate the adjacency probability while FCI and RFCI apply
hard thresholds over p-values to reject edges. Since the relationship between
p-values and adjacency probabilities may be intricate, SPOT is more adaptable
in capturing these dependencies. Second, despite having a similar focus on
variational inference, AVICI is not intended for use with causally
insufficient data, and the size of the graph may make it more challenging to
perform whole-graph inference. In contrast, SPOT efficiently resolves this
problem by edge decomposition.

Ablation Study. We also observe that the Transformer-based model (SPOT (w/ Ke
et al. (2022))) yields a lower AUROC (0.95) and AUPRC (0.77) as well as a
higher KL (0.07) than the ML4S-like cascade model (i.e., worse performance on
all criteria). This finding is consistent with our design consideration in
Sec. 3.2 that the cascade model is more suitable for the task of skeleton
posterior inference. We presume that the cascade model effectively decompose
the whole graph into edge predictions and reduce the complexity of the task,
which is crucial for large graphs.

Out-of-Distribution Setting. In addition, we also explore the effectiveness of
the domain adaption procedure in SPOT for out-of-distribution datasets. To do
so, we generate ten ADMGs from the Scale-Free (SF) model as test datasets and
evaluate the enhancement of our domain adaption strategy. On the SF graphs,
the original model trained on ER graphs suffers from a downgrade on KL (from
0.03 to 0.06). When augmented with the lightweight domain adaption step, the
KL is further reduced to 0.05 (16.7% enhancement).

###  4.3. Extension to Other Settings

Though the primary focus of SPOT lies in linear Gaussian data with latent
confounders, we also explore its effectiveness on other types of data. In
particular, we integrate SPOT with neural variational ADMG learning (Ashman et
al., 2022) (N-ADMG) which handles non-linear data and neural variational DAG
learning (Geffner et al., 2022) (DECI) and GFlowNets (Deleu et al., 2022) on
linear Gaussian data without latent confounders.

Table 4. Integration with neural variational ADMG learning (avg. on 100
sampled ADMGs).

| SHD | FDR | TPR
---|---|---|---
N-ADMG (Ashman et al., 2022) | 2.61 | 0.39 | 1.0
N-ADMG + SPOT | 2.07 -21% | 0.34 -13% | 1.0

Neural Variational ADMG Learning. Due to the scalability issue of N-ADMG, we
are unable to apply it to datasets with similar sizes as in Sec. 4.1. Instead,
we use the same dataset in its original paper (Ashman et al., 2022) with five
nodes and 2000 samples under the non-linear Gaussian setting. Here, we use the
kernel-based conditional independence test (Zhang et al., 2011) to compute the
required statistics in Sec. 3.2. However, it is worth noting that SPOT itself
is agnostic to the functional form of the SCM as long as the distribution is
faithful and the required statistics can be appropriately computed. Since
N-ADMG yields a distribution over ADMGs, we sample 100 ADMGs from the
posterior and report the average performance in Table 4. Note that, in this
setting, because the ADMG is identifiable, we compare the generated ADMG with
the ground-truth ADMG. We observe that SPOT significantly improves the
performance of N-ADMG on all criteria. In particular, it reduces the FDR by
13% and the SHD by 21%. In the dataset, N-ADMG is capable of identifying all
true edges (i.e., TPR=1.0). However, it also introduces many false positives.
SPOT effectively reduces the number of false positives and improves the
overall performance. In this regard, SPOT offers a principled way to
discourage potentially superfluous edges and thus alleviates the overfitting
issue.

Table 5. Integration with neural variational DAG learning (avg. on 100 sampled
DAGs).

| Skeleton F1 | Orientation F1 | SHD
---|---|---|---
DECI (Geffner et al., 2022) | 0.55 | 0.44 | 121.2
DECI + SPOT | 0.58 +5% | 0.48 +9% | 118.7 -2%

Neural Variational DAG Learning. Since MAG is a generalization of DAG, our
framework is naturally applicable to infer the skeleton posterior of DAGs and
thus can be integrated with neural variational DAG learning. To demonstrate
this potential extension of SPOT, we integrate it with DECI (Geffner et al.,
2022) and report the result in Table 5. Here, even though SPOT is originally
designed for MAGs, we can observe a non-trivial performance boost on DAGs. In
particular, SPOT improves the F1 score of skeleton and orientation by 5% and
9%, respectively, and reduces the SHD by 2%. This finding indicates a
potential avenue for future research to further improve the performance of
differentiable DAG learning algorithms with skeleton posterior. In addition to
DECI, we also tried to integrate SPOT with GFlowNets (Deleu et al., 2022)
which estimates the posterior of DAGs using generative flow networks. With the
same setting, we observe that SPOT reduces the expected SHD of the posterior
from 339.8 to 224.4, indicating a 34% improvement. We present the full results
in the appendix due to the space limit.

##  5\. Related Work

Causal Discovery with Latent Confounders. Causal discovery with latent
confounders involves constraint-based methods (Spirtes et al., 2000; Zhang,
2008; Colombo et al., 2012; Rohekar et al., 2021; Ma et al., 2023) and score-
based methods (Tsirlis et al., 2018; Chen et al., 2021; Claassen and Bucur,
2022). Constraint-based methods establish causal graphs based on conditional
independence (Ma et al., 2024), while score-based methods find maximal
likelihood estimations. Differentiable causal discovery utilizes continuous
optimization techniques (Bhattacharya et al., 2021). With the prosperity of
its practical applications (Ji et al., 2023a, c, b, d; Wang and Drton, 2023;
Ma et al., 2022b), SPOT advances scalable and accurate causal discovery with
latent confounders for real-world use.

Posterior Inference of Causal Structure. Traditional causal discovery provides
a maximum-likelihood point estimation, but reliability is often criticized for
small sample sizes (Cundy et al., 2021). Providing a posterior distribution is
more desired using MCMC (Viinikka et al., 2020), variational inference (Cundy
et al., 2021; Lorch et al., 2021, 2022), reinforcement learning (Yang et al.,
2023; Agrawal et al., 2018), generative flow network (Deleu et al., 2022), or
supervised learning (Ma et al., 2022a). All methods estimate DAGs’ posterior
distributions and is infeasible for MAGs. Among these methods, AVICI (Lorch et
al., 2022) is the most relevant to SPOT as it also uses amortized variational
inference. In a nutshell, it trains an inference model over samples from a
known (static) simulator and learn the DAGs from data. In the context of ADMG
learning, N-ADMG (Ashman et al., 2022) presents a variational inference method
to estimate the posterior of ADMGs on non-linear data. However, it is not
scalable well to large datasets, as shown in Sec. 4. SPOT trains an inference
model on these samples from static/dynamic simulators and alleviate the
additional assumption of the data.

The Role of Skeleton in Causal Discovery. Skeleton learning is crucial for
constraint-based methods (Yu et al., 2016) and affects overall performance.
Score-based DAG learning methods like MMHC use skeletons to reduce search
space (Tsamardinos et al., 2006). Ma et al. improve the performance of NOTEARS
via a more precise skeleton learned by ML4S (Ma et al., 2022a). SPOT
demonstrates the importance of skeleton learning in differentiable causal
discovery. Taking these evidences into consideration, it may become clear that
skeleton learning can serve as the backbone for general causal discovery
tasks.

##  6\. Conclusion

In this paper, we introduce a framework for differentiable causal discovery in
the presence of latent confounders with skeleton posterior. To this end, we
propose SPOT, which features a highly efficient variational inference
algorithm to estimate the underlying skeleton distribution from given
observational data. It also provides a stochastic optimization procedure to
seamlessly incorporate the skeleton posterior with the differentiable causal
discovery pipeline in an “out-of-the-box” manner. The results of our
experiments are highly encouraging and show that SPOT outperforms all existing
methods to a notable extent.

###### Acknowledgements.

The authors would like to thank the anonymous reviewers for their valuable
comments. The authors from HKUST were supported in part by a RGC CRF grant
under the contract C6015-23G.

##  A. Preliminary

To keep the paper self-contained, we provide additional definitions and
notations used in the main text.

###### Definition 0 (Directed Cycle).

Given an ADMG $G$ with a set of nodes $V$ and a set of edges $E$, a directed
cycle exists when there is a directed path from a node $V_{i}$ back to itself.

###### Definition 0 (Almost Directed Cycle).

Given an ADMG $G$ with a set of nodes $V$ and a set of edges $E$, an almost
directed cycle exists when there is a bidirected edge $V_{i}\leftrightarrow
V_{j}$ such that $V_{i}\in\bm{An}_{G}(V_{j})$. $\bm{An}_{G}(V_{j})$ denotes
the set of ancestors (by directed paths) of $V_{j}$ in $G$.

###### Definition 0 (Bow).

Given an ADMG $G$ with a set of nodes $V$ and a set of edges $E$, a bow exists
when there are two edges $V_{i}\rightarrow V_{j}$ and $V_{j}\leftrightarrow
V_{j}$.

###### Definition 0 (Bow-free ADMG).

An ADMG is bow-free if it does not contain any directed or a bow.

A path $(X,W_{1},\cdots,W_{k},Y)$ is said to be blocked by
$Z\subseteq\bm{X}\setminus\\{X,Y\\}$ if there exists a node
$W_{i}\in\\{W_{1},\cdots,W_{k}\\}$ such that a) $W_{i}$ is not a collider but
a member of $Z$, or b) $W_{i}$ is a collider but not an ancestor of any nodes
of $Z$. We now introduce m-separation.

###### Definition 0 (m-separation (Zhang, 2008)).

$X,Y$ are m-separated by $Z$ (denoted by $X{\Perp_{\mathcal{G}}}Y\mid Z$) if
all paths between $X,Y$ are blocked by $Z$.

##  B. Implementation of SCL Models in SPOT

In accordance with ML4S (Ma et al., 2022a), we implement the supervised causal
learning (SCL) model in SPOT as a series of cascade xgboost (Chen and
Guestrin, 2016) classifiers with default hyperparameters. For linear datasets,
we use Fisher-z test to calculate the conditional independence and constitute
the features for the classifier. For non-linear datasets, we use the kernel-
based conditional independence test (Zhang et al., 2011) instead. We use the
output of the last layer of the cascade xgboost classifiers as the posterior
probability of the presence of an edge.

Figure 5. Convergence of the model in Ke et al. (2022).

Since the model used in Ke et al. (2022) is not publicly available, we made
our best effort to replicate the model based on the description in the paper.
The experiments are conducted on a server with NVIDIA A6000 GPU and 256GB RAM.
The model is trained with 128 epochs and a batch size of 1. We use the Adam
optimizer with a learning rate of 0.0003. When trained with 128 epochs, the
metrics is seen to converge, as shown in Fig. 5. On the large-scale causal
graphs with 50–100 nodes, each dataset contains 1,000 observational samples.
However, the GPU (with 48GB memory which is the largest GPU memory available
to us) can at most handle 600 samples at a time. Hence, we drop the last 400
samples in each dataset. This issue, to a certain extent, indicates the
scalability limitation of end-to-end supervised causal learning models and
promotes the necessity of ML4S-like methods, which are much more scalable and
efficient.

##  C. Ablation Study on Sparsity Prior

In our optimization procedure, we enforce a heuristic in which sparsifying
proposals are always accepted. This design attempts to analogize a range of
the regularization term in score-based causal discovery algorithms where
sparse causal graphs are preferred over dense ones. To evaluate the impact of
this design, we conduct an ablation study by removing the sparsity prior from
the optimization procedure. The results are shown in Table 6.

Table 6. Ablation study on sparsity prior.

| Skeleton F1 | Head F1 | Tail F1
---|---|---|---
ABIC | 0.84 | 0.76 | 0.67
SPOT w/o Sparsity | 0.90 | 0.81 | 0.71
SPOT | 0.91 | 0.86 | 0.78

Overall, we observed that the strategy indeed provides a non-trivial
improvement to performance and believe that such sparsity consideration plays
an important role in “pruning” dense causal graphs.

##  D. Compute Time

We refrain from providing a table regarding compute time because different
algorithms are computed on different architectures, making it difficult to
draw meaningful conclusions. For instance, FCI can only use a single core,
whereas ABIC/SPOT can utilize up to 32 cores, and N-ADMG primarily relies on
GPU. Additionally, on some large datasets, ABIC may fail to terminate even
after running for one week. These issues complicate direct comparisons.
Following your suggestion and consistent with the main setup in our paper, we
report the compute times for ABIC and SPOT on the same computing architecture
below.

Table 7. Average compute time (in seconds) for ABIC and SPOT. Algorithm | Avg. Time (sec.)
---|---
ABIC | 3139
SPOT | 2883

Overall, we believe SPOT enables lower compute time as well as better
performance.

##  E. Comparison by Different Number of Nodes

Figure 6. Comparison of methods on different node sizes.

We report the full results of the comparison of different methods on different
methods, including FCI, RFCI, ABIC and SPOT. Aligned with our observations in
Sec. 4.1, SPOT manifests superior performance in all node sizes except on
20-node datasets.

##  F. Integration with GFlowNets

Table 8. Integration with GFlowNets (Deleu et al., 2022).

| F1 | Precision | Recall
---|---|---|---
GFlowNets (Deleu et al., 2022) | 0.22 | 0.14 | 0.45
GFlowNets + SPOT | 0.25 | 0.18 | 0.44

We also explore the integration of SPOT with GFlowNets (Deleu et al., 2022) to
further improve the performance of differentiable causal discovery. GFlowNets
is a generative flow network that learns the posterior distribution of DAGs.
We use the same experimental setup as in Sec. 4.3, and the results are shown
in Table 8. We observe that SPOT improves the performance of GFlowNets by
13.6% in F1 score and 28.6% in precision at the cost of a slight decrease in
recall. This result demonstrates the synergistical effect of SPOT in
differentiable causal discovery.

##  G. Real-world Data

Figure 7. Application on Sachs dataset.

We also explore the application of our proposed framework on Sachs (Sachs et
al., 2005), a real-world dataset with 853 entries on protein expressions
involved in human immune system cells. Fig. 7 presents the promising results
of SPOT applied to the Sachs dataset. Since every edge endpoint of Sachs is
unidentifiable, its PAG corresponds exactly to the skeleton. Therefore, we
evaluate the learned graph’s quality at the skeleton level. The discovered
skeleton in Fig.7 achieves an F1 score of 0.63. Although there is a decrease
in performance compared to the results in Sec. 4.1, it still significantly
outperforms ABIC (with an F1 score of 0.48) and provides satisfactory results.
The decrease in performance can be partially attributed to the non-Gaussian
nature of real-world data, which affects the fundamental premise of ABIC, the
backbone of our framework. Since SPOT itself (i.e., skeleton posterior
inference phase) is agnostic to Gaussianity, it delivers an impressive 31.2%
improvement over vanilla ABIC. We believe that the performance of SPOT has the
potential for even further improvement by incorporating differentiable methods
that handle non-Gaussian data in future work.

## References

  * (1)
  * Agrawal et al. (2018) Raj Agrawal, Caroline Uhler, and Tamara Broderick. 2018.  Minimal I-MAP MCMC for scalable structure discovery in causal DAG models. In _International Conference on Machine Learning_. PMLR, 89–98.
  * Ambrogioni et al. (2019) Luca Ambrogioni, Umut Güçlü, Julia Berezutskaya, Eva Borne, Yaǧmur Güçlütürk, Max Hinne, Eric Maris, and Marcel Gerven. 2019.  Forward amortized inference for likelihood-free variational marginalization. In _The 22nd International Conference on Artificial Intelligence and Statistics_. PMLR, 777–786.
  * Ashman et al. (2022) Matthew Ashman, Chao Ma, Agrin Hilmkil, Joel Jennings, and Cheng Zhang. 2022.  Causal Reasoning in the Presence of Latent Confounders via Neural ADMG Learning. In _The Eleventh International Conference on Learning Representations_.
  * Bernstein et al. (2020) Daniel Bernstein, Basil Saeed, Chandler Squires, and Caroline Uhler. 2020.  Ordering-based causal structure learning in the presence of latent variables. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 4098–4108.
  * Bhattacharya et al. (2021) Rohit Bhattacharya, Tushar Nagarajan, Daniel Malinsky, and Ilya Shpitser. 2021.  Differentiable causal discovery under unmeasured confounding. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 2314–2322.
  * Bühlmann et al. (2014) Peter Bühlmann, Jonas Peters, and Jan Ernest. 2014.  CAM: Causal additive models, high-dimensional order search and penalized regression.  _Annals of Stats._ (2014).
  * Chen et al. (2021) Rui Chen, Sanjeeb Dash, and Tian Gao. 2021.  Integer programming for causal structure learning in the presence of latent variables. In _International Conference on Machine Learning_. PMLR, 1550–1560.
  * Chen and Guestrin (2016) Tianqi Chen and Carlos Guestrin. 2016.  Xgboost: A scalable tree boosting system. In _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_. 785–794.
  * Claassen and Bucur (2022) Tom Claassen and Ioan G Bucur. 2022.  Greedy equivalence search in the presence of latent confounders. In _Uncertainty in Artificial Intelligence_. PMLR, 443–452.
  * Colombo et al. (2012) Diego Colombo, Marloes H Maathuis, Markus Kalisch, and Thomas S Richardson. 2012.  Learning high-dimensional directed acyclic graphs with latent and selection variables.  _Annals of Stats._ (2012).
  * Cundy et al. (2021) Chris Cundy, Aditya Grover, and Stefano Ermon. 2021.  Bcd nets: Scalable variational approaches for bayesian causal discovery.  _Advances in Neural Information Processing Systems_ 34 (2021), 7095–7110.
  * Dai et al. (2023) Haoyue Dai, Rui Ding, Yuanyuan Jiang, Shi Han, and Dongmei Zhang. 2023.  Ml4c: Seeing causality through latent vicinity. In _Proceedings of the 2023 SIAM International Conference on Data Mining (SDM)_. SIAM, 226–234.
  * Deleu et al. (2022) Tristan Deleu, António Góis, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. 2022.  Bayesian structure learning with generative flow networks. In _Uncertainty in Artificial Intelligence_. PMLR, 518–528.
  * Drton et al. (2009) Mathias Drton, Michael Eichler, and Thomas S Richardson. 2009.  Computing Maximum Likelihood Estimates in Recursive Linear Models with Correlated Errors.  _Journal of Machine Learning Research_ 10, 10 (2009).
  * Fletcher (1987) Roger Fletcher. 1987.  Practical methods of optimization.  _A Wiley Interscience Publication_ (1987).
  * Friedman et al. (1999) Nir Friedman, Moises Goldszmidt, and Abraham Wyner. 1999.  Data analysis with bayesian networks: a bootstrap approach. In _Proceedings of the Fifteenth conference on Uncertainty in artificial intelligence_. 196–205.
  * Geffner et al. (2022) Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Agrin Hilmkil, et al. 2022\.  Deep End-to-end Causal Inference. In _NeurIPS 2022 Workshop on Causality for Real-world Impact_.
  * Ji et al. (2023b) Zhenlan Ji, Pingchuan Ma, Zongjie Li, and Shuai Wang. 2023b.  Benchmarking and Explaining Large Language Model-based Code Generation: A Causality-Centric Approach.  _arXiv preprint arXiv:2310.06680_ (2023).
  * Ji et al. (2023a) Zhenlan Ji, Pingchuan Ma, and Shuai Wang. 2023a.  Perfce: Performance debugging on databases with chaos engineering-enhanced causality analysis. In _2023 38th IEEE/ACM International Conference on Automated Software Engineering (ASE)_. IEEE, 1454–1466.
  * Ji et al. (2023c) Zhenlan Ji, Pingchuan Ma, Shuai Wang, and Yanhui Li. 2023c.  Causality-aided trade-off analysis for machine learning fairness. In _2023 38th IEEE/ACM International Conference on Automated Software Engineering (ASE)_. IEEE, 371–383.
  * Ji et al. (2023d) Zhenlan Ji, Pingchuan Ma, Yuanyuan Yuan, and Shuai Wang. 2023d.  Cc: Causality-aware coverage criterion for deep neural networks. In _2023 IEEE/ACM 45th International Conference on Software Engineering (ICSE)_. IEEE, 1788–1800.
  * Ke et al. (2022) Nan Rosemary Ke, Silvia Chiappa, Jane X Wang, Jorg Bornschein, Anirudh Goyal, Melanie Rey, Theophane Weber, Matthew Botvinick, Michael Curtis Mozer, and Danilo Jimenez Rezende. 2022.  Learning to Induce Causal Structure. In _International Conference on Learning Representations_.
  * Kirkpatrick et al. (1983) Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. 1983.  Optimization by simulated annealing.  _Science_ (1983).
  * Lorch et al. (2021) Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. 2021.  Dibs: Differentiable bayesian structure learning.  _Advances in Neural Information Processing Systems_ 34 (2021), 24111–24123.
  * Lorch et al. (2022) Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. 2022.  Amortized Inference for Causal Structure Learning.  _arXiv:2205.12934_ (2022).
  * Ma et al. (2022a) Pingchuan Ma, Rui Ding, Haoyue Dai, Yuanyuan Jiang, Shuai Wang, Shi Han, and Dongmei Zhang. 2022a.  Ml4s: Learning causal skeleton from vicinal graphs. In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_. 1213–1223.
  * Ma et al. (2023) Pingchuan Ma, Rui Ding, Shuai Wang, Shi Han, and Dongmei Zhang. 2023.  XInsight: eXplainable Data Analysis Through The Lens of Causality.  _Proceedings of the ACM on Management of Data_ 1, 2 (2023), 1–27.
  * Ma et al. (2022b) Pingchuan Ma, Zhenlan Ji, Qi Pang, and Shuai Wang. 2022b.  Noleaks: Differentially private causal discovery under functional causal model.  _IEEE Transactions on Information Forensics and Security_ 17 (2022), 2324–2338.
  * Ma et al. (2024) Pingchuan Ma, Zhenlan Ji, Peisen Yao, Shuai Wang, and Kui Ren. 2024.  Enabling Runtime Verification of Causal Discovery Algorithms with Automated Conditional Independence Reasoning. In _Proceedings of the 46th IEEE/ACM International Conference on Software Engineering_. 1–13.
  * Microsoft (2024) Microsoft. 2024.  reliableAI.  https://github.com/microsoft/reliableAI.
  * Rohekar et al. (2021) Raanan Y Rohekar, Shami Nisimov, Yaniv Gurwicz, and Gal Novik. 2021.  Iterative causal discovery in the possible presence of latent confounders and selection bias.  _Advances in Neural Information Processing Systems_ 34 (2021), 2454–2465.
  * Sachs et al. (2005) Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan. 2005.  Causal protein-signaling networks derived from multiparameter single-cell data.  _Science_ (2005).
  * Salehkaleybar et al. (2020) Saber Salehkaleybar, AmirEmad Ghassami, Negar Kiyavash, and Kun Zhang. 2020.  Learning linear non-Gaussian causal models in the presence of latent variables.  _The Journal of Machine Learning Research_ 21, 1 (2020), 1436–1459.
  * Spirtes et al. (2000) Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. 2000.  _Causation, prediction, and search_.  MIT press.
  * Tsamardinos et al. (2006) Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. 2006.  The max-min hill-climbing Bayesian network structure learning algorithm.  _Machine learning_ 65 (2006), 31–78.
  * Tsirlis et al. (2018) Konstantinos Tsirlis, Vincenzo Lagani, Sofia Triantafillou, and Ioannis Tsamardinos. 2018.  On scoring maximal ancestral graphs with the max–min hill climbing algorithm.  _International Journal of Approximate Reasoning_ 102 (2018), 74–85.
  * Viinikka et al. (2020) Jussi Viinikka, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. 2020.  Towards scalable bayesian learning of causal dags.  _Advances in Neural Information Processing Systems_ 33 (2020), 6584–6594.
  * Vowels et al. (2022) Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden. 2022.  D’ya like dags? a survey on structure learning and causal discovery.  _Comput. Surveys_ 55, 4 (2022), 1–36.
  * Wang and Drton (2023) Y Samuel Wang and Mathias Drton. 2023.  Causal discovery with unobserved confounding and non-Gaussian data.  _Journal of Machine Learning Research_ 24, 271 (2023), 1–61.
  * Yang et al. (2023) Dezhi Yang, Guoxian Yu, Jun Wang, Zhengtian Wu, and Maozu Guo. 2023.  Reinforcement causal structure learning on order graph. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 37. 10737–10744.
  * Yu et al. (2016) Kui Yu, Jiuyong Li, and Lin Liu. 2016.  A review on algorithms for constraint-based causal discovery.  _arXiv preprint arXiv:1611.03977_ (2016).
  * Zhang (2008) Jiji Zhang. 2008.  On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias.  _Artificial Intelligence_ 172, 16-17 (2008), 1873–1896.
  * Zhang et al. (2011) Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. 2011.  Kernel-based conditional independence test and application in causal discovery. In _Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence_. 804–813.
  * Zheng et al. (2018) Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. 2018.  Dags with no tears: Continuous optimization for structure learning.  _Advances in neural information processing systems_ 31 (2018).
  * Zhu et al. (2019) Shengyu Zhu, Ignavier Ng, and Zhitang Chen. 2019.  Causal Discovery with Reinforcement Learning. In _International Conference on Learning Representations_.