# Differentiable Invariant Causal Discovery

Yu Wang, An Zhang, Xiang Wang, Yancheng Yuan, Xiangnan He, Tat-Seng Chua  Y.
Wang, X. Wang and X. He are with University of Science and Technology of
China. Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>A. Zhang, T. Chua are with Sea-NExt Joint Lab from National University of
Singapore. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>Y. Yuan is with The
Hong Kong Polytechnic University. Email<EMAIL_ADDRESS>Correspondonce to A. Zhang.

###### Abstract

Learning causal structure from observational data is a fundamental challenge
in machine learning. However, the majority of commonly used differentiable
causal discovery methods are non-identifiable, turning this problem into a
continuous optimization task prone to data biases. In many real-life
situations, data is collected from different environments, in which the
functional relations remain consistent across environments, while the
distribution of additive noises may vary. This paper proposes Differentiable
Invariant Causal Discovery (DICD), utilizing the multi-environment information
based on a differentiable framework to avoid learning spurious edges and wrong
causal directions. Specifically, DICD aims to discover the environment-
invariant causation while removing the environment-dependent correlation. We
further formulate the constraint that enforces the target structure equation
model to maintain optimal across the environments. Theoretical guarantees for
the identifiability of proposed DICD are provided under mild conditions with
enough environments. Extensive experiments on synthetic and real-world
datasets verify that DICD outperforms state-of-the-art causal discovery
methods up to 36% in SHD. Our code will be open-sourced.

###### Index Terms:

Causal Discovery, Invariant Learning, Causal Structure Learning, Causal Graph
Learning.

##  1 Introduction

Causal discovery (CD) is a fundamental problem in a variety of tasks, such as
understanding the generation process of data [1], and probing explainability
of models [2, 3]. It has tremendous impacts on various domains, like biology
[4, 5], and finance [6]. CD aims to learn the causal structure among a set of
variables from the observational data and represent the structure as a
directed acyclic graph (DAG). The acyclicity constraint frames CD as the
combinatorial optimization of discrete edges, which however, is NP-hard.
Recently, leading CD solutions [1, 7, 8, 9] gracefully convert the DAG
learning into the continuous optimization task. Specifically, the idea
stemming from NOTEARS [1] is to build a scoring function upon the adjacency
matrix over the variables and find an equivalent continuous constraint on
acyclicity.

NOTEARS [1] inspired the development of numerous differentiable causal
discovery algorithms that use gradient descent to find the optimal causal
graph [10, 11, 12]. When compared to traditional constraint-based causal
discovery, these methods have demonstrated superior performance and efficiency
in uncovering the true graph with a large amount of data. However, most of
them follow the paradigm of empirical risk minimization (ERM) [13] — first
imposing the scored DAG on the observations to reconstruct, and then
minimizing the empirical risks between the observational and reconstructed
data, so as to optimize the DAG. Despite the promising performance, we argue
that ERM is prone to capture data biases or shortcut [14, 15, 16], thus
derailing the structure learning of DAG. Specifically, the observations of
variables are often marred by some spurious correlations, such as the
annotator or selection biases in the data acquisition pipeline [14], thus
posing undesired entanglements of variables. ERM easily latches on these
correlations [14, 17, 15] to refine the DAG structure, as the following
example illustrates.

Consider the ground-truth DAG in Table I as the target being reconstructed by
the CD solutions, where each edge denotes a causal relationship between two
variable nodes. $X$ consists of three variables, which have different
relationships with $Y$: $A$ determines $Y$, $B$ is influenced by $Y$, and $C$
is irrelevant to $Y$. While the CD solutions are designed to identify the
causation edges, ERM does not need to learn the correct DAG to reach a low
reconstruct loss for fitting the observations. For example, instead of looking
at the true causation $B\leftarrow Y$, it is easy to capture the statistical
shortcut edge $B\rightarrow Y$, since $B$ is strongly correlated with $Y$ and
the reconstructed DAGs achieve even lower losses as compared to the ground-
truth DAG. This common problem is consistent with over-reconstruction [18] but
remains largely unexplored.

TABLE I: Examples that NOTEARS would find the wrong causal graph while multi-environment settings can help identify the true graph. We present the DAGs derived from minimal reconstruction loss in different environments. The graph on the left denotes the ground truth. The edges in red are wrong or spurious, while the dash one represents the corresponding coefficient is zero. We list the reconstruction loss value following with the triplets in brackets representing coefficients of $Y$ with $A,B$ and $C$. The loss values and parameters of the potential sub-optimal causal structures learned by NOTEARS in a single environment are highlighted in red. Detailed data generating processes are as follows. $X\sim\mathcal{N}(0,1)$, $A=X+z_{A}(\sim\mathcal{N}(0,1))$, $B=X+Y/2+z_{B}^{e}(\sim\mathcal{N}(0,(\sigma_{B}^{e})^{2}))$, $C=X/2+z_{C}^{e}(\sim\mathcal{N}(0,(\sigma_{C}^{e})^{2}))$, $Y=A/4+\epsilon_{Y}(\sim N(0,1))$. $\sigma_{B}^{e}$ and $\sigma_{C}^{e}$ varies across three environments, $e_{1},e_{2}$ and $e_{3}$, with $(\sigma_{B}^{e})^{2}=(\sigma_{C}^{e})^{2}=1,2,4$, respectively.  BCAXY |  | BCAXY | BCAXY | BCAXY
---|---|---|---|---
$e_{1}$ | 4.57 (0.24, 1.14, -1.32) | 5.07 (0.24, 0.32, 0.00) | 5.07 (0.24, 0.50, 0.00)
$e_{2}$ | 6.59 (0.24, 1.09, -1.19) | 7.14 (0.24, 0.23, 0.00) | 7.07 (0.24, 0.50, 0.00)
| $e_{3}$ | 10.60 (0.24, 1.07, -1.12) | 11.21 (0.24, 0.15, 0.00) | 11.07 (0.24, 0.50, 0.00)
|  | BCAXY | BCAXY | BCAXY
5.00 (0.25, 0.50, 0.00) | $e_{1}$ | 5.05 (0.24, 0.50, 0.10) | 4.57 (-0.23, 1.38, -1.54) | 5.12 (0.07, 0.29, 0.00)
7.00 (0.25, 0.50, 0.00) | $e_{2}$ | 7.06 (0.24, 0.50, 0.06) | 6.59 (-0.24, 1.36, -1.44) | 7.17 (0.14, 0.18, 0.00)
11.00 (0.25, 0.50, 0.00) | $e_{3}$ | 11.06 (0.24, 0.50, 0.03) | 10.59 (-0.24, 1.35, -1.39) | 11.21 (0.18, 0.11, 0.00)

In this study, we aim to design a paradigm that distinguishes the causation
edges from spuriously-correlated edges and obtains the faithful DAG. Although
learning causal structure from observational data is challenging, we draw
inspiration from invariant learning [14, 17, 15] and approximate the task by
searching the edges with invariant structural equations in multi-environment
settings. Across different environments, only factual causation edges remain
invariant, while edges with wrong causal directions hardly remain stable,
according to the assumptions in invariant learning [14]. Though reminiscent of
past ideas - _e.g.,_ causal structure learning in multi-domain, from
heterogeneous/nonstationary data [19, 20, 21, 22, 23, 24, 25] \- these methods
are only applicable in a linear system or restricted to conditional
independence tests, which suffer from high computation complexity as the
number of variables increases. To achieve effective differentiable causal
discovery in both linear and nonlinear scenarios, we reconsider Table I and
additionally exhibit the multi-environment information, which impacts the
distributions of additive noises. Applying ERM on these environments
separately results in different functions when there exist spurious
correlations in the causal graph. This inspires us to exclude such unstable
edges towards a robust DAG across environments.

Towards this end, we propose differentiable invariant causal discovery (DICD),
a novel scheme of DAG structure learning that incorporates the idea of
invariant learning — that is, learning an invariant DAG with the piece of
environment information. Specifically, a DAG generator module learns to
generate the environment-aware DAGs from individual environments. For each
DAG, the structure equation model (SEM) [26] describes the functions of
learned causation edges. We then exploit invariant risk minimization (IRM)
[14] to encourage the SEMs to be optimal across all environments, obtaining
DAGs regardless of environment changes. On synthetic and real-world datasets,
extensive experiments demonstrate the effectiveness of DICD to surpass current
state-of-the-art CD solutions.

Our main contributions are:

  * •

To the best of our knowledge, we are among the first class to adopt the
environment information into the differentiable causal discovery framework.

  * •

We propose a novel causal discovery solution, DICD, to incorporate invariant
learning in both linear and nonlinear settings. Experimental results on
synthetic and real-world datasets demonstrate that DICD could significantly
better reveal true correlations and eliminate the spurious ones compared to
prevalent methods.

  * •

We provide theoretical guarantees for the identifiability of proposed DICD in
linear systems under mild conditions, given certain assumptions about the
environments.

##  2 Preliminary

We first make some necessary and reasonable assumptions. Then, we give the
task formulation of causal discovery (CD). After that, we introduce the
continuous optimization paradigm via empirical risk minimization (ERM), as
well as the linear and nonlinear solutions. We list the notations mentioned in
our paper in Appendix Notations

Assumptions. (i) The structural equations are invariant across different
environments; (ii) The system is causally sufficient; (iii) Observational data
is generated from a structural equation model with independent additive noise;
(iv) The distribution shift of additive noise appears in different
environments.

The assumptions (i) - (iii) are crucial to our method and provide key insights
on how causal structure can be identified from observational data. Assumption
(i) states that the invariance principle is held across environments.
Assumption (ii) ensures that no hidden confounders exist in the system. In
other words, there is no systematic bias induced by hidden confounders.
Assumption (iii) implies that only independent additive noise is considered in
our paper, which is commonly used in causal discovery. Assumption (iv) defines
the multiple environments, which is consistent with the definition of
environment assumed in [20, 22] with heterogeneous data [21] or interventional
data [19]. A formal definition of different environments is given in
Definition 3.1. Technically speaking, assumption (iv) presents both challenges
and opportunities for causal discovery and provides a sufficient condition
that the causal structure becomes identifiable. A more formal Theorem 4.1 on
this argument will be presented later in Section 4.

Task Formulation of CD. Let
$\textbf{X}={[\textbf{x}_{1}|\cdots|\textbf{x}_{n}]}^{\top}\in\mathbb{R}^{n\times
d}$ denote the $n$ observational data of $d$ variables, which are generated
from a target directed acyclic graph (DAG). The target DAG is
$(\mathcal{V},\mathcal{D})$, where $\mathcal{V}$ represents the set of node
variables, denoted as $\\{X_{1},\cdots,X_{d}\\}$. And $\mathcal{D}$ is the set
of cause-effect edges between variables. The observational data is assumed to
be generated from the following SEM:

$X_{j}=F_{j}(Pa(X_{j}))+z_{j},j\in\\{1,\cdots,d\\},$ (1)

where $X_{j}$ is the $j$-th node variable, $F_{j}$ is the causal structure
function, $Pa(X_{j})$ is the set of the parents of $X_{j}$, and $z_{j}$ refers
to the additive noise with variance $\sigma_{j}^{2}$. Without loss of
generality, we assume that the noises are zero-mean. In real-life settings,
since the dataset may be obtained from various environments, the distribution
of additive noises $z_{j}$ may differ across different environments. The
causal structure functions $F_{j}$, on the other hand, are generally
invariant. The goal of CD is to learn a DAG to reconstruct the observations X.
The acyclicity restriction of DAG is the fundamental obstacle, because it
frames DAG learning as a NP-hard combinatorial optimization task.

Common Paradigm of ERM. Popular differentiable CD solutions, _e.g.,_ NOTEARS
[1] and its follow-up studies [10, 12], convert DAG learning into a continuous
optimization process to overcome the obstacle of the combinatorial
optimization problem. The primary idea is to build a scoring function upon the
adjacency matrix of variables and discover an equivalent continuous constraint
on acyclicity. To optimize the scoring function, they mostly adopt the
paradigm of ERM to minimize the empirical risks between the observational and
reconstructed data as:

$\displaystyle\min_{f}\mathbf{\mathop{\mathcal{L}}}(f)=\frac{1}{n}\sum_{i=1}^{n}l(\textbf{x}_{i},f(\textbf{x}_{i}))~{}~{}\text{s.t.}~{}~{}\mathcal{G}(f)\in\text{DAG},$

where $\textbf{x}_{i}$ is the $i$-th sample in the dataset. $l(\cdot,\cdot)$
is the reconstruction loss function, _i.e.,_ squared loss or negative log-
likelihood. $f=(f_{1},\cdots,f_{d})$ formulates the estimated structure
function [26] of variables, where $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$
is the estimated structure function of node variable $X_{i}$, and thus
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the structure function for all
nodes. $f_{i}(X_{1},\cdots,X_{d})$ is dependent on $X_{j}$, if
$X_{j}\in\text{Pa}(X_{i})$. The score-based solutions seek to learn $f$
conditioning on the DAG constraint: $\mathcal{G}(f)\in\text{DAG}$, where
$\mathcal{G}(f)$ refers to the graph corresponding to $f$.

Following prior studies [1, 10], we use the matrix
$\textbf{W}\in\mathbb{R}^{d\times d}$ to encode the graph $\mathcal{G}(f)$,
where each element $[\textbf{W}]_{ij}\neq 0$ indicates the existence of edge
$X_{i}\rightarrow X_{j}$. Then we define the matrix W as:

$\displaystyle[\textbf{W}(f)]_{ij}:=|\partial_{i}f_{j}|,$ (2)

where $\partial_{i}f_{j}=\frac{\partial f_{j}(X_{1},\cdots,X_{d})}{\partial
X_{i}}$ is $f_{j}$’s partial derivative _w.r.t._ $X_{i}$. Thus, $f_{j}$ does
not depend on $X_{i}$ if and only if $|\partial_{i}f_{j}|=0$. With Equation
(2), we can exploit linear and nonlinear functions $f$ to characterize
acyclicity in linear and nonlinear SEMs, respectively.

