# A Survey on Causal Discovery Methods for I.I.D. and Time Series Data

Uzma Hasan<EMAIL_ADDRESS>
Causal AI Lab
Department of Information Systems
University of Maryland, Baltimore County
Baltimore, MD, USA Emam Hossain<EMAIL_ADDRESS>
Causal AI Lab
Department of Information Systems
University of Maryland, Baltimore County
Baltimore, MD, USA Md Osman Gani<EMAIL_ADDRESS>
Causal AI Lab
Department of Information Systems
University of Maryland, Baltimore County
Baltimore, MD, USA

###### Abstract

The ability to understand causality from data is one of the major milestones
of human-level intelligence. Causal Discovery (CD) algorithms can identify the
cause-effect relationships among the variables of a system from related
observational data with certain assumptions. Over the years, several methods
have been developed primarily based on the statistical properties of data to
uncover the underlying causal mechanism. In this study, we present an
extensive discussion on the methods designed to perform causal discovery from
both independent and identically distributed (I.I.D.) data and time series
data. For this purpose, we first introduce the common terminologies used in
causal discovery literature and then provide a comprehensive discussion of the
algorithms designed to identify causal relations in different settings. We
further discuss some of the benchmark datasets available for evaluating the
algorithmic performance, off-the-shelf tools or software packages to perform
causal discovery readily, and the common metrics used to evaluate these
methods. We also evaluate some widely used causal discovery algorithms on
multiple benchmark datasets and compare their performances. Finally, we
conclude by discussing the research challenges and the applications of causal
discovery algorithms in multiple areas of interest.

##  1 Introduction

The identification of the cause-effect relationships among the variables of a
system from the corresponding data is called Causal Discovery (CD). A major
part of the causal analysis involves unfolding the _cause and effect
relationships_ among the entities in complex systems that can help us build
better solutions in health care, earth science, politics, business, education,
and many other diverse areas (Peyrot (1996), Nogueira et al. (2021)). The
_causal explanations_ precisely the causal factors obtained from a causal
analysis play an important role in decision-making and policy formulation as
well as in foreseeing the consequences of interventions without actually doing
them. Causal discovery algorithms enable the _discovery of the underlying
causal structure_ given a set of observations. The underlying causal structure
also known as a causal graph (CG) is a representation of the cause-effect
relationships between the variables in the data (Pearl (2009)). Causal graphs
represent the causal relationships with directed arrows from the cause to the
effect.

Figure 1: Causal Discovery: Identification of a causal graph from data.

Discovering the causal relations, and thereby, the estimation of their effects
would enable us to understand the underlying data generating mechanism (DGM)
better, and take necessary interventional actions. However, traditional
Artificial Intelligence (AI) applications rely solely on predictive models and
often ignore causal knowledge. Systems without the knowledge of causal
relationships often cannot make rational and informed decisions (Marwala
(2015)). The result may be devastating when correlations are mistaken for
causation. Because two variables can be highly correlated, and yet not have
any causal influence on each other. There may be a third variable often called
a latent confounder or hidden factor that may be causing both of them (see
Figure 2 (a)). Thus, _embedding the knowledge of causal relationships_ in
black-box AI systems is important to improve their explainability and
reliability (Dubois & Prade (2020), Ganguly et al. (2023)). In multiple fields
such as healthcare, politics, economics, climate science, business, and
education, the ability to understand causal relations can facilitate the
formulation of better policies with a greater understanding of the data.

(a) (b)

Figure 2: (a) Latent confounder $L$ causes both variables $S$ and $C$, and the
association between $S$ and $C$ is denoted by ? which can be mistaken as
causation. The graph in (b) is a causal graph depicting the causes and effects
of cancer (Korb & Nicholson (2010)).

The standard approach to discover the cause-effect relationships is to perform
randomized control trials (RCTs) (Sibbald & Roland (1998)). However, RCTs are
often infeasible to conduct due to high costs and some ethical reasons (Resnik
(2008)). As a result, over the last few decades, researchers have developed a
variety of methods to unravel causal relations from purely observational data
(Glymour et al. (2019), Vowels et al. (2021)). These methods are often based
on some assumptions about the data and the underlying mechanism. The _outcome_
of any causal discovery method is a causal graph or a causal adjacency matrix
where the cause and effect relations among the entities or variables are
represented. The structure of a causal graph is often similar to a _directed
acyclic graph (DAG)_ where directed edges from one variable to another
represent the cause-effect relationship between them. Figure 2 (b) represents
a causal graph showing the factors that are responsible for causing Cancer.
This type of structural representation of the underlying data-generating
mechanism is beneficial for understanding how the system entities interact
with each other.

There exists a wide range of approaches for performing causal discovery under
different settings or assumptions. Some approaches are designed particularly
for _independent and identically distributed (I.I.D.) data_ (Spirtes et al.
(2000b), Chickering (2002)) i.e. non-temporal data while others are focused on
_time series data_ (Runge et al. (2019), Hyvärinen et al. (2010)) or temporal
data. Since in real-world settings, both types of data are available in
different problem domains, it is essential to have approaches to perform
causal structure recovery from both of these. Recently, there has been a
growing body of research that considers _prior knowledge incorporation_ for
recovering the causal relationships (Mooij et al. (2020), Hasan & Gani (2022),
Hasan & Gani (2023)). Although there exist some surveys (see Table 1) on
causal discovery approaches (Heinze-Deml et al. (2018), Glymour et al. (2019),
Guo et al. (2020), Vowels et al. (2021), Assaad et al. (2022b)), none of these
present a comprehensive review of the different approaches designed for
structure recovery from both I.I.D. and time series data. Also, these surveys
do not discuss the approaches that perform causal discovery in the presence of
background knowledge. Hence, the goal of this survey is to provide an overview
of the wide range of existing approaches for performing causal discovery from
I.I.D. as well as time series data under different settings. Existing surveys
lack a combined overview of the approaches present for both I.I.D. and time
series data. So in this survey, we want to introduce the readers to the
methods available in both domains. We discuss prominent methods based on the
different approaches such as conditional independence (CI) testing, score
function usage, functional causal models (FCMs), continuous optimization
strategy, prior knowledge infusion, and miscellaneous ones. These methods
primarily differ from each other based on the primary strategy they follow.
Apart from introducing the different causal discovery approaches and
algorithms for I.I.D. and time series data, we also discuss the different
tools, metrics, and benchmark datasets used for performing CD and the
challenges and applications of CD in a wide range of areas.

Table 1: Comparison among the existing surveys for causal discovery approaches. A discussion on the different approaches can be found in section 3 and section 4. Survey |  Focused Approaches | I.I.D. Data | Time Series Data
---|---|---|---
Heinze-Deml et al. (2018) |  Constraint, Score, Hybrid & FCM-based approaches. | $\checkmark$ | $\times$
Glymour et al. (2019) |  Traditional Constraint-based, Score-based, & FCM-based approaches. | $\checkmark$ | $\times$
Guo et al. (2020) |  Constraint-based, Score-based, & FCM-based approaches. | $\checkmark$ | $\times$
Vowels et al. (2021) |  Continuous Optimization-based. | $\checkmark$ | $\times$
Assaad et al. (2022b) |  Constraint-based, Score-based, FCM-based, etc. approaches for time series data. | $\times$ | $\checkmark$
This study |  Constraint-based, Score-based, FCM-based, Hybrid-based, Continuous-Optimization-based, Prior-Knowledge-based, and Miscellaneous. | $\checkmark$ | $\checkmark$

To summarize, the structure of this paper is as follows: First, we provide a
brief introduction to the common terminologies in the field of causal
discovery (section 2). Second, we discuss the wide range of causal discovery
approaches that exist for both I.I.D. (section 3) and time-series data
(section 4). Third, we briefly overview the common evaluation metrics (section
5) and datasets (section 6) used for evaluating the causal discovery
approaches, and report the performance comparison of some causal discovery
approaches in section 7. Fourth, we list the different technologies and open-
source software (section 8) available for performing causal discovery. Fifth,
we discuss the challenges (section 9.1) and applications (section 9.2) of
causal discovery in multiple areas such as healthcare, business, social
science, economics, and so on. Lastly, we conclude by discussing the scopes of
improvement in future causal discovery research, and the importance of
causality in improving the existing predictive AI systems which can thereby
impact informed and reliable decision-making in different areas of interest
(section 10).

##  2 Preliminaries of Causal Discovery

In this section, we briefly discuss the important terminologies and concepts
that are widely used in causal discovery. Some common notations used to
explain the terminologies are presented in Table 2.

Table 2: Common notations. Notation | Description
---|---
$G$ | A graph or DAG or ground-truth graph
$G^{\prime}$ | An estimated graph
$X,Y,Z,W$ | Observational variables
$X$ — $Y$ |  An unoriented or undirected edge between $X$ and $Y$
$X$ → $Y$ |  A directed edge from $X$ to $Y$ where $X$ is the cause and Y is the effect
$X$ $\not\to$ $Y$ |  Absence of an edge or causal link between $X$ and $Y$
$X$ → $Z$ ← $Y$ |  V-structure or Collider where $Z$ is the common child of $X$ and $Y$
$\perp\\!\\!\\!\perp$ | Independence or d-separation
$X$ $\perp\\!\\!\\!\perp$ $Y$ $|$ $Z$ |  $X$ is d-separated from $Y$ given $Z$

###  2.1 Graphical Models

A graph G = (V, E) consists of a set of vertices (nodes) V and a set of edges
E where the edges represent the relationships among the vertices. Figure 3 (a)
represents a graph $G$ with vertices $V=[X,Y,Z]$ and edges
$E=[(X,Y),(X,Z),(Z,Y)]$. There can be different types of edges in a graph such
as directed edges (→), undirected edges (-), bi-directed edges
($\leftrightarrow$), etc. (Colombo et al. (2012)). A graph that consists of
only undirected edges (-) between the nodes which represent their adjacencies
is called a _skeleton graph_ $S_{G}$. This type of graph is also known as an
_undirected graph_ (Figure 3 (b)). A graph that has a mixture of different
types of edges is known as a _mixed graph_ $M_{G}$ (Figure 3 (c)). A _path_
$p$ between two nodes $X$ and $Y$ is a sequence of edges beginning from $X$
and ending at $Y$. A _cycle_ $c$ is a path that begins and ends at the same
vertex. A graph with no cycle $c$ is called an _acyclic graph_. And, a
directed graph in which the edges have directions (→) and has no cycle is
called a directed acyclic graph (DAG). In a DAG $G$, a directed path from $X$
to $Y$ implies that $X$ is an ancestor of $Y$, and $Y$ is a descendant of $X$.
The graph $G$ in Figure 3 (a) is a DAG as it is acyclic, and consists of
directed edges.

Figure 3: (a) A graph $G$, (b) its _skeleton_ graph $S_{G}$, (c) a _mixed
graph_ $M_{G}$ with directed & undirected edges.

There can be different kinds of DAGs based on the type of edges they contain.
A class of DAG known as partially directed acyclic graph (PDAG) contains both
directed ($\rightarrow$) and undirected (-) edges. The mixed graph of Figure 3
(c) is also a PDAG. A completed PDAG (CPDAG) consists of directed
($\rightarrow$) edges that exist in every DAG $G$ having the same conditional
dependencies, and undirected (-) edges that are reversible in $G$. An
extension of DAGs that retain many of the significant properties that are
associated with DAGs is known as ancestral graphs (AGs). Two different DAGs
may lead to the same ancestral graph (Richardson & Spirtes (2002a)). Often
there are hidden confounders and selection biases in real-world data.
Ancestral graphs can represent the data-generating mechanisms that may involve
latent confounders and/or selection bias, without explicitly modeling the
unobserved variables. There exist different types of ancestral graphs. A
maximal ancestral graph (MAG) is a mixed graph that can have both directed (→)
and bidirectional ($\leftrightarrow$) edges (Richardson & Spirtes (2002b)). A
partial ancestral graph (PAG) can have four types of edges such as directed
(→), bi-directed ($\leftrightarrow$), partially directed (o→), and undirected
($-$) (Spirtes (2001)). That is, edges in a PAG can have three kinds of
endpoints: $-$, o, or $>$. An ancestral graph without bi-directed edges
($\leftrightarrow$) is a DAG (Triantafillou & Tsamardinos (2016)).

###  2.2 Causal Graphical Models

A causal graphical model (CGM) or causal graph (CG) is a DAG $G$ that
represents a joint probability distribution $P$ over a set of random variables
$X=(X_{1},X_{2},…,X_{d})$ where $P$ is Markovian with respect to $G$. In a
CGM, the nodes represent variables $X$, and the arrows represent causal
relationships between them. The joint distribution $P$ can be factorized as
follows where $pa(x_{i},G)$ denotes the parents of $x_{i}$ in $G$.

$P(x_{1},…,x_{d})=\prod_{i=1}^{d}P(x_{i}|pa(x_{i},G))$ (1)

Causal graphs are often used to study the underlying data-generating mechanism
in real-world problems. For any dataset $D$ with variables $X$, causal graphs
can encode the cause-effect relationships among the variables using directed
edges ($\rightarrow$) from cause to the effect. Most of the time causal graphs
take the form of a DAG. In Figure 3 (a), $X$ is the cause that effects both
$Y$ and $Z$ (i.e. $Y$ $\leftarrow$ $X$ $\rightarrow$ $Z$). Also, $Z$ is a
cause of $Y$ (i.e. $Z$ $\rightarrow$ $Y$). The mechanism that enables the
estimation of a causal graph $G$ from a dataset $D$ is called _causal
discovery (CD)_ (Figure 1). The outcome of any causal discovery algorithm is a
causal graph $G$ where the directed edges ($\rightarrow$) represent the cause-
and-effect relationship between the variables $X$ in $D$. However, some
approaches have different forms of graphs (PDAGs, CPDAGs, ancestral graphs,
etc.) as the output causal graph. Table 3 lists the output causal graphs of
some common approaches which are discussed in section 3.

Table 3: List of some CD algorithms with their output causal graphs. A detailed discussion of the algorithms is in section 3. The cells with $\checkmark$ represent the type of graph produced by the corresponding algorithm. Algorithms | DAG | PDAG | CPDAG | MAG | PAG
---|---|---|---|---|---
PC |  |  | $\checkmark$ |  |
FCI |  |  |  |  | $\checkmark$
RFCI |  |  |  |  | $\checkmark$
GES |  |  | $\checkmark$ |  |
GIES |  | $\checkmark$ |  |  |
MMHC | $\checkmark$ |  |  |  |
LiNGAM | $\checkmark$ |  |  |  |
NOTEARS | $\checkmark$ |  |  |  |
GSMAG |  |  |  | $\checkmark$ |

####  2.2.1 Key Structures in Causal Graphs

There are three fundamental building blocks (key structures) commonly observed
in the graphical models or causal graphs, namely, _Chain_ , _Fork_ , and
_Collider_. Any graphical model consisting of at least three variables is
composed of these key structures. We discuss these basic building blocks and
their implications in dependency relationships below.

######  Definition 1 (Chain)

A chain $X\rightarrow Y\rightarrow Z$ is a graphical structure or a
configuration of three variables $X$, $Y$, and $Z$ in graph $G$ where $X$ has
a directed edge to $Y$ and $Y$ has a directed edge to $Z$ (see Figure 4 (a)).
Here, $X$ causes $Y$ and $Y$ causes $Z$, and $Y$ is called a mediator.

######  Definition 2 (Fork)

A fork $Y\leftarrow X\rightarrow Z$ is a triple of variables $X$, $Y$, and $Z$
where one variable is the common parent of the other two variables. In Figure
4 (b), the triple ($X$, $Y$, $Z$) is a fork where $X$ is a common parent of
$Y$ and $Z$.

######  Definition 3 (Collider/V-structure)

A v-structure or collider $X\rightarrow Z\leftarrow Y$ is a triple of
variables $X$, $Y$, and $Z$ where one variable is a common child of the other
two variables which are non-adjacent. In Figure 4 (c), the triple ($X$, $Y$,
$Z$) is a v-structure where $Z$ is a common child of $X$ and $Y$, but $X$ and
$Y$ are non-adjacent in the graph. Figure 4 (d) is also a collider with a
descendant $W$.

Figure 4: Fundamental building blocks in causal graphical models.

####  2.2.2 Conditional Independence in Causal Graphs

Testing for conditional independence (CI) between the variables is one of the
most important techniques to find the causal relationships among the
variables. Conditional independence between two variables $X$ and $Y$ results
when they are independent of each other given a third variable $Z$ (i.e. $X$
$\perp\\!\\!\\!\perp$ $Y$ $|$ $Z$). In the case of causal discovery, CI
testing allows deciding if any two variables are causally connected or
disconnected. An important criterion for CI testing is the d-separation
criterion which is formally defined below.

######  Definition 4 (d-separation)

(Pearl (1988)) A path $p$ in $G$ is blocked by a set of nodes $N$ if either

  1. i.

$p$ contains a chain of nodes $X\rightarrow Y\rightarrow Z$ or a fork
$X\leftarrow Y\rightarrow Z$ such that the middle node $Y$ is in $N$,

  2. ii.

$p$ contains a collider $X\rightarrow Y\leftarrow Z$ such that the collision
node $Y$ is not in $N$, and no descendant of $Y$ is in $N$.

If $N$ blocks every path between two nodes, then they are d-separated,
conditional on N, and thus are independent conditional on N.

In d-separation, _d_ stands for directional. The d-separation criterion
provides a set of rules to check if two variables are independent when
conditioned on a set of variables. The conditioning variable can be a single
variable or a set of variables. However, two variables with a directed edge
($\rightarrow$) between them are always dependent. The set of testable
implications provided by _d-separation_ can be benchmarked with the available
data $D$. If a graph $G$ might have been generated from a dataset $D$, then
_d-separation_ tells us which variables in $G$ must be independent conditional
on other variables. If every _d-separation_ condition matches a conditional
independence in data, then no further test can refute the model (Pearl
(1988)). If there is at least one path between two variables that is
unblocked, then they are _d-connected_. If two variables are _d-connected_ ,
then they are most likely dependent (except intransitive cases) (Pearl
(1988)). The d-separation or conditional independence between the variables in
the key structures (Figure 4) or building blocks of causal graphs follow some
rules which are discussed below:

  1. i.

_Conditional Independence in Chains:_ If there is only one unidirectional path
between variables $X$ and $Z$ (Figure 4 (a)), and $Y$ is any variable or set
of variables that intercept that path, then $X$ and $Z$ are conditionally
independent given $Y$, i.e. $X$ $\perp\\!\\!\\!\perp$ $Z$ $|$ $Y$.

  2. ii.

_Conditional Independence in Forks:_ If a variable $X$ is a common cause of
variables $Y$ and $Z$, and there is only one path between $Y$ and $Z$, then
$Y$ and $Z$ are independent conditional on $X$ (i.e. $Y$ $\perp\\!\\!\\!\perp$
$Z$ $|$ $X$) (Figure 4(b)).

  3. iii.

_Conditional Independence in Colliders:_ If a variable $Z$ is the collision
node between two variables $X$ and $Y$ (Figure 4(c)), and there is only one
path between $X$ and $Y$, then $X$ and $Y$ are unconditionally independent
(i.e. $X$ $\perp\\!\\!\\!\perp$ $Y$). But, they become dependent when
conditioned on $Z$ or any descendants of $Z$ (Figure 4(d)).

####  2.2.3 Markov Equivalence in Causal Graphs

A set of causal graphs having the same set of conditional independencies is
known as a Markov equivalence class (MEC). Two DAGs that are Markov equivalent
have the (i) same skeleton (the underlying undirected graph) and (ii)  same
v-structures (colliders)  (Verma & Pearl (2022)). That is, all DAGs in a MEC
share the same edges, regardless of the direction of those edges, and the same
colliders whose parents are not adjacent. Chain and Fork share the same
independencies, hence, they belong to the same MEC (Figure 5).

######  Definition 5 (Markov Blanket)

For any variable X, its Markov blanket (MB) is the set of variables such that
X is independent of all other variables given MB. The members in the Markov
blanket of any variable will include all of its parents, children, and
spouses.

Figure 5: Markov Equivalence in Chains and Fork.

Markov equivalence in different types of DAGs may vary. A partial DAG (PDAG)
a.k.a an essential graph (Perković et al. (2017)) can represent an equivalence
class of DAGs. Each equivalent class of DAGs can be uniquely represented by a
PDAG. A completed PDAG or CPDAG represents the union (over the set of edges)
of Markov equivalent DAGs, and can uniquely represent an MEC (Malinsky &
Spirtes (2016b)). More specifically, in a CPDAG, an undirected edge between
any two nodes $X$ and $Y$ indicates that some DAG in the equivalence class
contains the edge $X$→$Y$ and some DAG may contain $Y$→$X$. Figure 6 shows a
CPDAG and the DAGs ($G$ and $H$) belonging to an equivalence class.

Markov equivalence in the case of ancestral graphs works as follows. A maximal
ancestral graph (MAG) represents a DAG where all hidden variables are
marginalized out and preserves all conditional independence relations among
the variables that are true in the underlying DAG. That is, MAGs can model
causality and conditional independencies in causally insufficient systems
(Triantafillou & Tsamardinos (2016)). Partial ancestral graphs (PAGs)
represent an equivalence class of MAGs where all common edge marks shared by
all members in the class are displayed, and also, circles for those marks that
are uncommon are presented. PAGs represent all of the observed d-separation
relations in a DAG. Different PAGs that represent distinct equivalence classes
of MAGs involve different sets of conditional independence constraints. An MEC
of MAGs can be represented by a PAG (Malinsky & Spirtes (2016b)).

Figure 6: DAGs $G$ and $H$ belong to the same MEC. The leftmost graph is a
CPDAG of $G$ and $H$ with an undirected edge ($-$) between $X$ and $Z$, and
the rest of the edges same as in $G$ and $H$.

###  2.3 Structural Causal Models

Pearl (2009) defined a _class of models_ for formalizing structural knowledge
about the _data-generating process_ known as the structural causal models
(SCMs). The SCMs are valuable tools for reasoning and decision-making in
causal analysis since they are capable of representing the underlying causal
story of data (Kaddour et al. (2022)).

######  Definition 6 (Structural Causal Model)

Pearl (2009); A structural causal model is a 4-tuple $M=\langle
U,V,F,P(u)\rangle$, where

  1. i.

$U$ is a set of background variables (also called exogenous) that are
determined by factors outside the model.

  2. ii.

$V$ is a set $\\{V_{1},V_{2},\ldots,V_{n}\\}$ of endogenous variables that are
determined by variables in the model, viz. variables in $U\cup V$.

  3. iii.

$F$ is a set of functions $\\{f_{1},f_{2},\ldots,f_{n}\\}$ such that each
$f_{i}$ is a mapping from the respective domains of $U_{i}\cup{PA}_{i}$ to
$V_{i}$ and the entire set $F$ forms a mapping from $U$ to $V$. In other
words, each $f_{i}$ assigns a value to the corresponding $V_{i}\in V$,
$v_{i}\leftarrow f_{i}(pa_{i},u_{i}),$ for $i=1,2,\ldots n$.

  4. iv.

$P(u)$ is a probability function defined over the domain of $U$.

Each SCM $M$ is associated with a causal graphical model $G$ that is a DAG,
and a set of functions $f_{i}$. Causation in SCMs can be interpreted as
follows: a variable $Y$ is directly caused by $X$ if $X$ is in the function
$f$ of $Y$. In other words, each $f_{i}$ assigns a value to the corresponding
$V_{i}\in V$, $v_{i}\leftarrow f_{i}(pa_{i},u_{i}),$ for $i=1,2,\ldots n$. In
the SCM of Figure 7, $X$ is a direct cause of $Y$ as $X$ appears in the
function that assigns $Y$’s value. That is, if a variable $Y$ is the child of
another variable $X$, then $X$ is a direct cause of $Y$. In Figure 7, $U_{X}$,
$U_{Y}$ and $U_{Z}$ are the exogenous variables; $X$, $Y$ and $Z$ are the
endogenous variables, and $f_{X}$, $f_{Y}$ & $f_{Z}$ are the functions that
assign values to the variables in the system. Any variable is an exogenous
variable if $(i)$ it is an unobserved or unmeasured variable and $(ii)$ it
cannot be a descendant of any other variables. Every endogenous variable is a
descendant of at least one exogenous variable.

Figure 7: A Structural Causal Model (SCM) with causal graph $G$ and functions
$f_{X}$, $f_{Y}$ and $f_{Z}$ which denotes how the variables $X$, $Y$, and $Z$
are generated respectively.

###  2.4 Causal Assumptions

Often, the available data provide only partial information about the
underlying causal story. Hence, it is essential to make some assumptions about
the world for performing causal discovery (Lee & Honavar (2020)). Following
are the common assumptions usually made by causal discovery algorithms.

  1. i.

_Causal Markov Condition (CMC):_ The causal Markov assumption states that a
variable $X$ is independent of every other variable (except its descendants)
conditional on all of its direct causes (Scheines (1997)). That is, the CMC
requires that every variable in the causal graph is independent of its non-
descendants conditional on its parents (Malinsky & Spirtes (2016a)). In Figure
8, $W$ is the only descendant of $X$. As per the CMC, $X$ is independent of
$Z$ conditioned on its parent $Y$ ($X$ $\perp\\!\\!\\!\perp$ $Z$ $|$ $Y$).

Figure 8: Illustration of the causal Markov condition (CMC) among four
variables.

  2. ii.

_Causal Faithfulness Condition (CFC):_ The faithfulness assumption states that
except for the variables that are d-separated in a DAG, all other variables
are dependent. More specifically, for a set of variables $V$ whose causal
structure is represented by a DAG $G$, no conditional independence holds
unless entailed by the causal Markov condition (Ramsey et al. (2012)). That
is, the CFC a.k.a the Stability condition is a converse principle of the CMC.
CFC can be also explained in terms of d-separation as follows: For every three
disjoint sets of variables $X$, $Y$, and $Z$, if $X$ and $Y$ are not
d-separated by $Z$ in the causal DAG, then $X$ and $Y$ are not independent
conditioned on $Z$ (Ramsey et al. (2012)). The faithfulness assumption may
fail in certain scenarios. For example, it fails whenever there exist two
paths with equal and opposite effects between variables. It also fails in
systems with deterministic relationships among variables, and also, when there
is a failure of transitivity along a single path (Weinberger (2018)).

  3. iii.

_Causal Sufficiency:_ The causal sufficiency assumption states that there
exist no latent/hidden/unobserved confounders, and all the common causes are
measured. Thus, the assumption of causal sufficiency is satisfied only when
all the common causes of the measured variables are measured. This is a strong
assumption as it restricts the search space of all possible DAGs that may be
inferred. However, real-world datasets may have hidden confounders which might
frequently cause the assumption to be violated in such scenarios. Algorithms
that violate the causal sufficiency assumption may observe degradation in
their performance. The causal insufficiency in real-world datasets may be
overcome by leveraging domain knowledge in the discovery pipeline. The CMC
tends to fail for a causally insufficient set of variables.

  4. iv.

_Acyclicity:_ It is the most common assumption which states that there are no
cycles in a causal graph. That is, a graph needs to be acyclic in order to be
a causal graph. As per the acyclicity condition, there can be no directed
paths starting from a node and ending back to itself. This resembles the
structure of a directed acyclic graph (DAG). A recent approach (Zheng et al.
(2018)) has formulated a new function (Equation 2) to enforce the acyclicity
constraint during causal discovery in continuous optimization settings. The
weighted adjacency matrix $W$ in Equation 2 is a DAG if it satisfies the
following condition where $\circ$ is the Hadamard product, $e^{W\circ W}$ is
the matrix exponential of $W\circ W$, and $d$ is the total number of vertices.

$h(W)=tr(e^{W\circ W})-d=0$ (2)

  5. v.

_Data Assumptions:_ There can be different types of assumptions about the
data. Data may have linear or nonlinear dependencies and can be continuously
valued or discrete valued in nature. Data can be independent and identically
distributed (I.I.D.) or the data distribution may shift with time (e.g. time-
series data). Also, the data may belong to different noise distributions such
as Gaussian, Gumbel, or Exponential noise. Occasionally, some other data
assumptions such as the existence of selection bias, missing variables, hidden
confounders, etc. are found. However, in this survey, we do not focus much on
the methods with these assumptions.