Linear SEM. In the case of linear SEM [1], we formulate $f(\textbf{X})$ as a
linear matrix multiplication: $f(\textbf{X})=\textbf{X}\textbf{A}$ (with each
instance $f(\textbf{x}_{i})={\textbf{A}}^{\top}\textbf{x}_{i}$,
$i\in\\{1,\cdots.n\\}$), where $\textbf{A}\in\mathbb{R}^{d\times d}$ denotes
the coefficient matrix. This formulation frames the acyclicity in Equation (2)
as:

$[\textbf{W}(f)]_{ij}=|A_{ij}|.$ (3)

Nonlinear SEM. In the case of nonlinear SEM [10], we define
$f_{i}(\textbf{X})$ as a multilayer perceptron (MLP) with $h$ hidden layers
and an activation function $\sigma$ as:

$f_{i}(\textbf{X})=\text{MLP}(\textbf{X};\textbf{A}_{i}^{(1)},\cdots,\textbf{A}_{i}^{(h)})=\sigma(\cdots\sigma(\textbf{X}\textbf{A}_{i}^{(1)})\cdots)\textbf{A}_{i}^{(h)},$
(4)

where $\textbf{A}_{i}^{(l)}\in\mathbb{R}^{m_{l-1}\times m_{l}}$ is the
learnable weight matrix of the $l$-th hidden layer for the $i$-th node,
$m_{l}$ is the number of hidden units in the $l$-th layer, and $m_{0}=d$.
According to [10], we could define $\textbf{W}_{\theta}(f)$ as

$[\textbf{W}_{\theta}(f)]_{ij}=\left\lVert i\text{-th
column}(\textbf{A}_{j}^{(1)})\right\rVert_{2}.$ (5)

Then the acyclicity constraint could be latched on $\textbf{W}_{\theta}(f)$,
which serves as the replacement of the intractable $\textbf{W}(f)$ for
nonlinear setting.

##  3 Methodology

In this section, we first present differentiable invariant causal discovery
(DICD) to conduct causal structural learning over multiple environments. Then,
we give detailed formulations in both linear and nonlinear settings.

###  3.1 Differentiable Invariant Causal Discovery (DICD)

Despite the great success, the ERM paradigm easily captures spurious
correlations between variables by over-reconstructing the observations [18]
(See the toy example in Table I). It is essential to distinguish the causation
edges from the spuriously-correlated edges. Towards this end, we get access to
the multi-environment information and incorporate the idea of invariant
learning, so as to frame the causal discovery task as identifying environment-
invariant causation edges and discarding environment-dependent correlations.

First, we argue that the multi-environment, in various forms like explicit or
implicit metadata, is common in many real-world datasets and can partition the
data into different groups or domains. The environments in LFW dataset [27],
for instance, can divide the data into black-and-white and colorful
photographs. In ImageNet dataset [28], the data can be grouped by different
sources and years of images. Moreover, several benchmarks for the domain-
shifts problems have been provided by WILDS [29]. Among them, Camelyon17 [30]
includes the images collected from five hospitals that serve as different
environments, while Amazon [31] provides review texts from different reviewers
that can work as environments. Towards the end, the formal definition of the
different environments in our paper is given as follow:

######  Definition 3.1.

For two datasets generated from the same structure equation model as shown in
Equation (1). if there exists $i\in\\{1,\cdots,d\\}$ such that the
distribution of $z_{i}$ is different across these datasets and $X_{i}$ is not
the source node (see Definition 3.2) in the corresponding graph. Then we say
these two datasets are drawn from different environments.

######  Definition 3.2.

We define the node in the graph with no parents as the _Source Node_.

With the environment information, we utilize the invariant learning to conduct
multi-environment causal discovery — the function parameters of SEM from the
target DAG should remain optimal across all the environments. Guided by this
idea, our DICD consists of two modules: (1) Invariant structural model S,
which presents the structure of DAG, _i.e.,_ a binary adjacency matrix. By
“invariant”, we mean that S should be consistent across environments; and (2)
Optimal causal function $f$, which depicts the causal relation of variables.
By “optimal”, we mean that $f$ should be optimal over all the environments
once the structural model S is given. Considering the example in Table I
again, S is the DAG structure, while $f$ refers to the coefficients of edges.
For an invariant DAG structure, completely various optimal coefficients will
be learned in different environments. As a result, by utilizing the
environment information, incorrect DAGs with smaller reconstruction losses can
be excluded based on S and $f$.

Having environment-aware groups of observations, we define the empirical risk
within the environment $e\in\mathcal{E}$ as:

$\displaystyle\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ
f)=\frac{1}{n_{e}}\sum_{i=1}^{n_{e}}l(\textbf{X}^{e}_{i},(\textbf{S}\circ
f)(\textbf{X}^{e}_{i})),$

where $\circ$ refers to the composition of the two functions, and $n_{e}$ is
the number of samples in environment $e$. We then build two constraints on S
and $f$ across all environments, and establish the model of DICD:

$\displaystyle\min_{\textbf{S}\circ f}$
$\displaystyle~{}\sum_{e\in\mathcal{E}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ
f),$ (6) s.t.
$\displaystyle~{}~{}\mathcal{G}(\textbf{S})=\mathcal{G}(f)\in\text{DAG},$ (7)
$\displaystyle~{}~{}f=\arg\min_{\bar{f}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\bar{f}),\quad\forall
e\in\mathcal{E}.$ (8)

Equation (7) states the DAG constraint, where the DAGs represented by S and
$f$ are equivalent. In other words, $f$ latches on S’s DAG structure. Equation
(8) claims that $f$ refers to the optimal causal model fitting Equation (6)
across all environments. Following NOTEARS [10], the DAG constraints can
rewrite as hard DAG equation constraints:

$\displaystyle\mathcal{G}(\textbf{S})=\mathcal{G}(f),\quad
h(\textbf{W}(f))=0,$ (9)

where $\textbf{W}(f)$ is defined as Equation (2) in the linear setting and
should be replaced with $\textbf{W}_{\theta}(f)$ in Equation (5) for the
nonlinear setting. Besides
$h(\textbf{W})=\text{tr}(e^{\textbf{W}\circ\textbf{W}})-d$. Here, $\circ$
refers to the Hadamard product (_aka._ the element-wise product), and
tr($e^{\textbf{W}\circ\textbf{W}}$) is the trace of
$e^{\textbf{W}\circ\textbf{W}}$.

###  3.2 Linear SEM

When $f$ characterizes acyclicity in the linear SEM, we have
$f(\textbf{X})=\textbf{X}\textbf{A}$ and
$\partial_{i}f_{j}(\textbf{X})=A_{ij}$ (_cf._ Equation (3)). With the
invariant structure model S as the binary matrix, we can convert
$\textbf{S}\circ f$ as $\textbf{S}\circ\textbf{A}$, which is the Hadamard
product of S and A. However, the bilevel optimization formulation of the DICD
model is difficult to solve and is susceptible to failure due to over-
parametrization in our setting. We further simplify the learning of these two
matrices as the optimization of a new matrix
$\textbf{A}_{\textbf{S}}=\textbf{S}\circ\textbf{A}\in\mathbb{R}^{d\times d}$.
As such, we restrict the loss function in Equation (6) with the DAG constraint
in Equation (7) as:

$\displaystyle\min_{\textbf{A}_{\textbf{S}}}\sum_{e\in\mathcal{E}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{A}_{\textbf{S}}),\quad\text{s.t.}~{}h(\textbf{A}_{\textbf{S}})=0.$

This relaxed version only focuses on optimizing over $\textbf{A}_{\textbf{S}}$
during training. After learning $\textbf{A}_{\textbf{S}}$, we can simply reset
coefficient matrix $\textbf{A}=\textbf{A}_{\textbf{S}}$ and let S be the
binary indicator on $\textbf{A}_{\textbf{S}}$’s elements.

We then explore the tractable formulation for the optimality constraint across
environments in Equation (8) with the following theorem:

######  Theorem 3.1.

In the linear setting, we define a matrix variable as
$\textbf{B}\in\mathbb{R}^{d\times d}$. After replacing
$\textbf{S}\circ\textbf{A}$ with $\textbf{A}_{\textbf{S}}$, the constraint
Equation (7-8) satisfy the following necessary condition:

$\displaystyle\left\lVert\frac{\partial\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{A}_{\textbf{S}}\circ\textbf{B})}{\partial\textbf{B}}|_{\textbf{B}=\textbf{1}}\right\rVert^{2}_{2}=0,\quad\forall
e\in\mathcal{E},$ (10)

where 1 is the all-one matrix.

###### Proof.

The basic optimality constraint across environments for Equation (7-8) can be
described as:

$\displaystyle\textbf{A}^{*}=\arg\min_{\textbf{A}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}),\quad\forall
e\in\mathcal{E},$
$\displaystyle\text{s.t.}~{}\mathcal{G}(\textbf{S})=\mathcal{G}(\textbf{A})\in\text{DAG}.$
(11)

Then we introduce a new matrix variable B and insert it into Equation (3.2).
The optimal of B can be defined as:

$\displaystyle\textbf{B}^{*}=\arg\min_{\textbf{B}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*}\circ\textbf{B}),\quad\forall
e\in\mathcal{E}.$ (12)

We assert $\textbf{B}^{*}$ could be the all-one matrix 1, which indicates:

$\displaystyle\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*}\circ\textbf{1})\leq\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*}\circ\textbf{B}),\quad\forall\textbf{B}\in\mathbb{R}^{D\times
D},\forall e\in\mathcal{E}.$ (13)

Assuming there exists $\textbf{B}^{\prime}$ that satisfies
$\displaystyle\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*}\circ\textbf{B}^{\prime})<\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*}\circ\textbf{1})$,
then replacing $\textbf{A}^{*}\circ\textbf{B}^{\prime}$ with
$\textbf{A}^{\prime}$, we could have
$\displaystyle\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{\prime})<\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{S}\circ\textbf{A}^{*})$.
This result obviously contradicts with Equation (3.2). As such, Equation (13)
holds.

We replace $\textbf{S}\circ\textbf{A}$ with $\textbf{A}_{\textbf{S}}$. Now if
A is the optimal parameter of Equation (3.2) (_i.e.,_ $\textbf{A}^{*}$), then
$\textbf{A}_{\textbf{S}}$ becomes $\textbf{A}_{\textbf{S}}^{*}$, which yield
the following equation:

$\textbf{1}=\arg\min_{\textbf{B}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{A}_{\textbf{S}}^{*}\circ\textbf{B}),\quad\forall
e\in\mathcal{E}.$

According to the first-order optimality condition, we have:
$\displaystyle\frac{\partial\mathcal{L}^{e}(\textbf{A}_{S}^{*}\circ\textbf{A})}{\partial\textbf{B}}\big{|}_{\textbf{B}=\textbf{1}}=\textbf{0}$.
This concludes the proof. ∎

Based on Theorem 3.1, we could find that the intractable optimality constraint
across environments is now differentiable with the objective Equation (10).

We establish the DICD model in linear cases as follows:

$\displaystyle\min_{\textbf{A}_{\textbf{S}}}$
$\displaystyle\sum_{e\in\mathcal{E}}\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{A}_{\textbf{S}})+\lambda\sum_{e\in\mathcal{E}}\left\lVert\frac{\partial\mathbf{\mathop{\mathcal{L}}}^{e}(\textbf{A}_{\textbf{S}}\circ\textbf{B})}{\partial\textbf{B}}|_{\textbf{B}=\textbf{1}}\right\rVert^{2}_{2},$
s.t. $\displaystyle h(\textbf{A}_{\textbf{S}})=0.$ (14)

This objective function potentially incorporates the invariant structure and
the optimal coefficients into the single variable $\textbf{A}_{\textbf{S}}$,
thus only adds one more penalty for training compared with NOTEARS.

###  3.3 Nonlinear SEM

In nonlinear settings, We shall continue to concentrate on the scenario with
scalar-valued variables (_i.e.,_ $\textbf{X}\in\mathbb{R}^{n\times d}$) for
simplicity. However, the vector-valued variables (_i.e.,_
$\textbf{X}\in\mathbb{R}^{n\times d\times d_{x}}$) could be simply
incorporated into DICD.

NOTEARS-MLP [10] deploys the continuous DAG constraint for nonlinear SEM to
the first layer of all the MLPs, as shown in Equation (5). Inspired by this,
we establish the relationship between the structure matrix S and the first
layer of the MLPs. In other words, ensuring the parameters of MLP latch on the
structure S in order to satisfy $\mathcal{G}(\textbf{S})=\mathcal{G}(f)$
constraint. The precise solution to the equation
$\mathcal{G}(\textbf{S})=\mathcal{G}(f)$ is to maintain the following
restriction throughout the optimization process:

$\forall i,j,\quad\textbf{S}_{ji}=0\Longrightarrow||j\textrm{-th
column}(\textbf{A}_{i}^{(1)})||_{2}=0,$ (15)

where $\textbf{A}_{i}^{(1)}$ is the matrix parameter of the first layer in the
$i$-th MLP. The intuition behind the above equation is that setting
$j\textrm{-th column}(A_{i}^{(1)})$ to zero would block the correlation from
$j$ to $i$, corresponding to $\textbf{S}_{ji}=0$.

Though the first term in constraint Equation (9) could be formulated in
Equation (15). The optimization of Equation (6-8) is still intractable. We
need to construct an differentiable term to replace the constraint Equation
(15). To address this problem, we propose to simplify the optimization of S
and $f$ to optimizing a new function
$f_{\textbf{S}}=(f_{\textbf{S}1},\cdots,f_{\textbf{S}d})$, where
$f_{\textbf{S}i}=\textbf{S}\circ f_{i}$. The exact correlation is described in
the following equation:

$\displaystyle f_{\textbf{S}i}$ $\displaystyle=\textbf{S}\circ
f_{i}=\textrm{MLP}({Re([\textbf{S}]_{i},m_{1})}^{\top}\circ\textbf{A}_{i}^{(1)},\textbf{A}_{i}^{(2)},\cdots,\textbf{A}_{i}^{(h)})$
$\displaystyle=\textrm{MLP}(\textbf{A}_{\textbf{S}_{i}}^{(1)},\textbf{A}_{i}^{(2)},\cdots,\textbf{A}_{i}^{(h)}),$