forked edges, for tree=draw,align=left,l=1.7cm,edge=-latex [Causal Discovery
Algorithms
(for I.I.D. data) [Contraint-based [PC (3.1.1)
FCI (3.1.2)
Anytime FCI (3.1.3)
RFCI (3.1.4)
FCI with TBK∗ (3.1.5)
PC-stable (3.1.6)
PKCL∗ (3.1.7)
] ] [Score-based [GES (3.2.1)
FGS (3.2.2)
SGES (3.2.3)
RL-BIC (3.2.4)
A-star search (3.2.5)
Triplet A-star (3.2.6)
KCRL∗ (3.2.7)] ] [FCM-based [LiNGAM (3.3.1)
ANM (3.3.2)
PNL (3.3.3)
DirectLiNGAM∗ (3.3.4)
SAM (3.3.5)
CGNN (3.3.6)
CAM (3.3.7)
] ] [Gradient-based [NOTEARS⋄ (3.4.1)
GraN-DAG⋄ (3.4.2)
GAE (3.4.3)
DAG-GNN⋄ (3.4.4)
GOLEM⋄ (3.4.5)
DAG-NoCurl (3.4.6)
ENCO⋄ (3.4.7)] ] [Miscellaneous [FRITL (3.5.2)
HCM (3.5.3)
ETIO∗ (3.5.6)
JCI∗ (3.5.8)
Kg2Causal∗ (3.5.9)
Meta-RL (3.5.11)
LFCM (3.5.13)] ] ]

Figure 9: Taxonomy of some causal discovery approaches for I.I.D. data. The
approaches are classified based on their core contribution or the primary
strategy they adopt for causal structure recovery. The approaches that
leverage prior knowledge are marked by an $\ast$ symbol. Some of the gradient-
based optimization approaches that use a score function are indicated by a
$\diamond$ symbol. They are primarily classified as gradient-based methods
because of the use of gradient descent for optimization. However, they can be
a score-based method too as they compute data likelihood scores on the way.

##  3 Causal Discovery Algorithms for I.I.D. Data

Causal graphs are essential as they represent the underlying causal story
embedded in the data. There are two very common approaches to recovering the
causal structure from observational data, _i) Constraint-based_ (Spirtes et
al. (2000b), Spirtes (2001), Colombo et al. (2012)) and _ii) Score-based_
(Chickering (2002)). Among the other types of approaches, _functional causal
models (FCMs)-based_ (Shimizu et al. (2006), Hoyer et al. (2008)) approaches
and _hybrid_ approaches (Tsamardinos et al. (2006)) are noteworthy. Recently,
some _gradient-based_ approaches have been proposed based on neural networks
(Abiodun et al. (2018)) and a modified definition (Equation 2) of the
acyclicity constraint (Zheng et al. (2018), Yu et al. (2019)). Other
approaches include the ones that prioritize the use of _background knowledge_
and provides ways to incorporate prior knowledge and experts’ opinion into the
search process (Wang et al. (2020); Sinha & Ramsey (2021)). In this section,
we provide an overview of the causal discovery algorithms for I.I.D. data
based on the different types of approaches mentioned above. The algorithms
primarily distinguish from each other based on the core approach they follow
to perform causal discovery. We further discuss noteworthy similar approaches
specialized for non-I.I.D. or time series data in section 4.

###  3.1 Constraint-based

Testing for conditional independence (CI) is a core objective of constraint-
based causal discovery approaches. Conditional independence tests can be used
to recover the causal skeleton if the probability distribution of the observed
data is faithful to the underlying causal graph (Marx & Vreeken (2019)). Thus,
constraint-based approaches conduct CI tests between the variables to check
for the presence or absence of edges. These approaches infer the conditional
independencies within the data using the _d-separation criterion_ to search
for a DAG that entails these independencies, and detect which variables are
d-separated and which are d-connected (Triantafillou & Tsamardinos (2016)).
$X$ is conditionally independent of $Z$ given $Y$ i.e. $X$
$\perp\\!\\!\\!\perp$ $Z$ $|$ $Y$ in Figure 10 (a) and in Figure 10 (b), $X$
and $Z$ are independent, but are not conditionally independent given $Y$.
Table 4 lists different types of CI tests used by constraint-based causal
discovery approaches.

Figure 10: (a) $X$ $\perp\\!\\!\\!\perp$ $Z$ $|$ $Y$ and (b) $X$ and $Z$ are not conditionally independent given $Y$. Table 4: Types of conditional independence (CI) tests. Please refer to the study Runge (2018) for a detailed discussion on CI tests. |  Conditional Independence Test | Ref.
---|---|---
1. |  Conditional Distance Correlation (CDC) test | Wang et al. (2015)
2. |  Momentary Conditional Independence (MCI) | Runge et al. (2019)
3. |  Kernel-based CI test (KCIT) | Zhang et al. (2012)
4. |  Randomized Conditional Correlation Test (RCoT) | Strobl et al. (2019)
5. |  Generative Conditional Independence Test (GCIT) | Bellot & van der Schaar (2019)
6. |  Model-Powered CI test | Sen et al. (2017)
7. |  Randomized Conditional Independence Test (RCIT) | Strobl et al. (2019)
8. |  Kernel Conditional Independence Permutation Test | Doran et al. (2014)
9. |  Gaussian Processes and Distance Correlation-based (GPDC) | Rasmussen et al. (2006)
10. |  Conditional mutual information estimated with a k-nearest neighbor estimator (CMIKnn) | Runge (2018)

####  3.1.1 PC

The Peter-Clark (PC) algorithm (Spirtes et al. (2000b)) is one of the oldest
constraint-based algorithms for causal discovery. To learn the underlying
causal structure, this approach depends largely on conditional independence
(CI) tests. This is because it is based on the concept that two statistically
independent variables are not causally linked. The outcome of a PC algorithm
is a CPDAG. It learns the CPDAG of the underlying DAG in three steps: _Step 1
- Skeleton identification, Step 2 - V-structures determination, and Step 3 -
Edge orientations_. It starts with a fully connected undirected graph using
every variable in the dataset, then eliminates the unconditionally and
conditionally independent edges (skeleton detection), then it finds and
orients the v-structures or colliders (i.e. X → Y ← Z) based on the
d-separation set of node pairs, and finally orients the remaining edges based
on two aspects: i) availability of no new v-structures, and ii) not allowing
any cycle formation. The assumptions made by the PC algorithm include
acyclicity, causal faithfulness, and causal sufficiency. It is computationally
more feasible for sparse graphs. An implementation of this algorithm can be
found in the CDT repository
(https://github.com/ElementAI/causal_discovery_toolbox) and also, in the
gCastle toolbox (Zhang et al. (2021a)). A number of the constraint-based
approaches namely FCI, RFCI, PCMCI, PC-stable, etc. use the PC algorithm as a
backbone to perform the CI tests.

Figure 11:  Step-by-step workflow of the PC (Spirtes et al. (2000b))
algorithm.

####  3.1.2 FCI

The Fast Causal Inference (FCI) algorithm (Spirtes et al. (2000a)) is a
variant of the PC algorithm that can infer conditional independencies and
learn causal relations in the presence of many arbitrary latent and selection
variables. As a result, it is accurate in the large sample limit with a high
probability even when there exists hidden variables, and selection bias (Berk
(1983)). The first step of the FCI algorithm is similar to the PC algorithm
where it starts with a complete undirected graph to perform the skeleton
determination. After that, it requires additional tests to learn the correct
skeleton and has additional orientation rules. In the worst case, the number
of conditional independence tests performed by the algorithm grows
exponentially with the number of variables in the dataset. This can affect
both the speed and the accuracy of the algorithm in the case of small data
samples. To improve the algorithm, particularly in terms of speed, there exist
different variants such as the RFCI (Colombo et al. (2012)) and the Anytime
FCI (Spirtes (2001)) algorithms.

####  3.1.3 Anytime FCI

Anytime FCI (Spirtes (2001)) is a modified and faster version of the FCI
(Spirtes et al. (2000a)) algorithm. The number of CI tests required by FCI
makes it infeasible if the model has a large number of variables. Moreover,
when the FCI requires independence tests conditional on a large set of
variables, the accuracy decreases for a small sample size. The outer loop of
the FCI algorithm performs independence tests conditional on the increasing
size of variables. In the anytime FCI algorithm, the authors showed that this
outer loop can be stopped anytime during the execution for any smaller
variable size. As the number of variables in the conditional set reduces,
anytime FCI becomes much faster for the large sample size. More importantly,
it is also more reliable on limited samples since the statistical tests with
the lowest power are discarded. To support the claim, the authors provided
proof for the change in FCI that guarantees good results despite the
interruption. The result of the interrupted anytime FCI algorithm is still
valid, but as it cannot provide answers to most questions, the results could
be less informative compared to the situation if it was allowed to run
uninterrupted.

####  3.1.4 RFCI

Really Fast Causal Inference (RFCI) (Colombo et al. (2012)) is a much faster
variant of the traditional FCI for learning PAGs that uses fewer CI tests than
FCI. Unlike FCI, RFCI assumes that causal sufficiency holds. To ensure
soundness, RFCI performs some additional tests before orienting v-structures
and discriminating paths. It conditions only on subsets of the adjacency sets
and unlike FCI, avoids the CI tests given subsets of possible d-separation
sets which can become very large even for sparse graphs. As a result, the
number of these additional tests and the size of their conditioning sets are
small for sparse graphs which makes RFCI much faster and computationally
feasible than FCI for high-dimensional sparse graphs. Also, the lower
computational complexity of RFCI leads to high-dimensional consistency results
under weaker conditions than FCI.

####  3.1.5 FCI with Tiered Background Knowledge

Andrews et al. (2020) show that the Fast Causal Inference (FCI) algorithm
(Spirtes et al. (2000a)) is sound and complete with tiered background
knowledge (TBK). By _tiered background knowledge_ , it means any knowledge
where the variables may be partitioned into two or more mutually exclusive and
exhaustive subsets among which there is a known causal order. Tiered
background knowledge may arise in many different situations, including but not
limited to instrumental variables, data from multiple contexts and
interventions, and temporal data with contemporaneous confounding. The proof
that FCI is complete with TBK suggests that the algorithm is able to find all
of the causal relationships that are identifiable from tiered background
knowledge and observational data under the typical assumptions.

####  3.1.6 PC-stable

The independence tests in the original PC method are prone to errors in the
presence of a few samples. Additionally, because the graph is updated
dynamically, maintaining or deleting an edge incorrectly will affect the
neighboring sets of other nodes. As a result, the sequence in which the CI
tests are run will affect the output graph. Despite the fact that this order
dependency is not a significant issue in low-dimensional situations, it is a
severe problem in high-dimensional settings. To solve this problem, Colombo et
al. (2014) suggested changing the original PC technique to produce a stable
output skeleton that is independent of the input dataset’s variable ordering.
This approach, known as the stable-PC algorithm, queries and maintains the
neighbor (adjacent) sets of every node at each distinct level. Since the
conditioning sets of the other nodes are unaffected by an edge deletion at one
level, the outcome is independent of the variable ordering. They demonstrated
that this updated version greatly outperforms the original algorithm in high-
dimensional settings while maintaining the original algorithms’ low-
dimensional settings performance. However, this modification lengthens the
algorithm’s runtime even more by requiring additional CI checks to be done at
each level. The R-package pcalg contains the source code for PC-stable.

####  3.1.7 PKCL

Wang et al. (2020) proposed an algorithm, Prior-Knowledge-driven Local Causal
Structure Learning (PKCL), to discover the underlying causal mechanism between
bone mineral density (BMD) and its factors from clinical data. It first
discovers the neighbors of the target variables and then detects the
MaskingPCs to eliminate their effect. After that, it finds the spouse of
target variables utilizing the neighbors set. This way the skeleton of the
causal network is constructed. In the global stage, PKCL leverages the _Markov
blanket (MB)_ sets learned in the local stage to learn the global causal
structure in which prior knowledge is incorporated to guide the global
learning phase. Specifically, it learns the causal direction between feature
variables and target variables by combining the constraint-based and score-
based structure search methods. Also, in the learning phase, it automatically
adds casual direction according to the available prior knowledge.

###  3.2 Score-based

Score-based causal discovery algorithms search over the space of all possible
DAGs to find the graph that best explains the data. Typically, any score-based
approach has two main components: (i) a search strategy to explore the
possible search states or space of candidate graphs $G^{{}^{\prime}}$, and
(ii) a score function to assess the candidate causal graphs. The search
strategy along with a score function helps to optimize the search over the
space of all possible DAGs. More specifically, a score function
$S(G^{{}^{\prime}},D)$ maps causal graphs $G^{{}^{\prime}}$ to a numerical
score, based on how well $G^{{}^{\prime}}$ fits a given dataset $D$. A
commonly used score function to select causal models is the Bayesian
Information Criterion (BIC) (Schwarz (1978a)) which is defined below:

$\mathcal{S}(G^{{}^{\prime}},D)=-2\text{log}\
\mathcal{L}\\{G^{{}^{\prime}},D\\}+k\text{log}\ n,$ (3)

where $n$ is the number of samples in $D$, $k$ is the dimension of
$G^{{}^{\prime}}$ and $\mathcal{L}$ is the maximum-likelihood function
associated with the candidate graph $G^{{}^{\prime}}$. The lower the BIC
score, the better the model. BDeu, BGe, MDL, etc. (listed in Table 5) are some
of the other commonly used score functions. These objective functions are
optimized through a heuristic search for model selection. After evaluating the
quality of the candidate causal graphs using the score function, the score-
based methods output one or more causal graphs that achieve the highest score
(Huang et al. (2018b)). We discuss some of the well-known approaches in this
category below.

Figure 12: General components of a score-based causal discovery approach. Table 5: Some commonly used score functions for causal discovery. Please refer to the study Huang et al. (2018a) for a detailed discussion of the score functions. Score Function/Criterion | Ref.
---|---
Minimum description length (MDL) | Schwarz (1978b)
Bayesian information criterion (BIC) | Schwarz (1978a)
Akaike information criterion (AIC) | Akaike (1998)
Bayesian Dirichlet equivalence score (BDeU) | Buntine (1991)
Bayesian metric for Gaussian networks (BGe) | Geiger & Heckerman (1994)
Factorized normalized maximum likelihood (fNML) | Silander et al. (2008)

####  3.2.1 GES

Greedy Equivalence Search (GES) (Chickering (2002)) is one of the oldest
score-based causal discovery algorithms that perform a greedy search over the
space of equivalence classes of DAGs. Each search state is represented by a
CPDAG where some insert and delete operators allow for single-edge additions
and deletions respectively. Primarily GES works in two phases: i) Forwards
Equivalence Search (FES), and ii) Backward Equivalence Search (BES). In the
first phase, FES starts with an empty CPDAG (no-edge model), and greedily adds
edges by taking into account every single-edge addition that could be
performed to every DAG in the current equivalence class. After an edge
modification is done to the current CPDAG, a score function is used to score
the model. If the new score is better than the current score, only then the
modification is allowed. When the forward phase reaches a local maximum, the
second phase, BES starts where at each step, it takes into account all single-
edge deletions that might be allowed for all DAGs in the current equivalence
class. The algorithm terminates once the local maximum is found in the second
phase. Implementation of GES is available at the following Python packages:
Causal Discovery Toolbox or CDT (Kalainathan & Goudet (2019)) and gCastle
(Zhang et al. (2021a)). GES assumes that the score function is decomposable
and can be expressed as a sum of the scores of individual nodes and their
parents. A summary workflow of GES is shown in Figure 13.