where $m_{1}$ comes from $\textbf{A}_{i}^{(1)}\in\mathbb{R}^{m_{0}\times
m_{1}}$, and ${Re([\textbf{S}]_{i},m_{1})}^{\top}$ refers to
${[\underbrace{[\textbf{S}]_{i},\cdots,[\textbf{S}]_{i}}_{m_{1}\textrm{
times}}]}^{\top}$, with $[\textbf{S}]_{i}$ being the $i$-th column of S. We
denote $Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)}$ as
$\textbf{A}_{\textbf{S}i}^{(1)}$. Similar to the linear setting, once
$f_{\textbf{S}}$ is learned, we can simply set $f=f_{\textbf{S}}$ and S as the
adjacency matrix of $\mathcal{G}(f)$. Then $f$ and S will share the same graph
structure. Now we propose the following theorem for nonlinear setting:

######  Theorem 3.2.

In nonlinear system, we denote $\textbf{A}_{i}^{(l)},l\in\\{1,\cdots,d\\}$ as
the parameter of the $l$-th layer of MLP $f_{i}$, and
$\textbf{A}_{\textbf{S}i}$ is the parameter of the 1st layer of MLP
$f_{\textbf{S}i}$, $i\in\\{1,\cdots,d\\}$. Then after replacing
$Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)}$ with
$\textbf{A}_{\textbf{S}i}^{(1)}$, we define a matrix variable as
$\textbf{B}\in\mathbb{R}^{m_{1}\times d}$, the optimal condition Equations
(7-8) except for the DAG part satisfy the following necessary condition:

$\displaystyle\sum_{i=1}^{d}\Big{|}\Big{|}$
$\displaystyle\frac{\partial\mathcal{L}^{e}(\textrm{MLP}(\textbf{A}_{\textbf{S}i}^{(1)}\circ\textbf{B},\textbf{A}_{i}^{(2)},\cdots,\textbf{A}_{i}^{(n)}))}{\partial\textbf{B}}\Big{|}_{\textbf{B}=\textbf{1}}\Big{|}\Big{|}_{2}^{2}=0,\quad$
$\displaystyle\forall e\in\mathcal{E}.$ (16)

###### Proof.

The basic optimal condition described in Equation (7-8) without the DAG
constraint could be written as:

$\displaystyle\textbf{A}_{i}^{(1)*}$
$\displaystyle=\textrm{arg}\min_{A_{i}^{(1)}}\mathcal{L}^{e}(\textrm{MLP}(Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)},\cdots,\textbf{A}_{i}^{(n)})),$
$\displaystyle\forall e$
$\displaystyle\in\mathcal{E},\quad\textrm{s.t.}\,\,\mathcal{G}(\textbf{S})=\mathcal{G}(f),$
(17)

where $Re([\textbf{S}_{i}],m_{1})$ means
${[\underbrace{[\textbf{S}]_{i},\cdots,[\textbf{S}]_{i}}_{m_{1}\textrm{
times}}]}^{\top}$, with $[\textbf{S}]_{i}$ being the $i$-th column of S. Note
that the following equations with $e$ in this proof indicate all the
environments, we omit the condition $\forall e\in\mathcal{E}$ in the following
to simplify the proof.

Similar to the proof of Proposition 3.1, we define a new matrix
$\textbf{B}\in\mathbb{R}^{m_{1}\times d}$, and insert this term into the above
equation, then we can have:

$\textbf{B}^{*}=\textrm{arg}\min_{B}\mathcal{L}^{e}(\textrm{MLP}(Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)*}\circ\textbf{B},\cdots,\textbf{A}_{i}^{(n)})).$
(18)

We argue that $\textbf{B}^{*}$ could be 1, a matrix of dimensions $m_{1}\times
d$ filled with ones. This is expressed as:

$\textbf{1}=\textrm{arg}\min_{B}\mathcal{L}^{e}(\textrm{MLP}(Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)*}\circ\textbf{B},\cdots,\textbf{A}_{i}^{(n)})).$

Then replacing $Re([\textbf{S}]_{i},m_{1})\circ\textbf{A}_{i}^{(1)*}$ with
$\textbf{A}_{\textbf{S}i}^{*}$ would yield:

$\textbf{1}=\textrm{arg}\min_{B}\mathcal{L}^{e}(\textrm{MLP}(\textbf{A}_{\textbf{S}i}^{*}\circ\textbf{B},\textbf{A}_{i}^{(2)},\cdots,\textbf{A}_{i}^{(n)})).$

Again, by the first-order optimality condition, we have Equation (16). ∎

With Theorem 3.2, we could transform the intractable constraint in Equation
(8) into the differentiable term subject to the first layer of MLP. The
detailed formulation in the nonlinear setting can be expressed as:

$\displaystyle\min_{f_{\textbf{S}}}\sum_{e\in\mathcal{E}}L^{e}(f_{\textbf{S}})$
$\displaystyle+\lambda\sum_{e\in\mathcal{E}}\sum_{i=1}^{d}\Big{|}\Big{|}\frac{\partial\mathcal{L}^{e}(\textrm{MLP}(\textbf{A}_{\textbf{S}i}\circ\textbf{B},\textbf{A}_{i}^{(2)},\cdots,\textbf{A}_{i}^{(n)}))}{\partial\textbf{B}}\Big{|}_{\textbf{B}=\textbf{1}}\Big{|}\Big{|}_{2}^{2},$
$\displaystyle\textrm{ s.t. }h(\textbf{W}_{\theta}(f_{\textbf{S}}))=0.$ (19)

This objective function also potentially uses $f_{\textbf{S}}$ to both
incorporate $f$ and the invariant structure matrix S.

##  4 Theoretical Analysis

In this section, we aim to provide the sufficient conditions for
identifiability of DICD in the linear SEM systems. We will leave the
discussion about which assumptions can or cannot be further relaxed in future
work.

######  Theorem 4.1.

For linear SEMs systems with Gaussian additive noises as in Equation (1), if
for any $X_{i}\in\mathcal{V}$ that is not the source node, there exist two
environments $e_{1},e_{2}\in\mathcal{E}$, such that:

$\displaystyle Var(z_{i}^{e_{1}})\neq Var(z_{i}^{e_{2}}),$ (20)
$\displaystyle\forall
X_{j}\in\mathcal{V}\backslash\\{X_{i}\\},Var(z_{j}^{e_{1}})=Var(z_{j}^{e_{2}}),$
(21)

where $Var(\cdot)$ represents the variance, $z_{j}^{e_{1}}$ and
$z_{j}^{e_{2}}$ are the additive noise from $Pa(X_{j})$ to $X_{j}$ in $e_{1}$
and $e_{2}$, respectively. Then the causal structure is identifiable.

Here we emphasize that the graph and the coefficients that satisfy the
following sufficient conditions and Equation (6-8) will be exactly the true
causal graph and the true coefficients, corresponding to the identifiable
graph.

Theorem 4.1 indicates that our DICD is guaranteed to retrieve the true causal
graph in linear systems, when the diversity of environments is adequate.

The assumptions in Theorem 4.1 are restrictive mainly because Equation (21))
requires that the variances of the nodes other than $X_{i}$ to be the same
across environments $e_{1}$ and $e_{2}$. In our proof, Equation (21) is only
used for proving Lemma 4.2. However, we will show that, even if we don’t
assume Equation (21), we may still obtain Equation 32 in the proof of Lemma
4.2 below, which leads to the proof of Lemma 4.2.

With Equation (31), we know the sufficient and necessary condition of
$\hat{w}_{j_{0}}^{e_{1}}=\hat{w}_{j_{0}}^{e_{2}}$ is:

$\displaystyle\frac{\mathbb{E}_{e_{1}}[X_{i}X_{j_{0}}]-\sum_{k\in
Pa_{s}(i)\backslash\\{j_{0}\\}}\hat{w}^{e}_{k}\mathbb{E}_{e_{1}}[X_{k}X_{j_{0}}]}{\mathbb{E}_{e_{1}}[X_{j_{0}}^{2}]}$
$\displaystyle=\frac{\mathbb{E}_{e_{2}}[X_{i}X_{j_{0}}]-\sum_{k\in
Pa_{s}(i)\backslash\\{j_{0}\\}}\hat{w}^{e}_{k}\mathbb{E}_{e_{2}}[X_{k}X_{j_{0}}]}{\mathbb{E}_{e_{2}}[X_{j_{0}}^{2}]}.$

Since we do not have Equation (21)), every single term in the above equation
is not necessarily equivalent. The probability of the combinations of all
these non-equivalent terms being equivalent is very small. Thus during the
implementation, Lemma 4.2 is easy to be true. Then the other parts of the
proof of Theorem 4.1 will remain the same.

Next, we provide the proof skeleton of Theorem 4.1 and leave the full proof in
Appendix Proof of Theorem 4.1.

There are three main steps in our proof. First, we prove that if the causal
structure is correct, then the optimal coefficients are the ground truth
coefficients (Lemma 4.1). Second, we show that the stable graph (see
definition 4.1) with the minimal sum of reconstruction loss from all
environments is exactly the true graph (Theorem 4.2). In the end, we can prove
that our algorithm could yield the ground truth graph and coefficients with
adequate environments.

To begin with, we give the following lemma to show that our DICD could yield
the true coefficient in the linear SEM under the true causal structure.

######  Lemma 4.1.

Given the true graph $G_{0}$, the corresponding structure $\textbf{S}_{0}$ is
the adjacency matrix of $G_{0}$. $\forall e\in\mathcal{E}$, we denote the
optimal parameters for the coefficients in environment $e$ as:

$\hat{\textbf{W}}^{e}=arg\min_{\overline{\textbf{W}}}\mathcal{L}^{e}(\textbf{S}_{0}\circ\overline{\textbf{W}}),\quad\emph{s.t.
}\mathcal{G}(\textbf{S}_{0})=\mathcal{G}(\overline{\textbf{W}}).$ (22)

Then we have $\hat{\textbf{W}}^{e}=\textbf{W}_{0},\forall e\in\mathcal{E}$,
where $\textbf{W}_{0}$ is the ground truth coefficients.

The following definitions are necessary for the rest of this section.

######  Definition 4.1.

For a given graph $G$, if there exist $\textbf{S},\textbf{W}$ such that:

$\displaystyle\mathcal{G}(\textbf{S})=\mathcal{G}(\textbf{W})=G,$
$\displaystyle\textbf{W}=arg\min_{\overline{\textbf{W}}}\mathcal{L}^{e}(\textbf{S}\circ\overline{\textbf{W}}),\quad\forall
e\in\mathcal{E},$ (23)

then we call $G$ a stable graph.

######  Definition 4.2.

For two nodes $X_{1}$ and $X_{2}$ in a DAG, if $X_{2}$ is reachable from
$X_{1}$, then $X_{1}$ is a predecessor of $X_{2}$. We denote all the
predecessors of $X_{2}$ in graph $G_{0}$ as $Pre_{0}(X_{2})$.

With Definition 4.1, we further denote the parents of the variable $X_{i}$ in
$G_{0}$ as $Pa_{0}(X_{i})$ and the corresponding indexes as $Pa_{0}(i)$. We
propose the following Lemma 4.2 to demonstrate that the causal directions
between any two variables cannot violate each other in any stable graph and
the true causal graph, which could serve as the pre-conditions in Theorem 4.2
to prove it.

######  Lemma 4.2.

For any given stable graph $G_{s}$, if we assume the conditions in Theorem 4.1
hold, then $\forall X_{i}\in\mathcal{V}$ which is not a source node, we have
$X_{i}\notin Pre_{0}(X_{j})$ for any $X_{j}\in Pa_{s}(X_{i})$.

Finally, Theorem 4.1 can be implied by the following theorem:

######  Theorem 4.2.

For any given stable graph $G_{s}$ and the ground truth graph $G_{0}$, we
denote their corresponding structures as $\textbf{S}_{0}$ and
$\textbf{S}_{s}$, respectively. We further denote their corresponding
consistent optimal parameters as $\textbf{W}_{s}$ and $\textbf{W}_{0}$,
respectively. Then we have:

$\sum_{e\in\mathcal{E}}\mathcal{L}^{e}(\textbf{S}_{0}\circ\textbf{W}_{0})\leq\sum_{e\in\mathcal{E}}\mathcal{L}^{e}(\textbf{S}_{s}\circ\textbf{W}_{s}),$
(24)

and the equation holds only for $\textbf{W}_{0}=\textbf{W}_{s}$.

The proof of Theorem 4.2 is given in Appendix Proof of Theorem 4.1. Theorem
4.2 imples that, with the conditions in Equation (7) and Equation (8), The
ground truth graph will be yielded by Equation (6).

##  5 Experiments

In this section, we study the empirical performance of our proposed method. We
aim to answer the following research questions:

  * •

RQ1: How does DICD perform compared to the previous methods in both linear and
nonlinear settings?

  * •

RQ2: How do DICD and other baselines perform with various factors (_i.e.,_ the
number of environments, density of graph).

  * •

RQ3: How does DICD perform on real-world datasets compared with other
applicable baselines?

###  5.1 Experimental Settings

We now provide the detailed settings for our experiments. The descriptions for
the synthetic datasets and the real-world dataset are in Section 5.1.4 and
Section 5.1.5, respectively. For all the experiments in this paper, we all
generate 10 datasets for each graph setting and report the mean and standard
deviation.

####  5.1.1 Baselines

We select four state-of-the-art causal discovery methods for comparison:

  * •

CD-NOD [21] is a constrained-based causal discovery method designed for
heterogeneous datasets, _i.e.,_ datasets from different environments. CD-NOD
utilizes the independent changes across environments to determine the causal
orientations, and proposes constrained-based and kernel-based methods to find
the causal structure.

  * •

NOTEARS [1] is specifically designed for linear setting, and is also the
backbone of DICD in linear cases. NOTEARS estimates the true causal graph by
minimizing the reconstruction loss with the continuous acyclicity constraint.
We re-implement NOTEARS with replacing the L-BFGS-B iteration with Adam
gradient descent, which could yield compatible performance and more
importantly, could be deployed on GPU.

  * •

NOTEARS-MLP [10] is specifically designed for nonlinear setting, which also
serves as the foundation of DICD in nonlinear situations. NOTEARS-MLP
approximates the generative SEM model by MLP while only constraining the first
layer of the MLP with the continuous acyclicity constraint.

  * •

DAGGNN[7] formulates causal discovery with variational autoencoder, where the
encoder and decoder are all graph neural networks. Choosing the evidence lower
bound as the loss function and slightly modifying the acyclicity constraint,
DAGGNN could manage to recover the weighted adjacency matrix.

  * •

NOCURL [12] utilizes a two-step procedure: first find an initial cyclic
solution, then employ Hodge decomposition of graphs and learn an acyclic graph
by projecting the cyclic graph to the gradient of a potential function.

  * •

DARING [18] imposes explicit residual independence constraint with an
adversarial strategy. We choose the backbone as NOTEARS-MLP to conform with
the settings above.

####  5.1.2 Hyperparameter Settings

For linear settings, there are two hyper-parameters in total: $\lambda_{1}$
for the $l_{1}$-norm regularization term; $\lambda_{D}$ for the DICD penalty
term. We tune $\lambda_{1}$ in $\\{0.01,0.1\\}$ for NOTEARS and DICD. Besides,
we tune $\lambda_{D}$ in $\\{0.1,1\\}$ for DICD. Then for nonlinear settings,
there are three hyper-parameters in total:
$\lambda_{1},\lambda_{2},\lambda_{D}$, among which $\lambda_{1}$ and
$\lambda_{2}$ are for the $l_{1}$-norm and $l_{2}$-norm regularization terms,
respectively. We tune $\lambda_{1},\lambda_{2}$ both in $\\{0.01,0.1\\}$ and
$\lambda_{D}$ in $\\{0.1,1\\}$. The scheduler for $\lambda$ in Equation (14)
and Equation (19) is shown as follows:

$\lambda=\left\\{\begin{array}[]{cc}\displaystyle\frac{k}{K/3}\circ\lambda_{D}&k\leq
K/3\\\ \displaystyle\lambda_{D}&K/3\leq k\leq 2K/3\\\
\displaystyle\frac{K-k}{K/3}\circ\lambda_{D}&k\geq 2K/3\end{array}\right.,$

where $K$ is the estimated total step, and $k$ is the current iteration. The
intuition behind this scheduler is that we need to let the model fit the data
at the beginning, then as the training process goes, we need to enforce our
penalty to help find the true causal graph. Then for the last stage, the graph
structure has almost been inferred. We need to gradually remove our penalty to
let the structural causal function with the given structure fit the data.

TABLE II: Linear Setting, for ER and SF graphs of 10, 20, 50 nodes

| 10 nodes | 20 nodes | 50 nodes
---|---|---|---
ER4 | FDR | TPR | SHD | FDR | TPR | SHD | FDR | TPR | SHD
CD-NOD | 0.48$\pm$0.06 | 0.17$\pm$0.02 | 33.5$\pm$0.9 | 0.48$\pm$0.19 | 0.11$\pm$0.02 | 75.7$\pm$5.0 | 0.56$\pm$0.08 | 0.15$\pm$0.06 | 195.0$\pm$8.0
NOTEARS | 0.07$\pm$0.01 | 0.78$\pm$0.02 | 10.1$\pm$0.8 | 0.20$\pm$0.03 | 0.71$\pm$0.08 | 35.9$\pm$7.4 | 0.27$\pm$0.04 | 0.78$\pm$0.03 | 96.7$\pm$13.7
DAGGNN | 0.08$\pm$0.01 | 0.86$\pm$0.02 | 7.6$\pm$0.9 | 0.34$\pm$0.04 | 0.79$\pm$0.05 | 48.3$\pm$7.1 | 0.36$\pm$0.03 | 0.86$\pm$0.03 | 122.3$\pm$14.6
NOCURL | 0.14$\pm$0.01 | 0.81$\pm$0.01 | 9.8$\pm$0.4 | 0.19$\pm$0.02 | 0.92$\pm$0.01 | 22.3$\pm$2.4 | 0.31$\pm$0.02 | 0.94$\pm$0.01 | 89.5$\pm$9.5
DARING | 0.08$\pm$0.01 | 0.81$\pm$0.03 | 9.3$\pm$1.3 | 0.38$\pm$0.07 | 0.58$\pm$0.07 | 60.9$\pm$9.9 | 0.48$\pm$0.01 | 0.60$\pm$0.06 | 187.6$\pm$4.9
DICD | 0.03$\pm$0.01 | 0.88$\pm$0.02 | 4.9$\pm$1.0 | 0.16$\pm$0.03 | 0.89$\pm$0.04 | 19.7$\pm$5.8 | 0.26$\pm$0.02 | 0.89$\pm$0.06 | 82.0$\pm$5.3
SF4 | FDR | TPR | SHD | FDR | TPR | SHD | FDR | TPR | SHD
CD-NOD | 0.36$\pm$0.10 | 0.19$\pm$0.02 | 25.0$\pm$1.4 | 0.34$\pm$0.05 | 0.18$\pm$0.01 | 59.3$\pm$0.9 | 0.38$\pm$0.04 | 0.15$\pm$0.01 | 168.3$\pm$1.7
NOTEARS | 0.04$\pm$0.04 | 0.81$\pm$0.02 | 6.1$\pm$1.1 | 0.19$\pm$0.01 | 0.77$\pm$0.02 | 27.1$\pm$1.5 | 0.17$\pm$0.01 | 0.83$\pm$0.01 | 60.7$\pm$2.5
DAGGNN | 0.07$\pm$0.03 | 0.98$\pm$0.02 | 2.9$\pm$1.5 | 0.27$\pm$0.03 | 0.84$\pm$0.02 | 31.6$\pm$2.6 | 0.26$\pm$0.02 | 0.88$\pm$0.01 | 80.6$\pm$8.3
NOCURL | 0.06$\pm$0.02 | 0.86$\pm$0.02 | 4.8$\pm$1.1 | 0.25$\pm$0.01 | 0.86$\pm$0.01 | 28.2$\pm$1.3 | 0.26$\pm$0.08 | 0.93$\pm$0.05 | 71.8$\pm$10.6
DARING | 0.17$\pm$0.04 | 0.86$\pm$0.07 | 9.0$\pm$2.7 | 0.26$\pm$0.02 | 0.80$\pm$0.01 | 32.6$\pm$1.8 | 0.28$\pm$0.02 | 0.87$\pm$0.01 | 87.3$\pm$5.9
DICD | 0.05$\pm$0.04 | 0.98$\pm$0.02 | 2.3$\pm$1.7 | 0.16$\pm$0.05 | 0.81$\pm$0.09 | 22.1$\pm$8.4 | 0.18$\pm$0.03 | 0.91$\pm$0.01 | 53.7$\pm$8.7

####  5.1.3 Evaluation Protocols

We use the three most popular metrics in causal discover: false discovery rate
(FDR), true positive rate (TPR) and structural Hamming distance (SHD). Higher
TPR stand for better performances, while FDR and SHD should be lower to
represent the better strategies.

####  5.1.4 Synthetic Datasets

We conduct experiments on two synthetic datasets for linear and nonlinear
settings. Besides, we apply our method DICD on Colored MNIST [32] dataset to
explore its effectiveness on real-world datasets. As for the synthetic data,
the ground truth DAG is generated from two random graph models: Erdos-Renyi
(ER) and scale-free (SF), following [10]. For the overall experimental
comparison, we set the node degree as four. For the linear and nonlinear
setting, we construct the environment variable $E$ to simulate the effects of
the environments on the additive noises. For the linear setting, after
generating the graph, we randomly select $\lfloor 0.3*d\rfloor$ nodes in this
graph, and then build $\lfloor 0.3*d\rfloor$ new nodes. We call these new
variables as the environment variables and denote them as $E$, and they
satisfy the same distributions and are all independent from each other. Then
we create $\lfloor 0.3*d\rfloor$ edges from each environment node to the each
selected node. In this way, we could simulate different environments with
varying the distribution of the environment variables $E$.

Then given this new graph with $d+\lfloor 0.3*d\rfloor$nodes, we simulate
random edge weights to obtain a new matrix
$\textbf{W}\in\mathbb{R}^{(d+\lfloor 0.3*d\rfloor)\times(d+\lfloor
0.3*d\rfloor)}$. With W, we sample
$\textbf{X}=\textbf{W}^{T}\textbf{X}+z\in\mathbb{R}^{d+\lfloor 0.3*d\rfloor}$
with $z$ from Gaussian noise model to generate 10 random datasets
$\textbf{X}_{e}\in\mathbb{R}^{n\times(d+\lfloor 0.3*d\rfloor)}$. Then we
remove the column corresponding to the additional $\lfloor 0.3*d\rfloor$
environment variables $E$ to generate the final datasets
$\textbf{X}\in\mathbb{R}^{n\times d}$. We change the variances of the noises
of the environment variables $E$ to simulate different environments. In the
nonlinear setting, after generating the graph, we randomly select $\lfloor
0.5*d\rfloor$ nodes in the graph, and then create $\lfloor 0.5*d\rfloor$
environment nodes and also $\lfloor 0.5*d\rfloor$ edges from each environment
node to the selected node. Then given this new graph, we simulate the SEM
$\textbf{X}_{j}=F_{j}(\textbf{X}_{pa(j)})+z_{j}$ for all
$j\in\\{1\cdots,d+\lfloor 0.5*d\rfloor\\}$ in topological order. For the
environment variables, we vary the distributions of the noises in different
environments. Then for the other nodes, we set the distribution of the noises
as $\mathcal{N}(0,1)$. We choose $f_{j}$ to be Additive Noise Models with two-
layer MLPs. We will also remove the column corresponding to the variable $E$
to generate the datasets $\textbf{X}\in\mathbb{R}^{n\times d}$.

####  5.1.5 Real-world Dataset

For the real-world dataset, We sample 10000 images in total, and 2000 images
for each environment. We classify MNIST digits from 2 classes, where classes 0
and 1 indicate original digits (0,1,2,3,4) and (5,6,7,8,9). Then for each
environment, we select the ratio of class 0 being green, which is shown in the
”Ratio” column in Table III. Then the noise variances(i.e. noise scale) for
each environment are provided in the “Noise Scale” column in Table III.

TABLE III: Experimental Settings for Colored MNIST. Environment | Ratio | Noise Scale
---|---|---
$e_{1}$ | 0.16 | 10/255
$e_{2}$ | 0.32 | 20/255
$e_{3}$ | 0.48 | 30/255
$e_{4}$ | 0.64 | 40/255
$e_{5}$ | 0.80 | 50/255
TABLE IV: Nonlinear Setting, for ER and SF graphs of 10, 20, 50 nodes

| 10 nodes | 20 nodes | 50 nodes
---|---|---|---
ER4 | FDR | TPR | SHD | FDR | TPR | SHD | FDR | TPR | SHD
CD-NOD | 0.39$\pm$0.06 | 0.50$\pm$0.07 | 20.0$\pm$3.1 | 0.31$\pm$0.04 | 0.56$\pm$0.06 | 56.8$\pm$4.2 | 0.35$\pm$0.07 | 0.82$\pm$0.05 | 115.7$\pm$18.3
NOTEARS-MLP | 0.24$\pm$0.10 | 0.44$\pm$0.12 | 23.8$\pm$3.5 | 0.25$\pm$0.05 | 0.35$\pm$0.09 | 59.0$\pm$5.0 | 0.30$\pm$0.08 | 0.86$\pm$0.06 | 102.9$\pm$25.4
DAGGNN | 0.50$\pm$0.06 | 0.21$\pm$0.04 | 32.2$\pm$1.5 | 0.61$\pm$0.07 | 0.24$\pm$0.05 | 82.9$\pm$4.5 | 0.62$\pm$0.06 | 0.15$\pm$0.03 | 217.0$\pm$16.0
NOCURL | 0.38$\pm$0.02 | 0.38$\pm$0.05 | 27.4$\pm$1.5 | 0.56$\pm$0.07 | 0.34$\pm$0.07 | 80.2$\pm$7.6 | 0.69$\pm$0.07 | 0.28$\pm$0.07 | 258.2$\pm$25.2
DARING | 0.44$\pm$0.04 | 0.28$\pm$0.04 | 29.8$\pm$2.7 | 0.55$\pm$0.09 | 0.23$\pm$0.07 | 77.6$\pm$5.7 | 0.58$\pm$0.09 | 0.23$\pm$0.05 | 209.4$\pm$20.0
DICD | 0.20$\pm$0.08 | 0.57$\pm$0.10 | 19.2$\pm$4.3 | 0.26$\pm$0.06 | 0.69$\pm$0.07 | 40.3$\pm$6.1 | 0.26$\pm$0.04 | 0.88$\pm$0.03 | 84.4$\pm$11.7
SF4 | FDR | TPR | SHD | FDR | TPR | SHD | FDR | TPR | SHD
CD-NOD | 0.42$\pm$0.03 | 0.68$\pm$0.04 | 19.5$\pm$2.5 | 0.31$\pm$0.09 | 0.64$\pm$0.07 | 46.3$\pm$6.8 | 0.32$\pm$0.05 | 0.75$\pm$0.09 | 104.3$\pm$12.1
NOTEARS-MLP | 0.41$\pm$0.10 | 0.31$\pm$0.16 | 25.2$\pm$4.6 | 0.29$\pm$0.10 | 0.56$\pm$0.10 | 51.0$\pm$7.4 | 0.29$\pm$0.08 | 0.72$\pm$0.13 | 113.7$\pm$17.3
DAGGNN | 0.62$\pm$0.09 | 0.26$\pm$0.07 | 31.2$\pm$3.4 | 0.69$\pm$0.07 | 0.16$\pm$0.03 | 80.1$\pm$7.4 | 0.64$\pm$0.03 | 0.15$\pm$0.05 | 210.3$\pm$8.2
NOCURL | 0.54$\pm$0.06 | 0.40$\pm$0.10 | 26.2$\pm$2.0 | 0.70$\pm$0.04 | 0.27$\pm$0.05 | 87.2$\pm$3.0 | 0.70$\pm$0.02 | 0.22$\pm$0.01 | 240.8$\pm$9.3
DARING | 0.54$\pm$0.10 | 0.27$\pm$0.05 | 29.0$\pm$2.7 | 0.58$\pm$0.06 | 0.21$\pm$0.03 | 73.1$\pm$4.5 | 0.53$\pm$0.03 | 0.19$\pm$0.02 | 193.0$\pm$5.0
DICD | 0.34$\pm$0.06 | 0.71$\pm$0.15 | 16.2$\pm$4.3 | 0.27$\pm$0.06 | 0.68$\pm$0.15 | 37.9$\pm$7.1 | 0.29$\pm$0.04 | 0.80$\pm$0.03 | 99.3$\pm$9.0

###  5.2 Overall Performances (RQ1)

We present the overall performances of DICD and the baselines for fair
comparison. In the baselines, NOTEARS, DAGGNN, NoCurl, DARING are run on the
concatenated datasets from all the environments. CD-NOD is run with the
environment-id corresponding to each sample.

####  5.2.1 Linear Synthetic Data

In this experiment, we explore the improvements when introducing different
groups by comparing the DAG estimations against the ground truth structure. We
simulate {ER4, SF4} graphs with $d=\\{10,20,50,100\\}$ nodes. For each
environment, we generate $200$ samples. We evaluate our methods with datasets
from 5 environments, and the variances of Gaussian noise for the environment
variable $E$ in each environment are {0.2, 0.4, 0.6, 0.8, 1.0}. Table II and
Table VII summarizes the results when the number of nodes equals to
$\\{10,20,50,100\\}$. From these tables, we could have the following key
observations: (1) DICD has outperformed all other baselines across various
settings. More precisely, DICD achieves significant improvements over the
strongest baselines by up to 36% in SHD (10 nodes, ER graph). (2) Generally,
DICD has the lower FDR and higher TPR, which also coincides with the intuition
that DICD could eliminate spurious correlations and reveal the true ones.

####  5.2.2 Nonlinear Synthetic Data

We also conducte enormous experiments to demonstrate the effectiveness of DICD
in the nonlinear setting. Similarly to the linear setting, we simulate {ER4,
SF4} graphs with $d=\\{10,20,50,100\\}$ nodes. For each graph, we generate
data from two environments, 1000 samples for each environment. The Gaussian
noises for the environment variable $E$ in these two groups are
$\\{0.2,0.4\\}$. The results with $d=\\{10,20,50\\}$ are provided in Table IV
and results with $d=100$ are shown in Table VII. From these tables, we could
find: (1) DICD consistently outperforms other baselines in all eight settings
upon the most crucial metric SHD. The improvements on SHD over the best
baseline are up to 29% (50 nodes, ER graph). (2) DICD achieves compatible FDR
with NOTEARS but far higher TPR. This shows that in the nonlinear setting,
DICD is better at revealing the true causal correlations that might have been
missed by NOTEARS. (3) The over-reconstruction problem still exists in other
methods, while DICD has the potential to mitigate it, which could be the
reason for the performance improvements. (4) CD-NOD performs fairly well in
nonlinear cases, which means the multi-environment setting might be more
helpful when the correlations between variables are more complicated. However,
CD-NOD consumes in average more than 9 hours in the simplest setting (10
nodes), and more than 300 hours in the case of 50 nodes, which is far more
expensive than our algorithm. As shown in Table V, we only record the running
time of CD-NOD in nonlinear settings, since its running time is almost
unacceptable in these cases. Then we report the averaged running time for
different seeds and different graph types (ER or SF). From the table, we can
observe that CD-NOD is very expensive for nonlinear cases, while the running
time of DICD almost remain constant when the number of nodes get larger.

TABLE V: Running time comparison | 10 nodes | 20 nodes | 50 nodes
---|---|---|---
CD-NOD |  $>$9h |  $>$64h |  $>$300h
DICD | 15min | 15min | 15min
TABLE VI: Linear Setting, for ER and SF graphs of 100 nodes

| ER4 | SF4
---|---|---
Method | FDR | TPR | SHD | FDR | TPR | SHD
NOTEARS | 0.17$\pm$0.03 | 0.78$\pm$0.03 | 147.9$\pm$15.7 | 0.12$\pm$0.04 | 0.91$\pm$0.02 | 81.4$\pm$23.0
DAGGNN | 0.23$\pm$0.03 | 0.79$\pm$0.04 | 178.3$\pm$9.6 | 0.31$\pm$0.02 | 0.92$\pm$0.01 | 195.3$\pm$20.4
NOCURL | 0.25$\pm$0.01 | 0.94$\pm$0.00 | 136.0$\pm$5.3 | 0.14$\pm$0.04 | 0.98$\pm$0.01 | 69.8$\pm$23.2
DARING | 0.29$\pm$0.02 | 0.69$\pm$0.04 | 234.0$\pm$14.3 | 0.29$\pm$0.01 | 0.90$\pm$0.01 | 180.9$\pm$5.4
DICD | 0.20$\pm$0.03 | 0.87$\pm$0.01 | 133.9$\pm$16.2 | 0.12$\pm$0.03 | 0.97$\pm$0.01 | 61.9$\pm$15.6

TABLE VII: Nonlinear Setting, for ER and SF graphs of 100 nodes

| ER4 | SF4
---|---|---
Method | FDR | TPR | SHD | FDR | TPR | SHD
NOTEARS | 0.20$\pm$0.04 | 0.40$\pm$0.07 | 270.0$\pm$24.1 | 0.17$\pm$0.06 | 0.42$\pm$0.10 | 260.1$\pm$27.3
DAGGNN | 0.52$\pm$0.05 | 0.09$\pm$0.01 | 390.6$\pm$7.7 | 0.50$\pm$0.08 | 0.09$\pm$0.02 | 390.3$\pm$13.4
DARING | 0.41$\pm$0.05 | 0.15$\pm$0.02 | 367.0$\pm$8.8 | 0.42$\pm$0.04 | 0.15$\pm$0.02 | 367.3$\pm$6.9
NOCURL | 0.62$\pm$0.06 | 0.23$\pm$0.02 | 447.2$\pm$33.0 | 0.60$\pm$0.03 | 0.18$\pm$0.02 | 418.6$\pm$10.7
DICD | 0.18$\pm$0.05 | 0.54$\pm$0.10 | 226.6$\pm$23.1 | 0.18$\pm$0.03 | 0.53$\pm$0.06 | 228.0$\pm$16.3

###  5.3 Study of Various Factors (RQ2)

In this section, we discuss various factors that may affect the performances
of DICD and other methods. Due to limit of space, we discuss the effect of
environmental imbalance in Appendix 5.3.2.

####  5.3.1 The Effect of Environment Number

To explore how the number of environments affects the performances, we conduct
experiments on both linear and nonlinear settings. The total number of the
examples is set to 1000. We choose $d=10$ and $s_{0}=40$ with ER graph for
this case study. If the number of the environments is $N$, then the number of
samples for each environment is $\lfloor 1000/N\rfloor$. The noises for
different environments are sliced from the head of the array {0.2, 0.4, 0.6,
0.8, 1.0, 0.1, 0.3, 0.5, 0.7, 0.9}. For example, the noises are {0.2, 0.4}
when the number of environments is 2. In Figure 1, we can observe:

(a) Linear Setting

(b) Nonlinear Setting

Figure 1: SHD _w.r.t._ the number of environments.

(1) When the number of environments is greater than 4 in the linear setup,
more environments will not yield to better results, which means for simple
structure equation models, four environments may have already ruled out the
possibility of incorrect graphs. (2) In the nonlinear setting, it is amazing
to find that DICD could achieve better performance even there is only one
environment. The reason could be that when there are more parameters in the
function, the explicit constraint on the optimality of the parameters could
contribute much to learning the best graph, which is more likely to be the
ground truth graph. (3) More environments in nonlinear settings could degrade
the performance of all approaches; the reason for this could be that
heterogeneous noises can be particularly unfriendly when the relationships
between variables become quite intricate.

####  5.3.2 The Effect of the Imbalance between Different Environments

Since the imbalance of data is the major problem in machine learning, we aim
to explore how the imbalance of data size between different environments would
affect the performances of our method and other baselines. The total number of
samples in this dataset is 1000 for linear settings and 2000 for nonlinear
settings. We choose ER graphs with $d=20$ and $s_{0}=4d=80$ for this case
study. The noises for the variable $E$ in two environments are set to be {0.2,
0.4}. Then the ratio in Figure 2 means the percentage of the samples from the
first environment. From this figure, we could have the following observations:
(1) In most of the settings, DICD outperforms other methods consistently,
except when the data is higher imbalanced in linear setting, where there is
almost only one environment. (2) In both linear and nonlinear settings, the
balanced situation is the best for DICD, which means we have enough
information from every environment. (3) Even there is only one group, DICD
could make significant improvements against NOTEARS, which coincides with the
discovery in Section 5.3.1.

(a) Linear Setting

(b) Nonlinear Setting

Figure 2: SHD for different percentage of imbalance between two environments
in linear and nonlinear settings.

(a) Ground Truth

(b) NOTEARS

(c) DAGGNN

(d) NOCURL

(e) DICD

Figure 3: Revealed graph by different methods on Real World Data

####  5.3.3 The Effect of Density of Graph

(a) Linear Setting

(b) Nonlinear Setting

Figure 4: SHD for different density conditions.

We aim to discover how much the node degree (_i.e.,_ density) will affect our
algorithm compared to the other baselines. We choose SF graph with $d=20$ for
this case study. In the linear setting, we generate samples from five
environments, 200 samples for each. And for nonlinear settings, we draw
samples from two environments, 1000 samples for each. The x-axis represents
the mean degree of the nodes in the generated graph. For instance, Node degree
= 10 means there are 200 edges in total when generating the SF graph. from
Figure 4, we could observe: (1) DICD almost consistently outperforms the other
methods, except for few cases. (2) As the node degree (_i.e.,_ the density)
increases, especially in nonlinear settings, the improvements of DICD over
baselines get larger, which means DICD could better adapt to the denser
settings.

###  5.4 Real Data (RQ3)

Finally, we evaluate our method on the real-world data: ColoredMNIST. Our
method, as well as NOTEARS-MLP, could be simply extended to the setting of
vector-valued variable. DAG-GNN is already designed to work in this situation.
However, DARING do not discuss this situation in the paper. Though it seems
that DARING could also be generalized to vector-valued, efforts need to be
taken. Thus we omit it as comparison in this section. As stated earlier in
Section 5.1.5, we sample 10000 images from the MNIST dataset and generate the
colored version of MNIST. In the pictures, we have kept the digits to be
black, while the background has the color red or green. The images are resized
to 8*8, and then flattened to yield a 64-dim vector. The label vector and
color vector are set to be all zeros or all ones (according to the original
label) with dimension as 64. Then during searching (or training) we only
consider the possible correlations among the white backgrounds, colors,
noises, and digits. Thus we have four variables in the designed task. As shown
in Figure 3. We could find NOTEARS, DAGGNN and NoCurl all make mistakes, while
DICD gives reasonable predictions.

##  6 Related Work

Causal Discovery has caught enormous attention recently. We will mainly
discuss the differentiable score-based algorithms and the works considering
multi-environments.

Differentiable Score-based algorithms. Score-based causal discovery methods
aim to find the causal structure by optimizing a carefully defined score
function via various modelling methods. Though there are conventional methods
[33, 34, 35, 36, 37] applying various techniques such as hill-climbing[38] and
integer programming [39], the differentiable methods using gradient descent
show stronger power. NOTEARS [1] reformulate the causal discovery problem with
acyclicity constraint as a continuous program. DAG-GNN [7] proposes a variant
of the acyclicity constraint and solves the generalized linear SEM in a graph
autoencoder structure. NOTEARS-MLP [10] and GRAN-DAG [9] extend the continuous
acyclicity regularization into a neural network and achieve better results in
the nonlinear settings. RL-BIC [8] introduces RL to find the DAG with the best
BIC score. DARING [18] proposes to constrain the independence between the
residuals and adopts an adversarial training strategy. DAG-GAN [40] formulates
the problem of DAG structure learning from the perspective of distributional
optimization. [41] shows that applying soft sparsity and DAG constraints would
be enough. There are also other works that apply some alternatives of the
continuous acyclicity regularization. DAG-NOFEAR [42] propose the constraint
term that only depends on the absolute value of the adjacency matrix W,
instead of $\textbf{W}\odot\textbf{W}$ in $h(\textbf{W})$ in Equation (9),
which is correlated to $l1$ penalty and sparsity. NOCURL [12] aims at
eliminating the DAG constraint entirely, by showing that the set of weighted
adjacency matrices of DAGs are equivalent to the set of weighted gradients of
graph potential functions. NODAG [43] proposes to solve an $l_{1}$-penalized
optimization, while ENCO [44] provides convergence guarantees without
constraining the score function with respect to acyclicity. These methods all
have their own strategies to control the acyclicity.

Previous work considering multi-environments. Most work with multi-environment
settings are built on constraint-based methods [45, 46, 21, 47, 48]. Although
they achieved lots of improvements, there are some major limitations: 1) may
have overly strict domain definition [19]; 2) limited to linear cases only
[20, 23, 22]; 3) designed to solve much simpler problems than CD such as
causal direction identification [23, 24]; 4) may involve a large number of
independence tests and be very time-consuming [20, 21], which is also the
general problem of constrain-based methods. However, We define the
environments following [14, 20, 22], and our method is capable of solving the
general causal discovery problems in both linear and nonlinear cases. There is
another interesting and more general task called Federated Causal Discovery
[49, 50, 51, 52, 53], which aims to solve the problem about decentralized
datasets. Though our targets are different, the algorithms are somewhat
similar.