Figure 13: Different stages in the GES algorithm.

####  3.2.2 FGS

Fast Greedy Search (FGS) (Ramsey (2015)) is another score-based method that is
an optimized version of the GES algorithm (Chickering (2002)). This optimized
algorithm is based on the faithfulness assumption and uses an alternative
method to reduce scoring redundancy. An ascending list $L$ is introduced which
stores the score difference of arrows. After making a thorough search, the
first edge e.g. $X$ $\rightarrow$ $Y$ is inserted into the graph and the graph
pattern is reverted. For variables that are adjacent to $X$ or $Y$ with
positive score differences, new edges are added to $L$. This process in the
forward phase repeats until the $L$ becomes empty. Then the reverse phase
starts, filling the list $L$ and continuing until $L$ is empty. This study
considered the experiment where GES was able to search over 1000 samples with
50,000 variables in 13 minutes using a 4-core processor and 16GB RAM computer.
Following the new scoring method, FGS was able to complete the task with 1000
samples on 1,000,000 variables for sparse models in 18 hours using a
supercomputer having 40 processors and 384GB RAM at the Pittsburgh
Supercomputing Center. The code for FGS is available on GitHub as a part of
the Tetrad project: https://github.com/cmu-phil/tetrad.

####  3.2.3 SGES

Selective Greedy Equivalence Search (SGES) (Chickering & Meek (2015)) is
another score-based causal discovery algorithm that is a restrictive variant
of the GES algorithm (Chickering (2002)). By assuming perfect generative
distribution, SGES provides a polynomial performance guarantee yet maintains
the asymptotic accuracy of GES. While doing this, it is possible to keep the
algorithm’s large sample guarantees by ignoring all but a small fraction of
the backward search operators that GES considered. In the forward phase, SGES
uses a polynomial number of insert operation calls to the score function. In
the backward phase, it consists of only a subset of delete operators of GES
which include, consistent operators to preserve GES’s consistency over large
samples. The authors demonstrated that, for a given set of graph-theoretic
complexity features, such as maximum-clique size, the maximum number of
parents, and v-width, the number of score assessments by SGES can be
polynomial in the number of nodes and exponential in these complexity
measurements.

####  3.2.4 RL-BIC

RL-BIC is a score-based approach that uses _Reinforcement Learning (RL)_ and a
BIC score to search for the DAG with the best reward (Zhu et al. (2019)). For
data-to-graph conversion, it uses an _encoder-decoder architecture_ that takes
observational data as input and generates graph adjacency matrices that are
used to compute rewards. The reward incorporates a BIC score function and two
penalty terms for enforcing acyclicity. The _actor-critic RL algorithm_ is
used as a _search strategy_ and the final output is the causal graph that
achieves the best reward among all the generated graphs. The approach is
applicable to small and medium graphs of up to 30 nodes. However, dealing with
large and very large graphs is still a challenge for it. This study mentions
that their future work involves developing a more efficient and effective
score function since computing scores is much more time-consuming than
training NNs. The original implementation of the approach is available at:
https://github.com/huawei-noah/trustworthyAI.

Figure 14: Components of the RL-BIC (Zhu et al. (2019)) approach.

####  3.2.5 A* search

Xiang & Kim (2013) proposed a one-stage method for learning sparse network
structures with continuous variables using the A* search algorithm with lasso
in its scoring system. This method increased the computational effectiveness
of popular exact methods based on dynamic programming. The study demonstrated
how the proposed approach achieved comparable or better accuracy with
significantly faster computation time when compared to two-stage approaches,
including L1MB and SBN. Along with that, a heuristic approach was added that
increased A* lasso’s effectiveness while maintaining the accuracy of the
outcomes. In high-dimensional spaces, this is a promising approach for
learning sparse Bayesian networks.

####  3.2.6 Triplet A*

Lu et al. (2021) uses the _A* exhaustive search_ (Yuan & Malone (2013))
combined with an optimal BIC score that requires milder assumptions on data
than conventional CD approaches to guarantee its asymptotic correctness. The
optimal BIC score combined with the exhaustive search finds the MEC of the
true DAG if and only if the true DAG satisfies the optimal BIC Condition. To
gain scalability, they also developed an approximation algorithm for complex
large systems based on the A* method. This extended approach is named Triplet
A* which can scale up to more than 60 variables. This extended method is
rather general and can be used to scale up other exhaustive search approaches
as well. Triplet A* can particularly handle linear Gaussian and non-Gaussian
networks. It works in the following way. Initially, it makes a guess about the
parents and children of each variable. Then for each variable $X$ and its
neighbors $(Y,Z)$, it forms a cluster consisting of $X,Y,Z$ with their direct
neighbors and runs an exhaustive search on each cluster. Lastly, it combines
the results from all clusters. The study shows that empirically Triplet A*
outperforms GES for large dense networks.

####  3.2.7 KCRL

Prior Knowledge-based Causal Discovery Framework with Reinforcement Learning
a.k.a. KCRL (Hasan & Gani (2022)) is a framework for causal discovery that
utilizes prior knowledge as constraints and penalizes the search process for
violation of these constraints. This utilization of background knowledge
significantly improves performance by reducing the search space, and also,
enabling a faster convergence to the optimal causal structure. KCRL leverages
reinforcement learning (RL) as the search strategy where the RL agent is
penalized each time for the violation of any imposed knowledge constraints. In
the KCRL framework (Figure 15), at first, the observational data is fed to an
RL agent. Here, data-to-adjacency matrix conversion is done using an encoder-
decoder architecture which is a part of the RL agent. At every iteration, the
agent produces an equivalent adjacency matrix of the causal graph. A
comparator compares the generated adjacency matrix with the true causal edges
in the prior knowledge matrix $P_{m}$, and thereby, computes a penalty $p$ for
the violation of any ground truth edges in the produced graph. Each generated
graph is also scored using a standard scoring function such as BIC. A reward
$R$ is estimated as a sum of the BIC score $S_{BIC}$, the penalty for
acyclicity $h(W)$, and $\beta$ weighted prior knowledge penalty $\beta p$.
Finally, the entire process halts when the stopping criterion $S_{c}$ is
reached, and the best-rewarded graph is the final output causal graph.
Although originally KCRL was designed for the healthcare domain, it can be
used in any other domain for causal discovery where some prior knowledge is
available. Code for KCRL is available at https://github.com/UzmaHasan/KCRL.

$R=S_{BIC}+\beta p+h(W)$ (4)

Figure 15: The KCRL (Hasan & Gani (2022)) framework.

Another recent method called KGS (Hasan & Gani (2023)) leverages prior causal
information such as the presence or absence of a causal edge to guide a greedy
score-based causal discovery process towards a more restricted and accurate
search space. It demonstrates how the search space as well as scoring
candidate graphs can be reduced when different edge constraints are leveraged
during a search over equivalence classes of causal networks. It concludes that
any type of edge information is useful to improve the accuracy of the graph
discovery as well as the run time.

####  3.2.8 ILP-based structure learning

Bartlett & Cussens (2017) looked into the application of integer linear
programming (ILP) to the structure learning problem. To boost the
effectiveness of ILP-based Bayesian network learning, they suggested adding
auxiliary implied constraints. Experiments were conducted to determine the
effect of each constraint on the optimization process. It was discovered that
the most effective configuration of these constraints could significantly
boost the effectiveness and speed of ILP-based Bayesian network learning. The
study made a significant contribution to the field of structure learning and
showed how well ILP can perform under non-essential constraints.

###  3.3 Functional Causal Model-based

Functional Causal Model (FCM) based approaches describe the causal
relationship between variables in a specific functional form. FCMs represent
variables as a function of their parents (direct causes) together with an
independent noise term $E$ (see Equation 5) (Zhang et al. (2015)). FCM-based
methods can distinguish among different DAGs in the same equivalence class by
imposing additional assumptions on the data distributions and/or function
classes (Zhang et al. (2021b)). Some of the noteworthy FCM-based causal
discovery approaches are listed below.

$X=f(PA_{X})+E$ (5)

Figure 16: A functional causal model (FCM) with four variables.

####  3.3.1 L$i$NGAM

Linear Non-Gaussian Acyclic Model (LiNGAM) aims to discover the causal
structure from observational data under the assumptions that the data
generating process is linear, there are no unobserved confounders, and noises
have non-Gaussian distributions with non-zero variances (Shimizu et al.
(2006)). It uses the statistical method known as independent component
analysis (ICA) (Comon (1994)), and states that when the assumption of non-
Gaussianity is valid, the complete causal structure can be estimated. That is,
the causal direction is identifiable if the variables have a linear relation,
and the noise ($\varepsilon$) distribution is non-Gaussian in nature. Figure
17 depicts three scenarios where when $X$ and $\varepsilon$ are Gaussian (case
1), the predictor and regression residuals are independent of each other. For
the other two cases, $X$ and $\varepsilon$ are non-Gaussian, and we see that
for the regression in the anti-causal or backward direction ($X$ given $Y$),
the regression residual and the predictor are not independent as earlier. That
is, for the non-Gaussian cases, independence between regression residual and
predictor occurs only for the correct causal direction. There are 3 properties
of a LiNGAM. First, the variables $x_{i}={x_{1},x_{2},...,x_{n}}$ are arranged
in a causal order $k(i)$ such that the cause always preceedes the effect.
Second, each variable $x_{i}$ is assigned a value as per the Equation 6 where
$e_{i}$ is the noise/disturbance term and $b_{ij}$ denotes the causal strength
between $x_{i}$ and $x_{j}$. Third, the exogenous noise $e_{i}$ follows a non-
Gaussian distribution, with zero mean and non-zero variance, and are
independent of each other which implies that there is no hidden confounder.
Python implementation of the LiNGAM algorithm is available at
https://github.com/cdt15/lingam as well as in the gCastle package (Zhang et
al. (2021b)). Any standard ICA algorithm which can estimate independent
components of many different distributions can be used in LiNGAM. However, the
original implementation uses the FastICA (Hyvarinen (1999)) algorithm.

$x_{i}=\sum_{k(j)<k(i)}b_{ij}x_{j}+e_{i}$ (6)

Figure 17: Causal asymmetry between two variables having a linear relation
(Glymour et al. (2019)). Here, the causal direction is from $X$ to $Y$. A
total of three scenarios are depicted where both $X$ and $\varepsilon$ follow
the i) Gaussian, ii) Uniform, or iii) Super-Gaussian distribution for each of
the scenarios.

####  3.3.2 ANM

Hoyer et al. (2008) performs causal discovery with additive noise models
(ANMs) and provides a generalization of the linear non-Gaussian causal
discovery framework to deal with nonlinear functional dependencies where the
variables have an additive noise. It mentions that nonlinear causal
relationships typically help to break the symmetry between the observed
variables and help in the identification of causal directions. ANM assumes
that the data generating process of the observed variables is as per the
Equation 7 where a variable $x_{i}$ is a function of its parents and the noise
term $e_{i}$ which is an independent additive noise. An implementation of ANM
is available in the gCastle package (Zhang et al. (2021a)).

$x_{i}=f(PA_{x_{i}})+e_{i}$ (7)

####  3.3.3 PNL

Post-nonlinear (PNL) acyclic causal model with additive noise (Zhang &
Hyvärinen (2010)) is a highly realistic model where each observed continuous
variable is made up of additive noise-filled nonlinear functions of its
parents, followed by a nonlinear distortion. The influence of sensor
distortions, which are frequently seen in practice, is taken into account by
the second stage’s nonlinearity. A two-step strategy is proposed to separate
the cause from the effect in a two-variable situation, consisting of
restricted nonlinear ICA followed by statistical independence tests. The PNL
model was able to effectively separate causes from effects when applied to
solve the "CauseEffectPairs" task proposed by Mooij & Janzing (2010) in the
Pot-luck challenge. That is, it successfully distinguished the cause from the
effect, even if the nonlinear function of the cause is not invertible.

####  3.3.4 Direct-L$i$NGAM

Shimizu et al. (2011) proposed DirectLiNGAM, a direct method for learning a
linear non-Gaussian structural equation model (SEM) which is a direct method
to estimate causal ordering and connection strengths based on non-Gaussianity.
This approach estimates a causal order of variables by successively reducing
each independent component from given data in the model which is completed in
steps equal to the number of the variables in the model. Once the causal order
of variables is identified, their connection strengths are estimated using
conventional covariance-based methods such as least squares and maximum
likelihood approaches. If the data strictly follows the model i.e. if all the
model assumptions are met and the sample size is infinite, it converges to the
right solution within a small number of steps. If some prior knowledge on a
part of the structure is available, it suggests using those for more efficient
learning. Doing so will reduce the number of causal orders and connection
strengths to be estimated. Its implementation can be found at:
https://github.com/huawei-noah/trustworthyAI/tree/master/gcastle.

####  3.3.5 SAM

Kalainathan et al. (2018) proposed the algorithm known as Structural Agnostic
Modeling (SAM) that uses an _adversarial learning_ approach to find the causal
graphs. Particularly, it searches for an FCM using _Generative Adversarial
Neural-networks (GANs)_ and enforces the discovery of sparse causal graphs
through adequate regularization terms. A learning criterion that combines
distribution estimation, sparsity, and acyclicity constraints is used to
enforce the end-to-end optimization of the graph structure and parameters
through stochastic gradient descent. SAM leverages both conditional
independencies and distributional asymmetries in the data to find the
underlying causal mechanism. It aims to achieve an optimal complexity/fit
trade-off while modeling the causal mechanisms. SAM enforces the acyclicity
constraint of a DAG using the function in Equation 8 where, $A$ is the
adjacency matrix of the ground-truth graph $G$, and $d$ denotes the total
number of nodes in $G$. The latest implementation of SAM is available in the
CDT package (Kalainathan & Goudet (2019)). Also, an older version of SAM is
available at https://github.com/Diviyan-Kalainathan/SAM.

$\sum_{i=1}^{d}=\frac{\operatorname{tr}(A^{i})}{i!}=0$ (8)

####  3.3.6 CGNN

Causal Generative Neural Networks (CGNN) is an FCM-based framework that uses
_neural networks (NNs)_ to learn the joint distribution of the observed
variables (Goudet et al. (2018)). Particularly, it uses a generative model
that minimizes the _maximum mean discrepancy_ (MMD) between the generated and
observed data. CGNN has a high computational cost. However, it proposes an
approximate learning criterion to scale the computational cost to linear
complexity in the number of observations. This framework can also be used to
simulate interventions on multiple variables in the dataset. An implementation
of CGNN in Pytorch is available at
https://github.com/FenTechSolutions/CausalDiscoveryToolbox.

####  3.3.7 CAM

Causal Additive Model (CAM) is a method for estimating high-dimensional
additive structural equation models which are logical extensions of linear
structural equation models (Bühlmann et al. (2014)). In order to address the
difficulties of computation and statistical accuracy in the absence of prior
knowledge about underlying structure, the authors established consistency of
the maximum likelihood estimator and developed an effective computational
algorithm. The technique was demonstrated using both simulated and actual data
and made use of tools in sparse regression techniques. The authors also
discussed identifiability problems and the enormous size of the space of
potential models, which presents significant computational and statistical
accuracy challenges.

####  3.3.8 CAREFL

Causal Autoregressive Flows (CAREFL) uses _autoregressive flow models_ (Huang
et al. (2018c)) for causal discovery by interpreting the ordering of variables
in an autoregressive flow based on structural equation models (SEMs)
(Khemakhem et al. (2021)). In general, SEMs define a generative model for data
based on causal relationships. CAREFL shows that particularly affine flows
define a new class of causal models where the noise is modulated by the cause.
For such models, it proves a new causal identifiability result that
generalizes additive noise models. To learn the causal structure efficiently,
it selects the ordering with the highest test log-likelihood and reports a
measure of causal direction based on the likelihood ratio for non-linear SEMs.
Autoregressive flow models also enable CAREFL to evaluate interventional
queries by fixing the interventional variable while sampling from the flow.
Moreover, the invertible property of autoregressive flows facilitates
counterfactual queries as well. Code implementation of CAREFL is available at
https://github.com/piomonti/carefl.

###  3.4 Gradient-based

Some of the recent studies in causal discovery formulate the structure
learning problem as a continuous optimization task using the least squares
objective and an algebraic characterization of DAGs (Zheng et al. (2018), Ng
et al. (2020)). Specifically, the combinatorial structure learning problem has
been transformed into a continuous one and solved using gradient-based
optimization methods (Ng et al. (2019)). These methods leverage gradients of
an objective function with respect to a parametrization of a DAG matrix. Apart
from the usage of well-studied gradient-based solvers, they also leverage GPU
acceleration which has changed the nature of the task (Ng et al. (2020)).
Furthermore, to accelerate the task they often employ deep learning models
that are capable of capturing complex nonlinear mappings (Yu et al. (2019)).
As a result, they usually have a faster training time as deep learning is
known to be highly parallelizable on GPU, which gives a promising direction
for causal discovery with gradient-based methods (Ng et al. (2019)). In
general, these methods are more global than other approximate greedy methods.
This is because they update all edges at each step based on the gradient of
the score and as well as based on the acyclicity constraint.

####  3.4.1 NOTEARS

DAGs with NO TEARS (Zheng et al. (2018)) is a recent breakthrough in the field
of causal discovery that formulates the structure learning problem as a purely
continuous constrained optimization task. It leverages an algebraic
characterization of DAGs and provides a novel characterization of acyclicity
that allows for a smooth global search, in contrast to a combinatorial local
search. The full form of the acronym NOTEARS is Non-combinatorial Optimization
via Trace Exponential and Augmented lagRangian for Structure learning which
particularly handles linear DAGs. It assumes a linear dependence between
random variables and thus models data $D$ as a structural equation model. To
discover the causal structure, it imposes the proposed acyclicity function
(Equation 10) as a constraint combined with a weighted adjacency matrix $W$
with least squares loss. The algorithm aims to convert the traditional
combinatorial optimization problem into a continuous constrained optimization
task by leveraging an algebraic characterization of DAGs via the trace
exponential acyclicity function as follows:

$\underset{\text{subject to }G(W)\in
DAGs}{\mathop{\min}_{W\in\mathbb{R}^{d\times
d}}F(W)}\iff\underset{\text{subject to
}h(W)=0}{\mathop{\min}_{W\in\mathbb{R}^{d\times d}}F(W)},$ (9)

where $G(W)$ is a graph with $d$ nodes induced by the weighted adjacency
matrix $W$, $F:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}$ is a regularized
score function with a least-square loss $\ell$, and $h:\mathbb{R}^{d\times
d}\rightarrow\mathbb{R}$ is a smooth function over real matrices that enforces
acyclicity. Overall, the approach is simple and can be executed in about 50
lines of Python code. Its implementation in Python is publicly available at
https://github.com/xunzheng/notears. The acyclicity function proposed in
NOTEARS is as follows where $\circ$ is the Hadamard product and $e^{W\circ W}$
is the matrix exponential of $W\circ W$.

$h(W)=\operatorname{tr}(e^{W\circ W})-d=0$ (10)

####  3.4.2 GraN-DAG

Gradient-based Neural DAG Learning (GraN-DAG) is a causal structure learning
approach that uses _neural networks (NNs)_ to deal with non-linear causal
relationships (Lachapelle et al. (2019)). It uses a stochastic gradient method
to train the NNs to improve scalability and allow implicit regularization. It
formulates a _novel characterization of acyclicity_ for NNs based on NOTEARS
(Zheng et al. (2018)). To ensure acyclicity in non-linear models, it uses an
argument similar to NOTEARS and applies it first at the level of neural
network paths and then at the graph paths level. For regularization, GraN-DAG
uses a procedure called preliminary neighbors selection (PNS) to select a set
of potential parents for each variable. It uses a final pruning step to remove
the false edges. The algorithm works well mostly in the case of non-linear
Gaussian additive noise models. An implementation of GraN-DAG can be found at
https://github.com/kurowasan/GraN-DAG.

####  3.4.3 GAE

Graph Autoencoder (GAE) approach is a gradient-based approach to causal
structure learning that uses a _graph autoencoder framework_ to handle
nonlinear structural equation models (Ng et al. (2019)). GAE is a special case
of the causal additive model (CAM) that provides an alternative generalization
of NOTEARS for handling nonlinear causal relationships. GAE is easily
applicable to vector-valued variables. The architecture of GAE consists of a
variable-wise encoder and decoder which are basically multi-layer perceptrons
(MLPs) with shared weights across all variables $X_{i}$. The encoder-decoder
framework allows the reconstruction of each variable $X_{i}$ to handle the
nonlinear relations. The final goal is to optimize the reconstruction error of
the GAE with $l_{1}$ penalty where the optimization problem is solved using
the augmented Lagrangian method (Nemirovsky (1999)). The approach is
competitive in terms of scalability as it has a near-linear training time when
scaling up the graph size to 100 nodes. Also, in terms of time efficiency, GAE
performs well with an average training time of fewer than 2 minutes even for
graphs of 100 nodes. Its implementation can be found at the gCastle (Zhang et
al. (2021a)) repository.

####  3.4.4 DAG-GNN

DAG Structure Learning with Graph Neural Networks (DAG-GNN) is a graph-based
deep generative model that tries to capture the sampling distribution faithful
to the ground-truth DAG (Yu et al. (2019)). It leverages variational inference
and a parameterized pair of _encoder-decoders_ with specially designed _graph
neural networks (GNN)_. Particularly, it uses _Variational Autoencoders
(VAEs)_ to capture complex data distributions and sample from them. The
weighted adjacency matrix $W$ of the ground-truth DAG is a learnable parameter
with other neural network parameters. The VAE model naturally handles various
data types both continuous and discrete in nature. In this study, the authors
also propose a _variant of the acyclicity function_ (Equation 11) which is
more suitable and practically convenient for implementation with the existing
deep learning methods. In the acyclicity function, $d$ = the number of nodes,
$\alpha$ is a hyperparameter, and $I$ is an identity matrix. An implementation
of the DAG-GNN algorithm is available at https://github.com/fishmoon1234/DAG-
GNN.

$\operatorname{tr}[(I+\alpha W\circ W)^{d}]-d=0$ (11)

####  3.4.5 GOLEM

Gradient-based Optimization of DAG-penalized Likelihood for learning linear
DAG Models (GOLEM) is a _likelihood-based_ causal structure learning approach
with _continuous unconstrained optimization_ (Ng et al. (2020)). It studies
the asymptotic role of the sparsity and DAG constraints for learning DAGs in
both linear Gaussian and non-Gaussian cases. It shows that when the
optimization problem is formulated using a likelihood-based objective instead
of least squares (used by NOTEARS), then instead of a hard DAG constraint,
applying only soft sparsity and DAG constraints is enough for learning the
true DAG under mild assumptions. Particularly, GOLEM tries to optimize the
score function in Equation 12 w.r.t. the weighted adjacency matrix $B$
representing a directed graph. Here, $L(B;x)$ is the maximum likelihood
estimator, $R_{sparse}(B)$ is a penalty to encourage sparsity (i.e. fewer
edges), and $R_{DAG}(B)$ is the penalty that enforces DAGness on $B$.

$S(B;x)=L(B;x)+R_{sparse}(B)+R_{DAG}(B)$ (12)

In terms of denser graphs, GOLEM seems to outperform NOTEARS since it can
reduce the number of optimization iterations which makes it robust in terms of
scalability. With gradient-based optimization and GPU acceleration, it can
easily handle thousands of nodes while retaining high accuracy. An
implementation of GOLEM can be found at the gCastle (Zhang et al. (2021a))
repository.

####  3.4.6 DAG-NoCurl

DAG-NoCurl also known as DAGs with No Curl uses a two-step procedure for the
causal DAG search (Yu et al. (2021)). At first, it finds an initial cyclic
solution to the optimization problem and then employs the _Hodge
decomposition_ (Bhatia et al. (2012)) of graphs to learn an acyclic graph by
projecting the cyclic graph to the gradient of a potential function. The goal
of this study is to investigate how the causal structure can be learned
without any explicit DAG constraints by directly optimizing the DAG space. To
do so, it proposes the method DAG-NoCurl based on the graph Hodge theory that
implicitly enforces the acyclicity of the learned graph. As per the Hodge
theory on graphs (Lim (2020)), a DAG is a sum of three components: a curl-
free, a divergence-free, and a harmonic component. The curl-free component is
an acyclic graph that motivates the naming of this approach. An implementation
of the method can be found at the link https://github.com/fishmoon1234/DAG-
NoCurl.

####  3.4.7 ENCO

Efficient Neural Causal Discovery without Acyclicity Constraints (ENCO) _uses
both observational and interventional data_ by modeling a probability for
every possible directed edge between pairs of variables (Lippe et al. (2021)).
It formulates the graph search as an optimization of independent edge
likelihoods, with the edge orientation being modeled as a separate parameter.
This approach guarantees convergence when interventions on all variables are
available and do not require explicitly constraining the score function with
respect to acyclicity. However, the algorithm works on partial intervention
sets as well. Experimental results suggest that ENCO is robust in terms of
scalability, and is able to detect latent confounders. When applied to large
networks having 1000 nodes, it is capable of recovering the underlying
structure due to the benefit of its low-variance gradient estimators. The
source code of ENCO is available at this site:
https://github.com/phlippe/ENCO.

Figure 18: Graph optimization mechanism of ENCO.

####  3.4.8 MCSL

Masked Gradient-based Causal Structure Learning (MCSL) (Ng et al. (2022))
utilizes a reformulated structural equation model (SEM) for causal discovery
using gradient-based optimization that leverages the _Gumbel-Softmax approach_
(Jang et al. (2016)). This approach is used to approximate a binary adjacency
matrix and is often used to approximate samples from a categorical
distribution. MCSL reformulates the SEM with additive noises in a form
parameterized by the binary graph adjacency matrix. It states that, if the
original SEM is identifiable, then the adjacency matrix can be identified up
to super-graphs of the true causal graph under some mild conditions. For
experimentation, MCSL uses multi-layer perceptrons (MLPs), particularly having
4-layers as the model function which is denoted as MCSL-MLP. An implementation
of the approach can be found in the gCastle (Zhang et al. (2021a)) package.

####  3.4.9 DAGs with No Fears

Wei et al. (2020) provides an in-depth analysis of the NOTEARS framework for
causal structure learning. The study proposed a local search post-processing
algorithm that significantly increased the precision of NOTEARS and other
algorithms and deduced Karush-Kuhn-Tucker (KKT) optimality conditions for an
equivalent reformulation of the NOTEARS problem. Additionally, the authors
compared the effectiveness of NOTEARS and Abs-KKTS on various graph types and
discovered that Abs-KKTS performed better than NOTEARS in terms of accuracy
and computational efficiency. The authors concluded that this work improved
the understanding of optimization-based causal structure learning and may
result in further advancements in precision and computational effectiveness.
The code implementation is available at https://github.com/skypea/DAG_No_Fear.

###  3.5 Miscellaneous Approaches

Apart from the types of approaches mentioned so far, there are some other
causal discovery approaches that use some specialized or unique techniques to
search for the graph that best describes the data. There also exists some
methods that are specialized to handle latent or unobserved confounders. Also,
there are some approaches that are hybrid in nature, i.e. they are based on
the combination of constraint-based, score-based, FCM-based, gradient-based,
etc. causal discovery approaches. For example, some approaches integrate
conditional independence testing along with score functions to design a hybrid
approach for causal discovery. A detailed discussion can be found below.

####  3.5.1 MMHC