##  7 Conclusion and Future Work

Despite the great success in causal structure learning, today’s differentiable
causal discovery methods are still suffering from non-identifiability issue
and over-reconstruction problem from observation data. Utilizing the inherent
heterogeneity when environment partition is provided in advanced is not yet
discussed in the differentiable causal discovery task. In this paper, we
proposed a simple yet effective Differentiable Invariant Causal Discovery
(DICD) method to tackle the challenge by incorporating the multi-environment
information. Our main idea is that, given the true causal graph, the
parameters of SEM learned from different environments should remain
consistent. Theoretical guarantees for the identifiability of proposed DICD
are provided under certain assumptions about the environments. The extensive
experimental results demonstrate the effectiveness of DICD.

One limitation of DICD is the requirement of prior knowledge of the
environment. In the future study, we will explore the end-to-end causal
discovery technique with environment inference, _i.e.,_ directly infer
partitions of training data without access to environment label. We believe
that the idea of this work, _i.e.,_ invariant learning inspired causal
discovery, provides a potential research direction and will inspire more
valuable works for learning identifiable DAG from observation data.

## References

  * [1] X. Zheng, B. Aragam, P. Ravikumar, and E. P. Xing, “Dags with NO TEARS: continuous optimization for structure learning,” in _NeurIPS_ , 2018, pp. 9492–9503.
  * [2] P. Schwab and W. Karlen, “Cxplain: Causal explanations for model interpretation under uncertainty,” in _NeurIPS_ , 2019, pp. 10 220–10 230.
  * [3] M. N. Vu and M. T. Thai, “Pgm-explainer: Probabilistic graphical model explanations for graph neural networks,” in _NeurIPS_ , 2020.
  * [4] K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan, “Causal protein-signaling networks derived from multiparameter single-cell data,” _Science_ , vol. 308, no. 5721, pp. 523–529, 2005.
  * [5] R. Opgen-Rhein and K. Strimmer, “From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data,” _BMC systems biology_ , vol. 1, no. 1, pp. 1–10, 2007.
  * [6] A. D. Sanford and I. A. Moosa, “A bayesian network structure for operational risk modelling in structured finance operations,” _JORS_ , vol. 63, no. 4, pp. 431–444, 2012.
  * [7] Y. Yu, J. Chen, T. Gao, and M. Yu, “DAG-GNN: DAG structure learning with graph neural networks,” in _ICML_ , vol. 97, 2019, pp. 7154–7163.
  * [8] S. Zhu, I. Ng, and Z. Chen, “Causal discovery with reinforcement learning,” in _ICLR_ , 2020.
  * [9] S. Lachapelle, P. Brouillard, T. Deleu, and S. Lacoste-Julien, “Gradient-based neural DAG learning,” in _ICLR_. OpenReview.net, 2020.
  * [10] X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. P. Xing, “Learning sparse nonparametric dags,” in _AISTATS_ , vol. 108, 2020, pp. 3414–3425.
  * [11] R. Bhattacharya, T. Nagarajan, D. Malinsky, and I. Shpitser, “Differentiable causal discovery under unmeasured confounding,” in _AISTATS_. PMLR, 2021, pp. 2314–2322.
  * [12] Y. Yu and T. Gao, “Dags with no curl: Efficient dag structure learning,” in _NeurIPS 2020_ , 2020.
  * [13] V. Vapnik, “Principles of risk minimization for learning theory,” in _NeurIPS_ , 1991, pp. 831–838.
  * [14] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz, “Invariant risk minimization,” _CoRR_ , 2019.
  * [15] S. Sagawa, P. W. Koh, T. B. Hashimoto, and P. Liang, “Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization,” 2019.
  * [16] R. Cadène, C. Dancette, H. Ben-younes, M. Cord, and D. Parikh, “Rubi: Reducing unimodal biases in visual question answering,” _CoRR_ , vol. abs/1906.10169, 2019.
  * [17] D. Krueger, E. Caballero, J. Jacobsen, A. Zhang, J. Binas, D. Zhang, R. L. Priol, and A. C. Courville, “Out-of-distribution generalization via risk extrapolation (rex),” in _ICML_ , 2021, pp. 5815–5826.
  * [18] Y. He, P. Cui, Z. Shen, R. Xu, F. Liu, and Y. Jiang, “DARING: differentiable causal discovery with residual independence,” in _KDD_ , 2021, pp. 596–605.
  * [19] J. Peters, P. Bühlmann, and N. Meinshausen, “Causal inference by using invariant prediction: identification and confidence intervals,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 78, no. 5, pp. 947–1012, 2016.
  * [20] A. Ghassami, N. Kiyavash, B. Huang, and K. Zhang, “Multi-domain causal structure learning in linear systems,” in _NeurIPS_ , 2018, pp. 6269–6279.
  * [21] B. Huang, K. Zhang, J. Zhang, J. D. Ramsey, R. Sanchez-Romero, C. Glymour, and B. Schölkopf, “Causal discovery from heterogeneous/nonstationary data,” _JMLR_ , vol. 21, pp. 89:1–89:53, 2020.
  * [22] A. Ghassami, S. Salehkaleybar, N. Kiyavash, and K. Zhang, “Learning causal structures using regression invariance,” in _NIPS_ , 2017, pp. 3011–3021.
  * [23] B. Huang, K. Zhang, M. Gong, and C. Glymour, “Causal discovery and forecasting in nonstationary environments with state-space models,” in _ICML_ , 2019, pp. 2901–2910.
  * [24] R. Cai, J. Ye, J. Qiao, H. Fu, and Z. Hao, “FOM: fourth-order moment based causal direction identification on the heteroscedastic data,” _Neural Networks_ , vol. 124, pp. 193–201, 2020.
  * [25] J. Qiao, Y. Bai, R. Cai, and Z. Hao, “Causal discovery from multi-domain data using the independence of modularities,” _Neural Computing and Applicationis_ , vol. 34, no. 3, pp. 1939–1949, 2022.
  * [26] J. Pearl, M. Glymour, and N. P. Jewell, _Causal inference in statistics: A primer_. John Wiley & Sons, 2016.
  * [27] G. B. Huang, M. Mattar, T. Berg, and E. Learned-Miller, “Labeled faces in the wild: A database forstudying face recognition in unconstrained environments,” in _Workshop on faces in’Real-Life’Images: detection, alignment, and recognition_ , 2008.
  * [28] J. Deng, W. Dong, R. Socher, L. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _CVPR_ , 2009, pp. 248–255.
  * [29] P. W. Koh, S. Sagawa, S. M. Xie, M. Zhang, A. Balsubramani, W. Hu, M. Yasunaga, R. L. Phillips, I. Gao, T. Lee _et al._ , “Wilds: A benchmark of in-the-wild distribution shifts,” in _ICML_ , 2021, pp. 5637–5664.
  * [30] P. Bandi, O. Geessink, Q. Manson, M. Van Dijk, M. Balkenhol, M. Hermsen, B. E. Bejnordi, B. Lee, K. Paeng, A. Zhong _et al._ , “From detection of individual metastases to classification of lymph node status at the patient level: the camelyon17 challenge,” _TMI_ , 2018.
  * [31] J. Ni, J. Li, and J. McAuley, “Justifying recommendations using distantly-labeled reviews and fine-grained aspects,” in _EMNLP-IJCNLP_ , 2019\.
  * [32] L. Deng, “The mnist database of handwritten digit images for machine learning research,” _SPM_ , vol. 29, no. 6, pp. 141–142, 2012.
  * [33] D. Heckerman, D. Geiger, and D. M. Chickering, “Learning bayesian networks: The combination of knowledge and statistical data,” _Maching Learning_ , vol. 20, no. 3, pp. 197–243, 1995.
  * [34] I. Tsamardinos, L. E. Brown, and C. F. Aliferis, “The max-min hill-climbing bayesian network structure learning algorithm,” _Machine learning_ , vol. 65, no. 1, pp. 31–78, 2006.
  * [35] J. A. Gámez, J. L. Mateo, and J. M. Puerta, “Learning bayesian networks by hill climbing: efficient methods based on progressive restriction of the neighborhood,” _Data Mining and Knowledge Discovery_ , vol. 22, no. 1-2, pp. 106–148, 2011.
  * [36] S. Nie, D. D. Mauá, C. P. de Campos, and Q. Ji, “Advances in learning bayesian networks of bounded treewidth,” in _NeurIPS_ , 2014, pp. 2285–2293.
  * [37] M. Scanagatta, C. P. de Campos, G. Corani, and M. Zaffalon, “Learning bayesian networks with thousands of variables,” in _NeurIPS_ , 2015, pp. 1864–1872.
  * [38] C. Yuan and B. M. Malone, “Learning optimal bayesian networks: A shortest path perspective,” _JAIR_ , vol. 48, pp. 23–65, 2013.
  * [39] H. Manzour, S. Küçükyavuz, and A. Shojaie, “Integer programming for learning directed acyclic graphs from continuous data,” _CoRR_ , vol. abs/1904.10574, 2019.
  * [40] Y. Gao, L. Shen, and S.-T. Xia, “Dag-gan: Causal structure learning with generative adversarial nets,” in _ICASSP_. IEEE, 2021, pp. 3320–3324.
  * [41] I. Ng, A. Ghassami, and K. Zhang, “On the role of sparsity and DAG constraints for learning linear dags,” in _NeurIPS_ , 2020.
  * [42] D. Wei, T. Gao, and Y. Yu, “Dags with no fears: A closer look at continuous optimization for learning bayesian networks,” in _NeurIPS_ , 2020.
  * [43] G. Varando, “Learning dags without imposing acyclicity,” _CoRR_ , vol. abs/2006.03005, 2020.
  * [44] P. Lippe, T. Cohen, and E. Gavves, “Efficient neural causal discovery without acyclicity constraints,” _CoRR_ , vol. abs/2107.10483, 2021.
  * [45] P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman, _Causation, prediction, and search_. MIT press, 2000\.
  * [46] P. Spirtes, C. Meek, and T. S. Richardson, “Causal inference in the presence of latent variables and selection bias,” in _UAI_ , 1995, pp. 499–506.
  * [47] H. Zhang, K. Zhang, S. Zhou, J. Guan, and J. Zhang, “Testing independence between linear combinations for causal discovery,” in _AAAI_ , vol. 35, no. 7, 2021, pp. 6538–6546.
  * [48] R. Shanmugam, “Causality: Models, reasoning, and inference : Judea pearl; cambridge university press, cambridge, uk, 2000, pp 384, ISBN 0-521-77362-8,” _Neurocomputing_ , vol. 41, no. 1-4, pp. 189–190, 2001.
  * [49] E. Gao, J. Chen, L. Shen, T. Liu, M. Gong, and H. Bondell, “Federated causal discovery,” _CoRR_ , vol. abs/2112.03555, 2021.
  * [50] I. Ng and K. Zhang, “Towards federated bayesian network structure learning with continuous optimization,” _CoRR_ , vol. abs/2110.09356, 2021.
  * [51] R. E. Tillman, D. Danks, and C. Glymour, “Integrating locally learned causal structures with overlapping variables,” in _NIPS_. Curran Associates, Inc., 2008, pp. 1665–1672.
  * [52] S. Triantafilou, I. Tsamardinos, and I. G. Tollis, “Learning causal structure from overlapping variable sets,” in _AISTATS_ , ser. JMLR Proceedings, vol. 9. JMLR.org, 2010, pp. 860–867.
  * [53] S. Triantafillou and I. Tsamardinos, “Constraint-based causal discovery from multiple interventions over overlapping variable sets,” _J. Mach. Learn. Res._ , vol. 16, pp. 2147–2205, 2015.

|  Yu Wang Yu Wang is currently a first-year PhD student at University of
California, San Diego. His research interests lie in Machine Learning and
Natural Language Process. During his undergraduate stage in University of
Science and Technology of Chinan (USTC), he was awarded Baosteel Scholarship,
Huawei Scholarship, Excellent Student Scholarship - Gold, etc.
---|---
|  An Zhang is now a research fellow at National University of Singapore. She
received her Ph.D. degree from the Department of Statistics and Applied
Probability, National University of Singapore, in 2021. Her areas of interest
in research include explainable artificial intelligence, causal representation
learning, differentiable causal discovery, and graph neural networks. She has
publications appeared in several top conferences such as NeurIPS, ICLR, SIGIR.
---|---
|  Xiang Wang is now a professor at the University of Science and Technology
of China (USTC). He received his Ph.D. degree from National University of
Singapore in 2019. His research interests include recommender systems, graph
learning, and explainable deep learning techniques. He has published some
academic papers on international conferences such as NeurIPS, ICLR, KDD, WWW,
SIGIR. He serves as a program committee member for several top conferences
such as KDD, SIGIR, WWW, and IJCAI, and invited reviewer for prestigious
journals such as TKDE, TOIS, TNNLS.
---|---
|  Yancheng Yuan is a research assistant Professor of Department of Applied
Mathematics, The Hong Kong Polytechnic University. His research focuses on the
optimization theory, algorithm design and software development, mathematical
foundation of data science, and data-driven applications. He has published
papers in prestigious journals and conferences, including SIOPT, JMLR, IJAA,
OMS, ICML, WWW.
---|---
|  Xiangnan He is a professor at the University of Science and Technology of
China (USTC). His research interests span information retrieval, data mining,
and multi-media analytics. He has over 90 publications in top conferences such
as SIGIR, WWW, and MM, KDD, and journals including TKDE, TOIS, and TMM. His
work has received the Best Paper Award Honorable Mention in WWW 2018 and SIGIR
2016. He is in the Editorial Board of the AI Open journal, served as the PC
chair of CCIS 2019, the area chair of MM 2019, ECML-PKDD 2020, and the
(senior) PC member for top conferences including SIGIR, WWW, KDD, WSDM etc.
---|---
|  Tat-Seng Chua is the KITHCT Chair Professor at the School of Computing,
National University of Singapore. He was the Acting and Founding Dean of the
School during 1998-2000. Dr Chuas main research interest is in multimedia
information retrieval and social media analytics. In particular, his research
focuses on the extraction, retrieval and question-answering (QA) of text and
rich media arising from the Web and multiple social networks. He is the co-
Director of NExT, a joint Center between NUS and Tsinghua University to
develop technologies for live social media search. Dr Chua is the 2015 winner
of the prestigious ACM SIGMM award for Outstanding Technical Contributions to
Multimedia Computing, Communications and Applications. He is the Chair of
steering committee of ACM International Conference on Multimedia Retrieval
(ICMR) and Multimedia Modeling (MMM) conference series. Dr Chua is also the
General Co-Chair of ACM Multimedia 2005, ACM ICMR 2005, ACM SIGIR 2008, and
ACM Web Science 2015. He serves in the editorial boards of four international
journals. Dr. Chua is the co-Founder of two technology startup companies in
Singapore. He holds a PhD from the University of Leeds, UK.
---|---