Max-Min Hill Climbing (MMHC) is a hybrid causal discovery technique that
incorporates the concepts from both score-based and constraint-based
algorithms (Tsamardinos et al. (2006)). A challenge in causal discovery is the
identification of causal relationships within a reasonable time in the
presence of thousands of variables. MMHC can reliably learn the causal
structure in terms of time and quality for high-dimensional settings. MMHC is
a two-phase algorithm that assumes faithfulness. In the first phase, MMHC uses
Max-Min Parents and Children (MMPC) (Tsamardinos et al. (2003)) to initially
learn the skeleton of the network. In the second phase, using a greedy
Bayesian hill-climbing search, the skeleton is oriented. In the sample limit,
MMHC’s skeleton identification phase is reliable, but the orientation phase
offers no theoretical assurances. From the results of the experiments
performed, MMHC outperformed PC (Spirtes et al. (2000b)), Sparse Candidate
(Friedman et al. (2013)), Optimal Reinsertion (Moore & Wong (2003)), and GES
(Chickering (2002)) in terms of computational efficiency. Considering the
quality of reconstruction, MMHC performs better than all the above-mentioned
algorithms except for GES when the sample size is 1000. The authors also
proved the correctness of the results. The implementation of MMHC is available
at http://www.dsl-lab.org/supplements/mmhcpaper/mmhcindex.html as part of
Causal Explorer 1.3, a library of Bayesian network learning and local causal
discovery methods.

####  3.5.2 FRITL

To discover causal relationships in linear and non-Gaussian models, Chen et
al. (2021) proposed a hybrid model named FRITL. FRITL works in the presence or
absence of latent confounders by incorporating independent noise-based
techniques and constraint-based techniques. FRITL makes causal Markov
assumption, causal faithfulness assumption, linear acyclic non-Gaussianity
assumption, and one latent confounder assumption. In the first phase of FRITL,
the FCI algorithm is used to generate asymptotically accurate results.
Unfortunately, relatively few unconfounded direct causal relations are
normally determined by the FCI since it always reveals the presence of
confounding factors. In the second phase, FRITL identifies the unconfounded
causal edges between observable variables within just those neighboring
pairings that have been influenced by the FCI results. The third stage can
identify confounders and the relationships that cause them to affect other
variables by using the Triad condition (Cai et al. (2019)). If further causal
relationships remain, Independent Component Analysis (ICA) is finally applied
to a notably reduced group of graphs. The authors also theoretically proved
that the results obtained from FRITL are efficient and accurate. FRITL
produces results that are in close accord with neuropsychological opinion and
in exact agreement with a causal link that is known from the experimental
design when applied to real functional magnetic MRI data and the SACHS (Sachs
et al. (2005)) dataset.

Figure 19: Stages of the FRITL model.

####  3.5.3 HCM

Most of the causal discovery algorithms are applicable only to either discrete
or continuous data. However, in reality, we often have to work with mixed-type
data (e.g., shopping behavior of people) which don’t receive enough attention
in causal discovery. Li et al. (2022) proposed the approach _H ybrid Causal
Discovery on Mixed-type Data (HCM)_ to identify causal relationships with
mixed variables. HCM works under the causal faithfulness and causal Markov
assumption. HCM has three phases where in the first phase, the skeleton graph
is learned in order to limit the search space. To do this, they used the PC-
stable approach along with their proposed Mixed-type Randomized Causal
Independence Test (MRCIT) which can handle mixed-type data. They also
introduced a generalized score function called Cross-Validation based Mixed
Information Criterion (CVMIC). In the second phase, starting with an empty
DAG, they add edges to the DAG based on the highest CVMIC score. In order to
reduce false positives, the learned causal structure is pruned using MRCIT
once again in the final phase with a slightly bigger conditional set. They
compared their approach with other causal discovery approaches for mixed data
and showed HCM’s superiority. However, they didn’t consider any unobserved
confounders in the dataset which allows for further improvement. They made the
code available on the following GitHub site: https://github.com/DAMO-DI-
ML/AAAI2022-HCM.

Figure 20: Different phases of the method HCM.

####  3.5.4 SADA

One of the biggest limitations of the traditional causal discovery methods is
that these models cannot identify causal relations when the problem domain is
large or there is a small number of samples available. To solve this problem,
Cai et al. (2013) proposed a _Split-and-Merge_ causal discovery method named
SADA which assumes causal faithfulness. Even in situations when the sample
size is substantially less than the total number of variables, SADA can
reliably identify the causal factors. SADA divides the main problem into two
subproblems and works in three phases. Initially, SADA separates the variables
of the causal model into two sets $V_{1}$ and $V_{2}$ using a causal cut set
$C$ where all paths between $V_{1}$ and $V_{2}$ are blocked by $C$. This
partitioning is continued until the variables in each subproblem are less than
some threshold. In the next phase, any arbitrary causal algorithm is applied
to both subproblems and the causal graphs are generated. Here, they used
LiNGAM as the causal algorithm. Then these graphs are merged in the final
step. But to handle the conflicts while merging, they only kept the most
significant edge and eliminated the others whenever there existed multiple
causal paths between two variables in the opposite direction. They compared
the performance of SADA against baseline LiNGAM (without splitting and
merging), and the results showed that SADA achieved better performance in
terms of the metrics precision, recall, and F1 score.

####  3.5.5 CORL

Ordering-based Causal Discovery with Reinforcement Learning (CORL) formulates
the ordering search problem as a _multi-step Markov decision process_ (MDP) to
learn the causal graph (Wang et al. (2021)). It implements the ordering
generating process with an _encoder-decoder architecture_ and finally uses RL
to optimize the proposed model based on the reward mechanisms designed for
each order. A generated ordering is then processed using variable selection to
obtain the final causal graph. According to the empirical results, CORL
performs better than existing RL-based causal discovery approaches. This could
happen because CORL does not require computing the matrix exponential term
with O($d^{3}$) cost because of using ordering search. CORL is also good in
terms of scalability and has been applied to graphs with up to 100 nodes. The
gCastle package contains an implementation of CORL.

####  3.5.6 ETIO

ETIO is a versatile _logic-based_ causal discovery algorithm specialized for
business applications (Borboudakis & Tsamardinos (2016)). Its features include
i) the ability to utilize prior causal knowledge, ii) addressing selection
bias, hidden confounders, and missing values in data, and iii) analyzing data
from pre and post-interventional distribution. ETIO follows a _query-based
approach_ , where the user queries the algorithm about the causal relations of
interest. In the first step, ETIO performs several CI tests on the input
dataset. Particularly, it performs non-Bayesian tests that return p-values of
the null hypothesis of conditional independencies. Then it employs an
empirical Bayesian method that converts the p-values of dependencies and
interdependencies into probabilities. Later, it selects a consistent subset of
dependence, and prior knowledge constraints to resolve conflicts which are
ranked in order of confidence. Particularly, ETIO imposes an m-separation
constraint if a given independence is more probable than the corresponding
dependence. These imposed constraints are the ones that correspond to test
results, in order of probability, while removing conflicting test results.
Finally, it identifies all invariant features based on input queries using the
well-known declarative programming language, answer set programming (Gelfond &
Lifschitz (1988)).

####  3.5.7 $b$QCD

Discovering causal relationships from observational data has been a
challenging task, especially for the bivariate cases as it is difficult to
determine whether there actually exists a cause-effect relationship or whether
it is the effect of a hidden confounder. Tagasovska et al. (2020) proposed the
approach bivariate Quantile Causal Discovery (bQCD) to determine causal
relationships in bivariate settings. Although they made no assumptions on the
class of causal mechanisms, they did assume that there exists no confounder,
feedback, or selection bias. They utilized _quantile scoring_ in place of
Kolmogorov complexity (Kolmogorov (1963)), and used conditional quantiles,
pinball loss instead of conditional mean, and squared loss. The approach bQCD
performs almost similarly to the state-of-the-art techniques but it is much
more computationally inexpensive. Also, the usage of quantile conditioning
instead of mean conditioning makes bQCD more robust to heavy tails as the mean
is more susceptible to outliers than the quantile. Moreover, not making any
assumptions about the parametric class allows bQCD to be applied to a variety
of processes where baseline methods perform significantly poorly when the
assumptions do not hold. The source code of bQCD written in R is available on
this site: https://github.com/tagas/bQCD.

####  3.5.8 JCI

Joint Causal Inference (JCI) leverages prior knowledge by combining data from
multiple datasets from different contexts (Mooij et al. (2020)). Particularly,
JCI is a _causal modeling framework_ rather than a specific algorithm, and it
can be implemented using any causal discovery algorithm that can take into
account some background knowledge. The main idea of JCI is to first, consider
auxiliary context variables that describe the context of each data set, then,
pool all the data from different contexts, including the values of the context
variables, into a single data set, and finally apply standard causal discovery
methods to the pooled data, incorporating appropriate background knowledge on
the causal relationships involving the context variables. The framework is
simple and easily applicable as it deals with latent confounders, cycles (if
the causal discovery method supports this), and various types of interventions
in a unified way. The JCI framework also facilitates analysis of data from
almost arbitrary experimental designs which allow researchers to trade off the
number and complexity of experiments to be done with the reliability of the
analysis for the purpose of causal discovery.

Figure 21: Workflow of the JCI framework.

####  3.5.9 Kg2Causal

Kg2Causal (Sinha & Ramsey (2021)) uses a large-scale general-purpose
biomedical knowledge graph as a prior for data-driven causal discovery. With a
set of observed nodes in a dataset and some relationship edges between the
nodes derived from a knowledge graph, Kg2Causal uses the knowledge graph-
derived edges to guide the data-driven discovery of a causal graph. The main
ideas of this approach are first, mapping each variable in the dataset to a
node in the knowledge graph, and querying relationships between them; next,
extracting a subgraph containing the connected variables with edges between
them; and then this edge set is used as prior knowledge to guide an optimizing
scoring step for inferring the causal graph. An implementation of Kg2Causal is
available at https://github.com/meghasin/Kg2Causal in R language.

####  3.5.10 C-MCMC

Constrained MCMC (C-MCMC) introduces _prior knowledge_ into the _Markov chain
Monte Carlo (MCMC)_ algorithm for structure learning (Xu et al. (2015)).
C-MCMC uses the following _three types of prior knowledge_ : the existence of
parent nodes, absence of parent nodes, and distribution knowledge including
the conditional probability distribution (CPD) of edges and the probability
distribution (PD) of nodes. All prior knowledge should be given by domain
experts. Existence knowledge means that for any node $X_{i}$, a node-set
$pa(X_{i})$ includes all parent nodes of $X_{i}$. The absence of knowledge
means that for a node $X_{i}$, a node-set $pa(X_{i})$ does not include any
parent node of $X_{i}$. PD/CPD knowledge means that the PD of a node and the
CPD of an edge are known. Considering that the prior knowledge may not be
consistent and reliable, a confidence lambda is assigned by domain experts on
each of the prior knowledge that ranges from 0 to 1. This denotes the
certainty level of prior knowledge. A _lambda_ value of 1 indicates very high
confidence in this knowledge.

####  3.5.11 M$eta$-RL

Meta-RL is a _meta-learning algorithm_ in a Reinforcement Learning (RL)
setting where the agent learns to _perform interventions_ to construct a
causal graph (Sauter et al. (2022)). The goal is to be able to use previous
learning experiences during training to generalize in unseen environments.
This approach has some strong assumptions such as i) each environment is
defined by an acyclic SCM, ii) every observable variable can be intervened on,
iii) for each environment in the training set, the underlying SCM is given,
and iv) intervention can be performed on at most one variable at a time. Meta-
RL has two phases: i) Training, and ii) Application. The training phase starts
by randomly choosing an SCM from a set of environments. There are mainly two
sets of actions that an agent performs: _a) interventional actions_ , and _b)
structure actions_. In each step, any one action can be performed on the set
of variables to generate a PDAG. The _agent policy is updated_ via the
_interventional actions_ in each step. However, in case of the structural
actions (e.g. add, delete, or reverse), the agent policy only gets updated at
the end of the training procedure where a reward is sent to the agent. The
reward is computed by comparing the hamming distance of the generated PDAG to
the true causal structure when the training is completed. A _recurrent LSTM
layer_ enables the policy to remember samples from the post-interventional
distributions in the earlier steps. This should help to better identify causal
relations since the results of sequential interventions can be used to
estimate the distribution. Once trained, Meta-RL can then be applied to
environments that have a structure unseen during training. For training, 24
SCMs with 3 observable variables, and 542 SCMs with 4 observable variables
were created. Code to reproduce experiments or run Meta-RL is available at
https://github.com/sa-and/interventional_RL. One limitation of this approach
is that it needs modification in terms of scalability. Also, in real-world
scenarios, every variable might not be accessible for intervention.

Figure 22: Training phase of the Meta-RL algorithm (Sauter et al. (2022)).

####  3.5.12 Tabu search for SEM

Marcoulides et al. (1998) presents an approach to structural equation modeling
(SEM) specification search that makes use of Tabu search, a heuristic
optimization algorithm. Using a neighborhood of the current solution as its
focus, the tabu search technique avoids local optimality by examining the area
around the current solution. To prevent cycling, it assigns recently involved
attributes a tabu status. A number of definitions and parameters, such as the
neighborhood definition and the model selection criterion, are necessary to
implement the Tabu search procedure for SEM specification search. The authors
conclude that Tabu search is a promising strategy for SEM specification search
after demonstrating its efficacy in a number of example analyses.

####  3.5.13 LFCM

Latent Factor Causal Models (LFCMs) (Squires et al. (2022)) perform causal
discovery in the _presence of latent variables_. These models are motivated by
gene regulatory networks. LFCMs work in three stages where they discover: (i)
clusters of observed nodes, (ii) a partial ordering over clusters, and (iii)
finally, the entire structure over both observed and latent nodes. A graph $G$
is called a latent factor causal model (LFCM) if it satisfies the following
conditions: (a) Unique cluster assumption: Each observed node has exactly one
latent parent, (b) Bipartite assumption: There are no edges between pairs of
observed nodes or between pairs of latent nodes, (c) Triple-child assumption:
Each latent node has at least 3 observed children and (d) Double-parent
assumption. The other assumption of LFCMs is that it allows non-exogenous
latent variables. For cluster formation, LFCMs rely on t-separation (Sullivant
et al. (2010)). When two ordered pairs of variables [e.g. ($X_{i}$, $X_{j}$)
and ($X_{u}$, $X_{v}$)] are t-separated, then they belong to the same cluster.
LFCMs are a biologically motivated class of causal models with latent
variables. The limitations of LFCMs include their applicability to only a
linear Gaussian SEM, some major structural restrictions, and that it can fail
when the true graph violates the double parent assumption.