TABLE VIII: The specific confounder case in which NOTEARS[1] might make
mistakes. Detailed data generation processes are as follows. Environment 1:
$A=z_{A}(\sim\mathcal{N}(0,4)),B=A/2+z_{B}(\sim\mathcal{N}(0,1)),C=A+B/2+z_{C}(\sim\mathcal{N}(0,4))$.
Environment 2:
$A=z_{A}(\sim\mathcal{N}(0,4)),B=A/2+z_{B}(\sim\mathcal{N}(0,4)),C=A+B/2+z_{C}(\sim\mathcal{N}(0,1))$.
Note that only the variance of $z_{B}$ is different in two environments.

| BAC | BAC | BAC | BAC | BAC | BAC
---|---|---|---|---|---|---
$e_{1}$ | 6.00(1.00, 0.50, 0.50) | 6.05(0.00, 1.25, 0.40) | 8.97(0.00, 0.67, 0.40) | 5.67(0.00, 0.67, 1.50) | 8.97(0.00, 0.67, 0.40) | 5.00(1.00, 1.00, 0.50)
$e_{2}$ | 9.00(1.00, 0.50, 0.50) | 8.00(-0.75, 1.25, 1.00) | 11.22(-0.75, 0.61, 1.00) | 9.96(-0.29, 0.76, 0.90) | 11.56(-0.29, 0.76, 0.55) | 9.20(0.40, 1.00, 0.50)

## Basic Structure Example

In this section, we provide another example shown in Table VIII to demonstrate
our motivation introduced in Section 1. In this table, the blue graphs denote
the ground truth, while graphs with red edges are the wrong edges in the wrong
structures. Below the structures, the first number in front of the parenthesis
refers to the least reconstruction
loss($||\textbf{X}\textbf{W}-\textbf{X}||_{2}$) given the structure in the
same column of Table VIII. Given the structure, we can use the structure to
constrain which of the elements in W should be non-zero elements and others
are not. Then we can learn the optimal parameters for these non-zero elements.
As shown in the Table, the numbers in the parenthesis refer to the optimal
parameters for the edges between A-B, B-C, A-C, respectively. From this table,
we can find that with wrong causal structures which have the potential to
obtain lower reconstruction loss (_i.e.,_ over-reconstruction), the parameters
learned given these structures are not consistent across different
environments. Thus with the condition given by DICD, we can rule out wrong
graphs.

## Notations

We provide the meanings of all the notations used in this paper in Table IX.

TABLE IX: Notations and the corresponding Meanings

Notations | Meanings
---|---
n | number of the data samples
d | number of the variables (or nodes)
X |  $\mathbb{R}^{n\times d}$, observational data
$\textbf{x}_{j}$ |  $\mathbb{R}^{n}$, observational data corresponding to the $j$-th variable
$X_{j}$ |  $j$-th node variable
$F$ | ground truth structure function
$F_{j}$ | ground truth structure function corresponding to $X_{j}$
$f$ | estimated structure function
$f_{j}$ | estimated structure function corresponding to $X_{j}$
Pa($X_{j}$) | parent set of $X_{j}$
Pa(j) | the index set of the parents of $X_{j}$
$z_{j}$ | additive noise of $X_{j}$
$\sigma_{j}^{2}$ | variance of $z_{j}$
$\mathcal{G}(\cdot)$ | corresponding graph of $\cdot$
W |  $\mathbb{R}^{d\times d}$, the coefficients for all edges of $\mathcal{G}(f)$ in linear situations
$m_{l}$ | the number of hidden units in the $l$-th layer in MLP
S |  $\\{0,1\\}^{d\times d}$ invariant structure model

## Proof of Theorem 4.1

To begin with, we provide the notations used in this proof. We denote the set
of all variables as $\\{X_{1},\cdots,X_{d}\\}$, and the vector variable
${[X_{1},\cdots,X_{d}]}^{\top}$ is denoted as X in this section. The set of
parent nodes of $X_{i}$ is denoted as $Pa(X_{i})$, while the index set of
$Pa(X_{i})$ is denoted as $Pa(i)$.

First of all, we prove Lemma 4.1 to show that our DICD could yield the true
coefficient in linear systems under the true causal structure.

###### Proof.

We consider the linear regression coefficient $\hat{\textbf{W}}^{e}$ for every
environment $e$. The regression process in Equation (22) is equivalent to
performing the regression for every node $X_{m},m\in\\{1,\cdots,d\\}$ from its
parents $Pa(X_{m})$.

We incorporate the environment information into the expectation operator
$\mathbb{E}_{e}$. The coefficients of the node $X_{m}$ is denoted as:

$[\hat{\textbf{W}}^{e}]_{m}=arg\min_{\overline{\textbf{w}}}\mathbb{E}_{e}||X_{m}-{(\overline{\textbf{w}}\circ[\textbf{S}_{0}]_{m})}^{\top}\textbf{X}||,$
(25)

where $[\hat{\textbf{W}}^{e}]_{m}$ represents the $m$-th column of
$[\hat{\textbf{W}}^{e}]$ and $[\textbf{S}_{0}]_{m}$ is the $m$-th column of
$[\textbf{S}_{0}]$. Then we can further express $[\hat{\textbf{W}}^{e}]_{m}$
in Equation (25) as:

$[\hat{\textbf{W}}^{e}]_{m}=\mathbb{E}_{e}[X_{m}\Big{(}([\textbf{S}_{0}]_{m}\circ\textbf{X}){[(\textbf{S}_{0}]_{m}\circ\textbf{X})}^{\top}\Big{)}^{-1}[\textbf{S}_{0}]_{m}\circ\textbf{X}].$

For simplicity, we denote $[\textbf{S}_{0}]_{m}\circ\textbf{X}$ as
$\textbf{X}_{\textbf{S}_{0}m}$.

We extract the $m$-th column of $\textbf{W}_{0}$ as $[\textbf{W}_{0}]_{m}$,
and the true generation process for $X_{m}$ could be expressed as:

$X_{m}={[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}+z_{m},$ (26)

where $z_{m}\sim\mathcal{N}(0,(\sigma_{m}^{e})^{2})$. Then as
$[\textbf{W}_{0}]$ is the ground truth coefficients, $[\textbf{W}_{0}]$ should
be consistent with the graph $G_{0}$, as well as the structure
$\textbf{S}_{0}$. Thus we have:
$\textbf{W}_{0}=\textbf{S}_{0}\circ\textbf{W}_{0}$, which indicates
${[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}={[\textbf{S}_{0}]_{m}\circ[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}={[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}_{\textbf{S}_{0}m}$.
Thus the generation process in Equation (26) could be rewritten as:

$X_{m}={[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}_{\textbf{S}_{0}m}+z_{m}.$ (27)

Then we can substitute Equation (27) into Equation (25) to obtain:

$[\hat{\textbf{W}}^{e}]_{m}=\mathbb{E}_{e}[({[\textbf{W}_{0}]_{m}}^{\top}\textbf{X}_{\textbf{S}_{0}m}+z_{m})(\textbf{X}_{\textbf{S}_{0}m}{(\textbf{X}_{\textbf{S}_{0}m})}^{\top})^{-1}\textbf{X}_{\textbf{S}_{0}m}].$

Since $z_{m}$ is an independent additive noise, we can remove $z_{m}$ in the
above equation and yield:

$\displaystyle[\hat{\textbf{W}}^{e}]_{m}$
$\displaystyle=\mathbb{E}_{e}[({[\textbf{w}_{0}]_{m}}^{\top}\textbf{X}_{\textbf{S}_{0}m})(\textbf{X}_{\textbf{S}_{0}m}{(\textbf{X}_{\textbf{S}_{0}m})}^{\top})^{-1}\textbf{X}_{\textbf{S}_{0}m}]$
$\displaystyle=\mathbb{E}_{e}[(\textbf{X}_{\textbf{S}_{0}m}{(\textbf{X}_{\textbf{S}_{0}m})}^{\top})^{-1}\textbf{X}_{\textbf{S}_{0}m}{\textbf{X}}^{\top}_{\textbf{S}_{0}m}]{[\textbf{W}_{0}]}^{\top}_{m}={[\textbf{W}_{0}]}^{\top}_{m}.$
(28)

Note that ${[\textbf{W}_{0}]}^{\top}_{m}\textbf{X}_{\textbf{S}_{0}m}$ is a
scalar. We adopt the rule
$(\textbf{a}\cdot\textbf{b})\textbf{c}=(\textbf{c}\cdot\textbf{b}^{T})\cdot\textbf{a}$,
where $\textbf{a},\textbf{b},\textbf{c}$ are vectors of the same dimension in
the second line. Moreover, since ${[\textbf{W}_{0}]}^{\top}_{m}$ is
independent of the environment $e$ and the variable vector X, it can be moved
out of the expectation operation.

Based on Equation (28), the obtained coefficients in environment $e$,
${[\hat{\textbf{W}}^{e}]}^{\top}_{m}$, is exactly the true causal coefficients
${[\textbf{W}_{0}]}^{\top}_{m}$. Since $m$ is any number in
$\\{1,\cdots,d\\}$, we have $\hat{\textbf{W}}^{e}=\textbf{W}_{0}$. The lemma
is proved. ∎

With Definition 4.1, we further denote the parents of the variable $X_{i}$ in
$G_{0}$ as $Pa_{0}(X_{i})$ and the corresponding indexes as $Pa_{0}(i)$.

Lemma 4.2 demonstrates that the causal directions between any two variables
cannot violate each other in any stable graph and the true causal graph. And
we give the proof as following:

###### Proof.

To prove the lemma, we only need to show that given any
$X_{i}\in\\{X_{1},\cdots,X_{d}\\}$ without being the source node, if
$\exists~{}X_{j}$ such that:

$X_{j}\in Pa_{s}(X_{i}),\quad X_{i}\in Pre_{0}(X_{j}).$ (29)

Then $G_{s}$ cannot be a stable graph.

To achieve this, we need to incorporate the information of the environment
$e$. Given $j$ satisfying Equation (29), we aim to find $j_{0}$ such that:

$X_{j_{0}}\in Pa_{s}(X_{i}),\quad z_{j_{0}}\perp Pa_{s}(X_{i})\backslash
X_{j_{0}}.$

Suppose there exists $j$ satisfying Equation (29), we will find $j_{0}$ with
the following process:
(1) Set $j_{0}$ to be $j$;
(2) If $\exists~{}X_{k}\in Pa_{s}(X_{i})$ such that there exists the path from
$X_{j_{0}}$ to $X_{k}$, then we set $j_{0}$ to be $k$. In this step, the
condition $X_{i}\in Pre_{0}(X_{j_{0}})$ still holds. If there is no such
$X_{k}\in Pa_{s}(X_{i})$ satisfying the above condition, we will exit the
iteration.

In this iteration, we always have $X_{i}\in Pre_{0}(X_{j_{0}})$. Also, with
the above process, we know that $\forall~{}X_{k}\in
Pa_{s}(X_{i})\backslash\\{X_{j_{0}}\\},X_{j_{0}}\notin Pre_{0}(X_{k})$, which
means:

$z_{j_{0}}^{e}\perp X_{k}|\forall k\in Pa_{s}(i)\backslash\\{j_{0}\\}.$ (30)

Now since the result of the regression from $Pa_{s}(X_{i})$ to $X_{i}$ is
exactly the $m$-th column of $\hat{\textbf{W}}^{e}$, we denote the result of
the regression based on Equation (23) as:

$\displaystyle[\hat{\textbf{W}}^{e}]_{m}$
$\displaystyle=arg\min_{{[w_{1},\cdots,w_{d}]}^{\top}}\mathbb{E}_{e}\parallel
w_{j_{0}}X_{j_{0}}$ $\displaystyle+\sum_{k\in
Pa_{s}(m)\backslash\\{j\\}}w_{k}X_{k}-X_{i}\parallel_{2}^{2},$
$\displaystyle\emph{s.t. }w_{j}=0,\quad j\in\\{1,\cdots,d\\}\backslash
Pa_{s}(m).$

Then the optimal value for $w_{j_{0}}$ is:

$\hat{w}_{j_{0}}^{e}=\frac{\mathbb{E}_{e}[X_{i}X_{j_{0}}]-\sum_{k\in
Pa_{s}(i)\backslash\\{j_{0}\\}}\hat{w}^{e}_{k}\mathbb{E}_{e}[X_{k}X_{j_{0}}]}{\mathbb{E}_{e}[X_{j_{0}}^{2}]}.$
(31)

According to the conditions of theorem 4.1, there exist two distinct
environments $e_{1},e_{2}\in\mathcal{E}$ such that
$Var({z_{j_{0}}^{e_{1}}})\neq Var({z_{j_{0}}^{e_{2}}})$ and for any other node
$X_{i}$, _i.e.,_ $i\neq j_{0}$, we have
$Var(z_{i}^{e_{1}})=Var(z_{i}^{e_{2}})$.

Hence in the right term in Equation (31), with $X_{i}\perp z_{j_{0}}$
(according to the condition $X_{i}\in Pre_{0}(X_{j_{0}})$), and $\forall k\in
Pa_{s}(i)\backslash\\{j_{0}\\},z_{j_{0}}^{e}\perp X_{k}$ in Equation (30), we
could have:

$\displaystyle\mathbb{E}_{e_{1}}[X_{i}X_{j_{0}}]$
$\displaystyle=\mathbb{E}_{e_{2}}[X_{i}X_{j_{0}}],$ $\displaystyle\quad\forall
k\in Pa_{s}(X_{i})\backslash\\{j_{0}\\}$
$\displaystyle,\mathbb{E}_{e_{1}}[X_{k}X_{j_{0}}]=\mathbb{E}_{e_{1}}[X_{k}X_{j_{0}}].$

If there exists $k\in Pa_{s}(i)\backslash\\{j_{0}\\}$ such that
$\hat{w}_{k}^{e}$ is different across $e_{1}$ and $e_{2}$, then $G_{s}$ is
already not stable, which raises the contradiction. Otherwise, if
$\hat{w}_{k}^{e_{1}}=\hat{w}_{k}^{e_{2}},\forall k\in
Pa_{s}(i)\backslash\\{j_{0}\\}$, then since $E_{e_{1}}[X_{j_{0}}^{2}]\neq
E_{e_{2}}[X_{j_{0}}^{2}]$, we will have

$\hat{w}_{j_{0}}^{e_{1}}\neq\hat{w}_{j_{0}}^{e_{2}},$ (32)

indicating $G_{s}$ is not a stable graph. The contradiction still exists.

Thus we can have the conclusion: If $G_{s}$ is a stable graph, then $\forall
X_{j}\in Pa_{s}(X_{i})$, we have $X_{i}\notin Pre_{0}(X_{j})$. The lemma is
proved. ∎

Now we give the proof of Theorem 4.2. Since it’s written from Theorem 4.1 with
Definition 4.1, then Theorem 4.1 will also be proved.

###### Proof.

Note that in this case $\textbf{W}_{0}$ is the optimal parameters for $G_{0}$,
and thus it is also the true coefficients for data generation. (as shown in
Lemma 4.1.)

To prove Theorem 4.1, we will split the regression loss in Equation (22) to
the specific loss value in every environment and every node except source
nodes. In other words, we illustrate that for any environment
$e\in\mathcal{E}$,
$\mathcal{L}^{e}(\textbf{S}_{0}\circ\textbf{W}_{0})\leq\mathcal{L}^{e}(\textbf{S}_{s}\circ\textbf{W}_{s})$,
which could directly lead to Equation (24). Hence, in the following part, with
a slight abuse of notation, we will omit the subscript $e$. We will show that
$\forall~{}X_{i}\in\mathcal{V}$, where $\mathcal{V}$ is the node set, the
regression loss from $Pa(X_{i})$ to $X_{i}$ in $G_{s}$ is larger or equal than
in $G_{0}$, with the equaling condition holds when:

$\textrm{i-th column($\textbf{W}_{0}$) = i-th column($\textbf{W}_{s}$)}.$ (33)

Given any $X_{i}\in\mathcal{V},i\in\\{1,\cdots,d\\}$, we denote all the
parents of $X_{i}$ in $G_{s},G_{0}$ as $Pa_{s}(X_{i}),Pa_{0}(X_{i})$ and the
corresponding indexes of its parents as $Pa_{s}(i),Pa_{0}(i)$, respectively.
Then the true generation process for the variable $X_{i}$ could be denoted as:

$X_{i}=\sum_{k\in Pa_{0}(i)}[\textbf{W}_{0}]_{ki}X_{k}+z_{i};\quad
z_{i}\sim\mathcal{N}(0,\sigma_{i}^{2}).$ (34)

When performing regression from $Pa_{s}(X_{i})$ to $X_{i}$, we need to
constrain $\mathcal{G}(\textbf{W}_{s})=G_{s}$, which means

$[\textbf{W}_{s}]_{ki}=0|\forall k\in\\{1,\cdots,d\\}\backslash Pa_{s}(i).$

We denote the minimal loss and the optimal elements as

$\displaystyle L_{i}^{s}=\min_{\\{\overline{w}_{k}|k\in Pa_{s}(i)\\}}$
$\displaystyle\mathbb{E}\parallel X_{i}-\sum_{k\in
Pa_{s}(i)}\overline{w}_{k}X_{k}\parallel_{2}^{2},$ (35)
$\displaystyle\\{w_{k}^{*}|k\in Pa_{s}(i)\\}$
$\displaystyle=arg\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel X_{i}$ $\displaystyle-\sum_{k\in
Pa_{s}(i)}\overline{w}_{k}X_{k}\parallel_{2}^{2},$

where $\\{w_{k}^{*}|k\in Pa_{s}(i)\\}$ is a set and its elements are exactly
the corresponding elements of $[\textbf{W}_{s}]_{i}$, _i.e.,_ the non-zero
elements of $i$-th column of $\textbf{W}_{s}$.

Then we define three sets of indexes as follows:

$\displaystyle\mathcal{I}_{u}=Pa_{0}(i)\cap Pa_{s}(i),$
$\displaystyle\quad\mathcal{I}_{v}=Pa_{0}(i)\backslash Pa_{s}(i),$
$\displaystyle\quad\mathcal{I}_{k}=Pa_{s}(i)\backslash Pa_{0}(i).$

Then after substituting Equation (34) into Equation (35), we have:

$\displaystyle L_{i}^{s}$ $\displaystyle=\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel\sum_{u\in\mathcal{I}_{u}}([\textbf{W}_{0}]_{ui}-\overline{w}_{u})X_{u}$
$\displaystyle+\sum_{v\in\mathcal{I}_{v}}[\textbf{W}_{0}]_{vi}X_{v}+\sum_{k\in\mathcal{I}_{k}}\overline{w}_{k}X_{k}+z_{i}\parallel_{2}^{2}.$

According to Lemma 4.1, $\forall~{}X_{j}\in Pa_{s}(X_{i})$, we have
$X_{i}\notin Pre_{0}(X_{j})$, which means $X_{j}\perp z_{i}$. Thus we have
$z_{i}\perp Pa_{s}(X_{i})$. By definition, we have $z_{i}\perp Pa_{0}(X_{i})$,
which yield $z_{i}\perp\\{Pa_{s}(X_{i})\cup Pa_{0}(X_{i})\\}$. Then we can
rewrite the above equation as:

$\displaystyle L_{i}^{s}$ $\displaystyle=\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel\sum_{u\in\mathcal{I}_{u}}([\textbf{W}_{0}]_{ui}-\overline{w}_{u})X_{u}$
$\displaystyle+\sum_{v\in\mathcal{I}_{v}}[\textbf{W}_{0}]_{vi}X_{v}+\sum_{k\in\mathcal{I}_{k}}\overline{w}_{k}X_{k}\parallel_{2}^{2}+\sigma_{i}^{2}.$

Since the reconstruction loss for $X_{i}$ in the true graph $G_{0}$, denoted
as $L_{i}^{0}$ is exactly the variance of the additive noise added on $X_{i}$,
i.e., $\sigma_{i}^{2}$. Thus here we already have $L_{i}^{s}\geq
L_{i}^{0}=\sigma_{i}^{2}$. In the following part, We will explore the
conditions that make the equation hold exactly.

Note that $\forall j\in Pa_{0}(i)$, we have $[\textbf{W}_{0}]_{ji}\neq 0$.
Then for all the nodes in
$\mathcal{I}_{u}\cup\mathcal{I}_{v}\cup\mathcal{I}_{k}$, there must exists a
node denoted as
$X_{end}\in\mathcal{I}_{u}\cup\mathcal{I}_{v}\cup\mathcal{I}_{k}$ that has no
successors in $\mathcal{I}_{u}\cup\mathcal{I}_{v}\cup\mathcal{I}_{k}$
according to $G_{0}$. This means the noises corresponding to $X_{end}$ will be
independent from all the other nodes in
$\mathcal{I}_{u}\cup\mathcal{I}_{v}\cup\mathcal{I}_{k}$. Then we discuss three
situations for $X_{end}$ to be in $\mathcal{I}_{u}$, $\mathcal{I}_{v}$, and
$\mathcal{I}_{k}$, respectively:

  * •

If $X_{end}\in\mathcal{I}_{u}$, we denote the index for $X_{end}$ as $u_{1}$.
Then we have:

$\displaystyle L_{i}^{s}$ $\displaystyle-
L_{i}^{0}\geq\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel\sum_{u\in\mathcal{I}_{u}\backslash\\{u_{1}\\}}([\textbf{W}_{0}]_{ui}-\overline{w}_{u})X_{u}$
$\displaystyle+\sum_{v\in\mathcal{I}_{v}}[\textbf{W}_{0}]_{vi}X_{v}+\sum_{k\in\mathcal{I}_{k}}\overline{w}_{k}X_{k}\parallel_{2}^{2}+([\textbf{W}_{0}]_{u_{1}i}-\overline{w}_{u_{1}})^{2}\sigma_{u_{1}}^{2}.$

Then the minimum operation $\displaystyle\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}$ in Equation (35) will directly set the optimal value of
$\overline{w}_{u_{1}}$, _i.e.,_ ${w}^{*}_{u_{1}}$ to be
$[\textbf{W}_{0}]_{u_{1}i}$.

  * •

If $X_{end}\in\mathcal{I}_{v}$, we denote the index for $X_{end}$ as $v_{1}$.
Then we have:

$\displaystyle L_{i}^{s}$ $\displaystyle-
L_{i}^{0}\geq\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel\sum_{u\in\mathcal{I}_{u}}([\textbf{W}_{0}]_{ui}-\overline{w}_{u})X_{u}$
$\displaystyle+\sum_{v\in\mathcal{I}_{v}\backslash\\{v_{1}\\}}[\textbf{W}_{0}]_{vi}X_{v}+\sum_{k\in\mathcal{I}_{k}}\overline{w}_{k}X_{k}\parallel_{2}^{2}+[\textbf{W}_{0}]_{v_{1}i}^{2}\sigma_{v_{1}}^{2}.$

In this case, we will have the strict inequality $L_{i}^{s}>L_{i}^{0}$, which
means the reconstruction loss for $G_{s}$ will be strictly larger than
$G_{0}$.

  * •

If $X_{end}\in\mathcal{I}_{k}$, we denote the index for $X_{end}$ as $k_{1}$.
Then we have:

$\displaystyle L_{i}^{s}$ $\displaystyle-
L_{i}^{0}\geq\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}\mathbb{E}\parallel\sum_{u\in\mathcal{I}_{u}}([\textbf{W}_{0}]_{ui}-\overline{w}_{u})X_{u}$
$\displaystyle+\sum_{v\in\mathcal{I}_{v}}[\textbf{W}_{0}]_{vi}X_{v}+\sum_{k\in\mathcal{I}_{k}\backslash\\{k_{1}\\}}\overline{w}_{k}X_{k}\parallel_{2}^{2}+\overline{w}_{k_{1}}^{2}\sigma_{k_{1}}^{2}.$

Then the minimum operation $\displaystyle\min_{\\{\overline{w}_{k}|k\in
Pa_{s}(i)\\}}$ in Equation (35) will directly set the optimal value of
$\overline{w}_{k_{1}}$, _i.e.,_ ${w}^{*}_{k_{1}}$ to be 0.

If $X_{end}\in\mathcal{I}_{v}$, then we can already have
$L_{i}^{s}>L_{i}^{0}$. For the other two situations, we will remove
$X_{end}\in\mathcal{I}_{u}$ or $\mathcal{I}_{k}$ to form the new set
$\mathcal{I}_{u}^{\prime}$ and $\mathcal{I}_{k}^{\prime}$, then find the new
$X_{end}$ in
$\mathcal{I}_{u}^{\prime}\cup\mathcal{I}_{v}\cup\mathcal{I}_{k}^{\prime}$ to
perform the above process again. According to this discussion, we could easily
give the necessary and sufficient conditions for $L_{i}^{s}$ to be equal to
$L_{i}^{0}$:

  * •

$\mathcal{I}_{v}=\emptyset$. If not, then there must be one step of removing
the current $X_{end}$ such that $X_{end}\in\mathcal{I}_{v}$, which will
directly lead to $L_{i}^{s}>L_{i}^{0}$.

  * •

$\forall u\in\mathcal{I}_{u},{w}^{*}_{u}=[\textbf{W}_{0}]_{ui}$.

  * •

$\forall k\in\mathcal{I}_{k},{w}^{*}_{k}=0$.

Now we have the optimal values: $\\{w_{u}|u\in\mathcal{I}_{u}\\}$ and
$\\{w_{k}|k\in\mathcal{I}_{k}\\}$. Besides, according to the definition of
$\mathcal{I}_{u}$ and $\mathcal{I}_{k}$, and $\mathcal{I}_{v}=\emptyset$, we
have: $\mathcal{I}_{u}\cup\mathcal{I}_{k}=Pa_{s}(j)$ and
$\mathcal{I}_{u}=Pa_{0}(i)$. Since $\mathcal{I}_{k}\cap Pa_{0}(i)=\emptyset$,
we have: $\forall k\in\mathcal{I}_{k},[\textbf{W}_{0}]_{ki}=0$.

Since the structure of $\textbf{W}_{s}$ and $\textbf{W}_{0}$ need to be
constrained to be the same as $G_{s}$ and $G_{0}$, respectively, we must have:

$[\textbf{W}_{0}]_{ji}=[\textbf{W}_{s}]_{ji}=0,\quad\forall
j\in\\{1,\cdots,d\\}\backslash(\mathcal{I}_{u}\cup\mathcal{I}_{k}).$

Then according to the above conditions, we have:

$\textrm{i-th column($\textbf{W}_{0}$) = i-th column($\textbf{W}_{s}$)}.$

which is Equation (33). Since $X_{i}$ is any variable except source node in
$\mathcal{V}$, we could have $\textbf{W}_{0}=\textbf{W}_{s}$. (Note that for
the columns corresponding to the source node, all the elements in these
columns should be zero such that $\mathcal{G}(\textbf{W}_{s})=G_{s}$ and
$\mathcal{G}(\textbf{W}_{0})=G_{0}$).

We conclude the proof of Theorem 4.2.

∎