Figure 23: The graph $G$ on the left is a latent factor causal model (LFCM),
and the graph on the right is the latent graph $L(G)$ for $G$ (Squires et al.
(2022)).

There are some other noteworthy methods that are specialized to handle latent
variables or unobserved confounders. Kocaoglu et al. (2017) presented a non-
parametric algorithm for learning a causal graph in the presence of hidden
variables. The study took a stage-by-stage approach, first to learn the
induced graph between observational variables, and then use it to discover the
existence and location of the latent variables. The authors further proposed
an algorithm to discover the latent structure between variables depending on
the adjacency. To identify ancestral relationships and transitive closure of
the causal graph, the algorithm employed a pairwise independence test under
interventions. Then Kocaoglu et al. (2019) addressed causal discovery by
linking conditional independencies from observed data to graphical constraints
using the d-separation criterion. It broadened the application of this
strategy to scenarios involving numerous experimental and observational
distributions. The authors proposed that CIs and d-separation constraints are
just a subset of broader constraints obtained from comparing various
distributions, which is especially useful in the context of do-calculus for
soft interventions. They introduced the notion interventional equivalence
class of causal graphs with latent variables, which linked graphical
structures to groups of interventional distributions that adhered to do-
calculus. Two causal graphs are interventionally equivalent if they produce
identical interventional distributions that can not be distinguished by
invariances.

Sometimes complex systems require the knowledge of both observations and
experiments for recovering the underlying causal relationships. Utilizing both
observational and interventional data from various domains, the authors in Li
et al.  proposed a novel approach for identifying causal structures in semi-
Markovian systems with latent confounders. They made a link between learning
from interventional data within a single domain and learning from
observational data across domains. They introduced the idea of S-Markov, a
property connecting multi-domain distributions to pairs of causal graphs and
interventional targets, to navigate the complexities of observational and
experimental data. A new causal discovery algorithm called S-FCI was
introduced that builds on the S-Markov property and is capable of effectively
learning from a mixture of observational and interventional data from various
domains.

Jaber et al. (2020) integrates soft experimental and observational data to
find the structure in non-Markovian systems with latent variables. They
introduced the idea of $\Psi$-Markov in this context when the intervention
targets were unidentified. This idea links a causal graph G and a list of
interventional targets I to the causal invariances found in both observational
and interventional data distributions. They also introduced a graphical method
for evaluating equivalence between causal graphs with various interventional
targets and an algorithm for learning the equivalence class.

Recently, Bellot et al. (2022) studied structure learning in discrete models
with arbitrary latent dependencies, and proposed a new score based on the
asymptotic expansion of the marginal likelihood to capture both equality and
inequality constraints in observational data. Furthermore, it claims to be the
first score-based method to learn causal models with latent variables.

The different methods discussed in this section so far use a variety of
strategies to perform causal discovery under diverse settings and assumptions.
Therefore, we present a comparative analysis of some of the common methods in
Table LABEL:table-IID-comparison based on their assumptions, output causal
graph, techniques used, advantages, and disadvantages. This comparative
analysis will help readers to find the similar and dissimilar methods, and
also help in deciding which method could be appropriate for performing causal
discovery given the data and its assumptions.

Table 6: Comparison among some causal discovery algorithms for I.I.D. data. Methods | Assumptions |  | Outcome
---
Technique used | Advantages | Disadvantages
PC | Faithfulness, Sufficiency | CPDAG | Conditional Independence (CI) Tests |  | Computationally more
---
feasible for sparse graphs.
| Lacks scalability as less
---
feasible for denser graphs
FCI |  | Causal Markov
---
condition, faithfulness
PAG |  | CI tests, Skeleton
---
finding step same as PC
| Handles latent and
---
selection variables.
| In the worst case, the no. of
---
CI tests performed grows
exponentially with the no.
of variables.
RFCI | CMC, faithfulness | PAG | Conditional Independence tests |  | Faster variant of FCI, uses
---
fewer CI tests, computationally feasible
for high-dimensional sparse graphs
| Performs some additional tests
---
before orienting v-structures
and discriminating paths
GES |  | Decomposable
---
score function
CPDAG | Score-based greedy search |  | Run time faster compared to
---
the constraint-based methods
| Search space can grow
---
exponentially with the growing
no. of variables
FGS | Weak Faithfulness | DAG | Score-based greedy search |  | Faster variant of GES, reduces scoring
---
redundancy, enables parallelization
| Require high power
---
computing to run
RL-BIC |  | Decomposable score
---
function, acyclicity
DAG |  | BIC score-based
---
reinforcement learning search
| Model gets feedback to update
---
/correct its search strategy
| Scalable only up to a few
---
(around 30) variables
Triplet A* | Acyclicity | DAG |  | A* search combined
---
with BIC score
| Can handle both linear Gaussian &
---
non-Gaussian networks, scales up to
more than 60 variables
Has complexity issues
KCRL |  | Decomposable score-function,
---
unbiased prior knowledge
DAG | Score-based RL search strategy |  | Considers prior
---
knowledge constraints
Lacks scalability
LiNGAM |  | Linear DGP, no unobserved
---
confounders, non-Gaussian noises
DAG |  | FCM-based, uses independent
---
component analysis
| Determines the direction of
---
every causal arrow, does not r
require the faithfulness assumption.
| Estimated results may
---
vary in case of mixed data
(categorical values)
SAM | Acyclicity | DAG | Searches for an FCM using GANs | Good for sparse causal graphs |  | Not suitable for
---
dense graphs
Tabu search |  | Finite search space,
---
well-defined objective function
DAG |  | Iteratively modifies the model
---
and evaluates its fit to the data
by examining a neighborhood of the
current solution
| Considers avoiding local optima
---
and also, cycle avoidance
| Quality of the results depends
---
on the quality of the initial
solution and the choice
of parameters
CAM | Sufficiency, Acyclicity | DAG |  | FCM-based, uses sparse
---
regression techniques and decouples
order search among the variables
| Establishes consistency of the
---
maximum likelihood estimator
for low and high-dimensional cases
| Faces performance and
---
computational challenges
as the number of
variables increases
NOTEARS |  | Acyclicity, linear dependence
---
between variables
DAG |  | Models data as a SEM,
---
uses a regularized score-function
with a least-square loss
| Simple method,
---
easy to implement
| Works well mostly for
---
continuous data
GAE |  | Structure learning under
---
additive noise models
DAG |  | Gradient-based approach that uses a
---
graph autoencoder framework
| Can scale up to 100 nodes,
---
Lower run time, easily applicable
to vector-valued variables, & handles
non-linear relations well
| May not work well for
---
linear causal relations
DAG-GNN | Faithfulness, sufficiency | DAG |  | Gradient-based, uses blackbox
---
stochastic optimization solvers
to solve the sub problem of
maximizing the ELBO
| Can handle both discrete and
---
vector-valued variables, and is
capable of capturing complex
nonlinear mappings
| Assumes acyclic causal
---
relationships, which may not
always be the case
in real-world scenarios
MMHC |  | Faithfulness, sufficiency,
---
decomposable score-function
DAG |  | Hybrid: uses both score and
---
constraint-based techniques such as
MMPC to initially learn the skeleton of
the network, then uses greedy
Bayesian hill-climbing search
| Good computational efficiency,
---
scalability, and applicable for
high-dimensional settings
| Scales up better only with
---
large number of samples
FRITL |  | Faithfulness, non-Gaussianity,
---
latent confounder assumption
PAG |  | Hybrid, uses FCI, Triad
---
condition and ICA
| Works in presence or
---
absence of latent confounders
| Not generalizable in
---
non-linear gaussian cases
Kg2Causal |  | Presence of a knowledge graph
---
DAG |  | Prior knowledge constraints
---
added in a score-based method
| Leverages information from
---
existing literature or domain
Requires a knowledge graph

##  4 Causal Discovery Algorithms for Time Series Data

Time series data arise when observations are collected over a period of time.
So far the methods that we have discussed are specialized for causal discovery
from I.I.D. or time-independent data. However, often, real-world data in
different domains can be a time series (non-I.I.D. data). For this type of
data, there are different specialized causal discovery approaches based on CI
testing, SEM/FCMs, Granger causality (Granger (1969)), or deep neural
networks. In this section, first, we provide a brief introduction to some of
the common terminologies related to time-series data and temporal causal
discovery. Then, we discuss the notable causal discovery approaches for time-
series data.

######  Definition 7 (Time Series Data)

Time series data is a collection of observations measured over consistent
intervals of time. The observation of a time series variable $X^{j}$ at time
$t$ is denoted by $X^{j}_{t}$.

Examples of time series data include retail sales, stock prices, climate data,
heart rate of patients, brain activity recordings, temperature readings, etc.
Any time series data may have the following _properties_ :

  1. i

_Trend :_ When the data show a long-term rise or fall, a trend is present.
Such long-term increases or decreases in the data might not be always linear.
The trend is also referred to as changing direction when it might switch from
an upward trend to a downward trend.

  2. ii

_Seasonality :_ It refers to the seasonal characteristics of time series data.
Seasonality exists when the data regularly fluctuates based on different time
spans (e.g. daily/weekly/ monthly/quarterly/yearly). An example is temperature
data, where it is mostly observed that the temperature is higher in the
summer, and lower in the winter. Any analysis related to time series usually
takes advantage of the seasonality in data to develop more robust models.

  3. iii

Autocorrelation: Autocorrelation or self-correlation is the degree of
similarity between a given time series and a lagged version of itself over
successive time intervals. Time series data is usually autocorrelated i.e.,
the past influences the present and future (Lawton et al. (2001)).

  4. iv

Stationarity & Non-stationarity: Stationarity means that the joint probability
distribution of the stochastic process does not change when shifted in time. A
time series is stationary if it has causal links such that for variables
$X^{i}$ and $X^{j}$, if $X^{i}$ → $X^{j}$ at any timestamp $t$, then $X^{i}$ →
$X^{j}$ also holds for all $t^{\prime}$ $\not=$ $t$. This condition does not
hold for a non-stationary time series where $X^{i}$ → $X^{j}$ at a particular
time $t$ need not necessarily be true at any other time stamp $t^{\prime}$.

Let $X^{j}_{1:t}$ = {$X^{1}_{1:t}$, $X^{2}_{1:t}$, …, $X^{n}_{1:t}$} be a
multivariate time series with $n$ variables and $t$ time steps. At any
particular timestamp $t$, the state of the $n$ variables can be represented as
$X^{j}_{t}$ = {$X^{1}_{t}$, $X^{2}_{t}$, …, $X^{n}_{t}$}. The past of a
variable $X^{j}_{t}$ is denoted by $X^{j}_{1:t-1}$. The parent set of a
variable includes all the nodes with an edge towards it. The goal of any
temporal causal discovery approach is to discover the causal relationships
between the time series variables. Any time series causal graph may have the
following _types of causal relationships/edges_ : (i) Instantaneous edges, and
(ii) Lagged edges.

Figure 24: Types of causal relationships: Instantaneous edges (red), Lagged
edges (blue), and Changing modules (green).

######  Definition 8 (Instantaneous Causal Effect)

When the delay between cause and effect is 0 timesteps, i.e. causal effects
are of the form $X^{i}_{t}$ $\rightarrow$ $X^{j}_{t}$ or $X^{i}_{t}$
$\rightarrow$ $X^{i}_{t}$ (self-causation), then it is known as an
instantaneous or contemporaneous causal relationship/effect (Nauta et al.
(2019)).

######  Definition 9 (Lagged Causal Effect)

When the delay between cause and effect is at least 1 or more timesteps (i.e.
causal effects of the form $X^{i}_{t-}$ $\rightarrow$ $X^{j}_{t}$ or
$X^{i}_{t-}$ $\rightarrow$ $X^{i}_{t}$), then it is known as a lagged causal
relationship/effect. That is, a lagged causal effect occurs when a variable
causes another variable or itself with a time lag = 1 or more.

In Figure 24, the red-colored edges represent the instantaneous causal effect
(relationships among the variables at the same time step), and the blue edges
represent the lagged causal effect. The green edges represent a special form
of temporal causal relationships known as the _changing modules_ (CM). The CMs
represent the direct effect of a time stamp on a variable (e.g. $t\rightarrow
X^{1}_{t}$ in Figure 24). Details on CM are available in Ferdous et al.
(2023b).

The causal graphs produced by different temporal causal discovery algorithms
vary based on the details of the relationships they represent. Any temporal
causal discovery algorithm may produce any of the following two types of
temporal causal graph as its outcome: a full-time causal graph or a summary
causal graph (see Figure 25).

######  Definition 10 (Full-time Causal Graph)

A full-time causal graph represents both the instantaneous ($X^{i}_{t}$
$\rightarrow$ $X^{j}_{t}$ or $X^{i}_{t}$ $\rightarrow$ $X^{i}_{t}$) and time-
lagged ($X^{i}_{t-}$ $\rightarrow$ $X^{j}_{t}$ or $X^{i}_{t-}$ $\rightarrow$
$X^{i}_{t}$) causal edges where the lag between a cause and effect is
specified in the graph. A full-time casual graph may sometimes present the
changing modules as well.

Figure 25: Full-time causal graph (Left) & Summary causal graph (Right)

Figure 25 (left) represents a full-time causal graph where both instantaneous
relations (e.g. $X^{1}_{t}$ $\rightarrow$ $X^{2}_{t}$, $X^{1}_{t-1}$
$\rightarrow$ $X^{2}_{t-1}$), and lagged relations (e.g. $X^{2}_{t-1}$
$\rightarrow$ $X^{3}_{t}$, $X^{3}_{t-1}$ $\rightarrow$ $X^{3}_{t}$) among the
variables are depicted.

######  Definition 11 (Summary Causal Graph)

A summary causal graph is a reduced version of a full-time causal graph where
each lagged node represents the entire past ($X^{j}_{t-}$) of its
corresponding instantaneous node ($X^{j}_{t}$), and the exact time lag between
the cause and effect is not specified in the graph.

In the following subsections, we describe briefly some of the notable causal
discovery algorithms that focus on time series data. Figure 26 presents a
taxonomy of some of the discussed approaches.

forked edges, for tree=draw,align=left,l=1.7cm,edge=-latex [Causal Discovery
Algorithms
(for time series data) [Contraint-based [tsFCI (4.1.1)
PCMCI (4.1.2)
LPCMCI∗ (4.1.4)
CDANs (4.1.6)] ] [FCM-based [VarLiNGAM (4.2.1)
TiMINo (4.2.2)] ] [Gradient-based [DYNOTEARS⋄ (4.3.1)
NTS-NOTEARS⋄∗ (4.3.2)] ] [Granger Causality-based [GVAR (4.4.1)
NAVAR (4.4.2)
ACD (4.4.3)] ] [Miscellaneous [oCSE (4.5.1)
TCDF (4.5.2)
NBCB (4.5.3)
PCTMI (4.5.4)
] ] ]

Figure 26: Taxonomy of some of the discussed causal discovery approaches for
time series data. The approaches are classified based on their core
contribution or the primary strategy they adopt for causal structure recovery.
The approaches that can leverage prior knowledge are marked by an $\ast$
symbol. Some of the gradient-based approaches that use a score function are
indicated by a $\diamond$ symbol. They are primarily classified as such as
they use gradient descent for optimization. However, they can be a score-based
method too as they compute data likelihood scores on the way.

###  4.1 Constraint-based

####  4.1.1 $ts$FCI

The algorithm time series FCI or tsFCI (Entner & Hoyer (2010)) adapts the Fast
Causal Inference (Spirtes et al. (2000a)) algorithm (developed for the causal
analysis of non-temporal variables) to infer causal relationships from time
series data. It works in two phases: (i) an adjacency phase, and (ii) an
orientation phase. It makes use of temporal priority and consistency
throughout time to orient edges and restrict conditioning sets. It provides a
window causal graph, and an advantage is that it _can detect lagged hidden
confounders_. However, a disadvantage is that it cannot model cyclic
contemporaneous causation, and also instantaneous relationships. A code
package that implements tsFCI is available at
https://sites.google.com/site/dorisentner/publications/tsfci.

####  4.1.2 PCMCI

A problem with large-scale time series data is that although adding more
variables makes causal analysis more interpretable, if the additional
variables don’t have a significant effect on the causal model, this, in turn,
makes the analysis less powerful, and original causal relations may also be
overlooked. Moreover, at large dimensions, certain nonlinear tests even lose
their ability to limit false positive rates (FPRs). Runge et al. (2019)
proposed a two-stage algorithm PCMCI that can overcome this problem. In
Step-1, the model selects conditions using $PC_{1}$ (a variant of the skeleton
discovery part of the PC algorithm) to remove irrelevant variables which solve
the issue of low power in the causal discovery process. In Step-2, the
momentary conditional independence (MCI) test is used which helps to reduce
the FPR even when the data is highly correlated. The MCI test measures if two
variables are independent or not given their parent sets (see Equation 13).

$X^{i}_{t-\tau}\perp\\!\\!\\!\perp
X^{j}_{t}|P_{A}(X^{j}_{t}),P_{A}(X^{i}_{t-\tau})$ (13)

PCMCI assumes that the data is stationary, has time-lagged dependencies, and
also assumes causal sufficiency. Even when the stationary assumption is
violated (probably by obvious confounders), PCMCI still provides a more robust
performance than Lasso regression or the PC algorithm. However, for highly
predictable systems where little new information is produced at each time
step, PCMCI is not a good fit. Python implementation of PCMCI is available in
the Tigramite package (https://github.com/jakobrunge/tigramite).

Figure 27: Steps involved in the PCMCI method for time series causal
discovery.

####  4.1.3 PCMCI+

PCMCI+ (Runge (2020)) is an extension of the PCMCI algorithm to discover
contemporary or instantaneous causal links. PCMCI+ also assumes causal
sufficiency like the PCMCI algorithm. It is also a two-stage algorithm where
in the first stage, irrelevant edges from the causal model are eliminated.
Unlike PCMCI, the edges are removed separately for lagged and contemporary
conditioning sets where the contemporary phase employs more CI tests than the
lagged phase. In the second stage, PCMCI+ employs the notion of momentary
conditional independence (MCI) to improve the selection of conditioning sets
for the various CI tests, improving their autocorrelation calibration, and
boosting their detection power. The results show that when there is high
autocorrelation in the data, PCMCI+ can achieve better performance in terms of
higher recall, lower false positives, and faster execution compared to the PC
algorithm. For lower autocorrelation, PCMCI+ performs almost similarly to PC.
Implementation of PCMCI+ is also available in the Tigramite package
(https://github.com/jakobrunge/tigramite).

####  4.1.4 LPCMCI

Latent PCMCI (LPCMCI) is a constraint-based causal discovery algorithm to
determine causal relationships from large-scale time series data (Gerhardus &
Runge (2020)). This is another extension of the PCMCI algorithm as it can
discover causal relationships even in the presence of latent confounders.
Moreover, it gives the flexibility to use the model when the data is linear or
nonlinear, and also when the data has lagged or contemporary conditioning
sets. The authors identified that when the CI tests have a low effect size,
existing techniques like FCI suffer from low recall in the presence of
autocorrelation. They demonstrated that this issue can be solved by including
causal parents in the conditioning sets. By utilizing the orientation rules,
these parents can be identified as early as in the edge removal stage. The
results show that the proposed LPCMCI method can achieve higher recall than
the baseline model SVAR-FCI. However, LPCMCI cannot differentiate all members
of the Markov class, and also, when the faithfulness assumption doesn’t hold,
LPCMCI might lead to an incorrect conclusion. Along with PCMCI and PCMCI+, the
Python code of LPCMCI is also available in the Tigramite GitHub package.

####  4.1.5 CD-NOD

Many existing approaches assume that the causal model is static, and
therefore, there will be a fixed joint distribution of the observed data.
However, these methods fail when the underlying data changes over time, and
causal parameters vary during the period. Huang et al. (2020) proposed a
causal discovery method that assumes that the parameter of the causal model
can change over time or different datasets, and they named the method CD-NOD,
Constraint-based Causal Discovery from Heterogeneous/Nonstationary Data. The
proposed method can determine causal direction by taking advantage of
distribution shifts, and these distribution changes, in the presence of
stationary confounders, are helpful for causal discovery. The distribution
shifts can be either time or domain indexes and are denoted by a surrogate
variable $C$. Broadly, CD-NOD has two phases where in the first phase it
recovers the causal skeleton $S_{G}$, and in the second phase it orients the
edges as per some orientation rules. Given that the causal model offers a
concise summary of how the joint distribution changes, they demonstrated that
distribution shift contains important information for causal discovery.
Recently, researchers discovered that this idea could help solve machine
learning problems of domain adaptation and forecasting in nonstationary
situations (Schölkopf et al. (2012); Zhang et al. (2013)). The conducted
experiments in this study demonstrate the changes of causal influence between
the different states of brain functions, and the empirical results show that
CD-NOD has improved precision and F1 score. However, they didn’t consider that
the causal directions might flip, and the power of conditional independence
tests might reduce because of the distribution shifts. The algorithm’s source
code is available in the following link: https://github.com/Biwei-
Huang/Causal-Discovery-from-Nonstationary-Heterogeneous-Data.

Figure 28: Illustration of CD-NOD’s phase-1.

####  4.1.6 CDANs

Ferdous et al. (2023b) introduces a constraint-based causal discovery approach
called CDANs for autocorrelated and non-stationary time series data that
handles high dimensionality issues. The method identifies both lagged and
instantaneous causal edges along with changing modules that vary over time. By
optimizing the conditioning sets in a constraint-based search, and also
considering lagged parents instead of conditioning on the entire past, it
tries to address the high dimensionality problem. CDANs first detect the
lagged adjacencies, then identify the changing modules and instantaneous
adjacencies, and finally determine the causal direction. The code to implement
this method is available at https://github.com/hferdous/CDANs. An extended
version of this study is presented in Ferdous et al. (2023a), where the method
called eCDANs is introduced that is capable of detecting lagged and
instantaneous causal relationships along with temporal changes. The method
eCDANs addresses high dimensionality by optimizing the conditioning sets while
conducting CI tests and identifies the changes in causal relations by
introducing a proxy variable to represent time dependency.

###  4.2 Functional Causal Model (FCM)-based

####  4.2.1 VarLiNGAM

VarLiNGAM (Hyvärinen et al. (2010)) combines the non-Gaussian instantaneous
models with autoregressive models and shows that a non-Gaussian model is
identifiable without prior knowledge of network structure. It estimates both
instantaneous and lagged causal effects in models that are an example of
structural vector autoregressive (SVAR) models. These models are a combination
of structural equation models (SEM) and vector autoregressive (VAR) models.
VarLiNGAM also shows that taking instantaneous influences into account can
change the values of the time-lagged coefficients to a great extent. Thus,
neglecting instantaneous influences can lead to misleading interpretations of
causal effects. It also assesses the significance of the estimated causal
relations. An implementation of this method is available at:
https://lingam.readthedocs.io/en/latest/tutorial/var.html.

####  4.2.2 T$i$MINo

Time-series Models with Independent Noise (TiMINo) (Peters et al. (2013))
studies a class of restricted structural equation models (SEMs) for time-
series data that include nonlinear and instantaneous effects. It assumes
$X_{t}$ to be a function of all direct causes and some noise variable, the
collection of which is supposed to be jointly independent. The algorithm is
based on unconditional independence tests and is applicable to multivariate,
linear, nonlinear, and instantaneous interactions. If the model assumptions
are not satisfied by the data, TiMINo remains mostly undecided instead of
making wrong causal decisions. While methods like Granger causality are built
on the asymmetry of time direction, TiMINo additionally takes into account
identifiability emerging from restricted SEMs. This leads to a straightforward
way of dealing with unknown time delays in different time series. An
implementation of TiMINo is available in this repository:
https://github.com/ckassaad/causal_discovery_for_time_series.

###  4.3 Gradient-based

####  4.3.1 DYNOTEARS

Pamfil et al. (2020) proposed the Dynamic NOTEARS (DYNOTEARS) which is a
structure learning approach for dynamic data that simultaneously estimates
contemporaneous (intra-slice) and time-lagged (inter-slice) relationships
between variables in a time-series. DYNOTEARS revolves around minimizing a
penalized loss subject to an acyclicity constraint. The optimization finds the
conditional dependencies that are best supported by the data. It leverages
insight from the approach NOTEARS (Zheng et al. (2018)) which uses an
algebraic characterization of acyclicity in directed graphs for static data.
The assumptions made by DYNOTEARS include that the structure of the network is
fixed through time, and is identical for all time series in the data. This
approach is scalable to high-dimensional datasets. An implementation of this
approach is available in the CausalNex library
(https://github.com/quantumblacklabs/causalnex), and also at
https://github.com/ckassaad/causal_discovery_for_time_series.

####  4.3.2 NTS-NOTEARS

NTS-NOTEARS (Sun et al. (2021)) is a causal discovery method for time series
data that uses 1-D convolutional neural networks (CNNs) to capture linear,
nonlinear, lagged, and instantaneous relations among variables in a time
series data along with ensuring the acyclicity property of a DAG. It extends
the continuous optimization-based approach NOTEARS for learning nonparametric
instantaneous DAGs, and adapts the acyclicity constraint from that approach.
It assumes that there are no latent confounders in the data, and the
underlying data-generating process is fixed and stationary over time. NTS-
NOTEARS is faster than other constraint-based methods because of the use of
nonlinear conditional independence tests. It incorporates prior knowledge into
the learning process to promote the use of optimization constraints on
convolutional layers for better casual discovery. Its implementation is
available at: https://github.com/xiangyu-sun-789/NTS-NOTEARS/.

###  4.4 Granger Causality (GC)-based

Granger (1969) investigated the causal relationships between the variables in
a time series data which is known as Granger Causality (GC). It is based on
the basic assumption that _causes precede their effects_. The author defines
GC as follows: A time series variable $X^{i}$ causes $X^{j}$, if the
probability of $X^{j}$ conditional on its own past, and the past of $X^{i}$
(besides the set of the available information) does not equal the probability
of $X^{j}$ conditional on its own past alone. The GC test can’t be performed
directly on non-stationary data. The non-stationary data needs to be
transformed into stationary data by differencing it, either using first-order
or second-order differencing. Granger Causality can be used when there are no
latent confounders, and also, no instantaneous effects exist, i.e., no
variable causes another variable at the same time stamp.

####  4.4.1 GVAR

Generalized Vector AutoRegression (GVAR) (Marcinkevičs & Vogt (2021)) is a
framework for inferring multivariate Granger causality under nonlinear
dynamics based on autoregressive modeling with self-explaining neural
networks. It allows the detection of signs of Granger-causal effects and
inspection of their variability over time in addition to relational inference.
It focuses on two aspects: first, inferring Granger-causal relationships in
multivariate time series under nonlinear dynamics, and second, inferring signs
of Granger-causal relationships. A reproducible code of the approach is
available at: https://github.com/i6092467/GVAR.

####  4.4.2 NAVAR

Bussmann et al. (2021) proposed the approach Neural Additive Vector
AutoRegression (NAVAR) which is a causal discovery approach for capturing
nonlinear relationships using _neural networks_. It is particularly trained
using deep neural networks that extract the (additive) Granger causal
influences from the time evolution in multivariate time series. NAVAR assumes
an additive structure where the predictions depend linearly on independent
nonlinear functions of the individual input variables. These nonlinear
functions are modeled using neural networks. The additive structure of NAVAR
allows scoring and ranking the causal relationships. Currently, NAVAR is
implemented with MLPs and LSTMs as the backbone using Python which is
available at: https://github.com/bartbussmann/NAVAR. However, more complex
architectures such as dilated CNNs and transformers can also be used to model
NAVAR.

####  4.4.3 ACD

Most causal discovery algorithms applied for time-series analysis find a
causal graph for the data, and then refit the model whenever new samples do
not fit with the underlying causal graph. But in many cases, samples share
connections among them, for example, the brain activity of different regions
at different times. When the algorithms fit a new model, this dynamic nature
between the samples is lost, and can no longer identify the actual causal
relation. To solve this problem, Löwe et al. (2022) proposed the Amortized
Causal Discovery (ACD) technique which can identify the causal relations when
samples are from different causal graphs but share common dynamics. ACD
consists of an encoder and a decoder. The encoder predicts the causal graph’s