Score-based Causal Discovery from Heterogeneous Data
Authors: Chenwei Ding
Biwei Huang
Mingming Gong
Kun Zhang
Tongliang Liu
Dacheng Tao
Venue:  (2021)
Keywords: causal discovery
heterogeneous data
structure learning

Abstract: Causal discovery has witnessed significant progress over the past decades. Most algorithms in causal discovery consider a single domain with a fixed distribution. However, it is commonplace to encounter heterogeneous data (data from different domains with distribution shifts). Applying existing methods on such heterogeneous data may lead to spurious edges or incorrect directions in the learned graph. In this paper, we develop a novel score-based approach for causal discovery from heterogeneous data. Specifically, we propose a Multiple-Domain Score Search (MDSS) algorithm, which is guaranteed to find the correct graph skeleton asymptotically. Furthermore, benefiting from distribution shifts, MDSS enables the detection of more causal directions than previous algorithms designed for single domain data. The proposed MDSS can be readily incorporated into off-the-shelf search strategies, such as the greedy search and the policy-gradient-based search. Theoretical analyses and extensive experiments on both synthetic and real data demonstrate the efficacy of our method.

Under review as a conference paper at ICLR 2021
SCORE-BASED CAUSAL DISCOVERY FROM HETEROGENEOUS DATA
Anonymous authors Paper under double-blind review
ABSTRACT
Causal discovery has witnessed significant progress over the past decades. Most algorithms in causal discovery consider a single domain with fixed distribution. However, it is commonplace encounter heterogeneous data (data from different domains with distribution shifts). Applying existing methods on such heterogeneous data may lead to spurious edges or incorrect directions the learned graph. In this paper, we develop a novel score-based approach for causal discovery from heterogeneous data. Specifically, we propose a Cultiple-Domainn Score Search (MDSS) algorithm, which is guaranteed to find the correct graph skeleton asymptotically. Furthermore, benefiting from distribution shifts, MDSS enables the detection of more causal directions than previous algorithms designed for single domain data. The proposed MDSS can be readily incorporated into off-the-shelf search strategies, such as the greedy search and the policy- sericc--rrdient-bbeg based search. Theoretical analyses and extensive experiments on both synthetic and real data demonstrate the efficacy of our method.
INTRODUCTION
Discovering causal relations among variables is a fundamental problem in various fields such as economics, biology, drug testing, and commercial decision making. Because conducting randomized controlled trials is usually expensive or even infeasible, discovering causal relations from observational data, i.e. causal discovery (Pearl, 2000; Spirtes al., 2000), has received much attention over the past few decades. Early causal discovery algorithms can be roughly categorized into two types: constraint-based ones (e.g. PC (Spirtes et al., 2000)) and score-based ones (e.g. GES (Chickering, 2002)). In general, these methods cannot uniquely identify the causal graph but are guaranteed to output Markov equivalence class. Since the seminal work by Shimizu et al. (2006), several methods have been developed, achieving identifiability of the whole causal structure by making use of constrained Functional Causal Models (FCMs), including the linear non-Gaussian model (Shimizu et al., 2006), the nonlinear additive noise model (Hoyer et al., 2009), and the post-nonlinear model (Zhang & Hyvärinen, 2009). Recently, Zheng et al. (2018) proposed score-based method that formulates the causal discovery problem as continuous optimization with a structural constraint that ensures acyclicity. Based on the continuous structural constraint, several researchers further proposed to model the causal relations by neural networks (NNs) (Lachapelle et al., 2019; Yu et al., 2019; Zheng et al., 2019). Another recent work Zhu & Chen (2019) used reinforcement learning (RL) for causal discovery, where the RI. agent searches over the graph space and outputs a graph that fits the data best.
The above approaches are designed for data from single domain with a fixed causal model, with the limitation that many of the edge directions cannot be determined without strong functional constraints. In addition, the sample size of data from one domain is usually not large enough guarantee small statistical estimation errors. One way to improve statistical reliability is to combine datasets from multiple domains, such as P-value meta-analyses (Lee, 2015; Marot et al., 2009). The idea of combining multiple-domain data is commonly seen in learning with mixture of Bayesion networks (Thiesson et al., 1998). While mixture of Bayesion networks are usually used for density estimation, the purpose of causal analysis from multiple-domain data is completely different, it aims at discovering the underlying causal graphs for all domains. Regarding causal analysis from multiple-domain data, a challenge is the data heterogeneity problem: the data distribution may vary across domains. For example, in fMRI hippocampus signal analysis. the connection strength among different brain regions may change across different subjects (domains). Due to the distribution shift, directly pooling the data from multiple domains may lead to spurious edges. To tackle the issue, different ways have been investigated, including using sliding windows (Calhoun et al., 2014 4), online change point detection (Adams& & MacKay, 2007), online undirected graph learning (Talih

---

Under review as a conference paper at ICLR 2021
& Hengartner, 2005), locally stationary structure tracker (Kummerfeld & Danks, 2013), and regime aware learning (Bendtsen, 2016). However, these methods may suffer from high estimation variance due to sample scarcity, large type II errors, and large number of statistical tests. Huang al. (2015) recovers causal relations with changing modules by making use of certain types of smoothness of the change, while it does not explicitly locate the changing causal modules. Other similar methods, including Xing et al. (2010); Song et al. (2009), can be reduced to online parameter learning because the causal directions are given.
By utilizing the invariance property (Hoover, 1990; Tian & Pearl, 2001; Peters al., 2016) and the more general independent change mechanism (Pearl, 2000), recently, Ghassami et al. (2018) developed two methods: identical boundaries (IB) and minimal changes (MC), for causal discovery from multi-domain data. However, the proposed methods 1) assume causal sufficiency (i.e., all common causes of variables are measured), which is usually not held in real cireumstances, 2) are designed for linear systems only, 3) and are not capable of identifying causal directions from more than ten domains. Huang et al. (2019) proposed a more general approach called CD-NOD for both linear and nonlinear heterogeneous data, by extending the PC algorithm to tackle the heterogeneity issue. However, inheriting the drawbacks of constraint-based methods, CD-NOD involves a multiple testing problem and is time-consuming due to large number of independence tests.
To overcome the limitations of existing works, we propose a Multiple-Domal Score Search (MDSS) method for causal discovery from heterogencous data, which enjoys the following properties. (1) To avoid spurious edges when combing multi-domain data, MDSS searches over the space of augmented graphs, which includes additional domain index surrogate variable characterize the distribution shift. (2) The changing causal modules can be immediately identified from the recovered augmented graph. (3) Benefiting from causal invariance and the independent change mechanism, MDSS uses novel Multiple-Doman Score (MDS) to help identify more causal directions beyond those in the Markov equivalence class from -sstibuttoon-shifeed shifted data. (4) MDSS can be readily incorporated into off-the-shelf search strategies and is time-efficient and applicable to both linear and nonlinear data. (5) Theoretically, we show that MDSS is guaranteed to find the correct graph skeleton asymptotically, and further identify more causal directions than other traditional score-based and constraint-based algorithms. Empirical studies on both synthetic and real data prove the efficacy of our method.
METHODOLOGY
In this section, we start from a brief introduction to causal discovery and distribution shifts (Section 2.1), and then in Section 2.2 and 2.3, introduce our proposed Mulriple-Domain Score Search (MDSS). In Section 2.2, MDSS starts with a predefined graph search algorithm to learn the skeleton of the causal graph, with the linear Bayesian information criterion (BIC) score or nonlinear generalized score (GS (Huang et al., 2018)) on the augmented causal system. Then in Section 2.3, MDSS further identifies causal directions with Kaaltiple-Domaii Score (MDS) based on the identified skeleton of the graph from Section 2.2. Both theoretically and empirically, we show that MDSS can identify more directions compared to algorithms that are designed for i...d or stationary data.
2.1 BACKGROUND IN CAUSAL DISCOVERY AND DISTRIBUTION SHIFTS
The basic causal discovery problem can be formulated as follows: Suppose there are d observable random variables, i.e. V = (Vi,.... Va). Each random variable satisfies the following generating process: V = fi (PA', where ft is a function to model the causal relation between V/2 and its parents PA', and 6, is a noise variable with non-zero variance. All the noise variables are independent of each other. The task of causal discovery is to recover the causal adjacency matrix B given the observed data matrix X E Txd where By 1 indicates that V, is a parent of Vy, and 7 is the sample size.
We denote the underlying causal graph over V as Go. For each V,, call P(V]|PA)) its causal module. For single domain, the joint probability can be factorized as P(V) = I P(V, PA*). Suppose there are n domains with distribution shifts (i.e. P(V) changes across domains), which implies that some causal modules change across domains. The changes may be caused by the variation of functional models, causal strength, or noise variance. Furthermore, we have the following assumptions.
Assumption 1. The changes of causal modules can be represented as functions of domain index C, denoted by g(C),
2

---

Under review as a conference paper at ICLR 2021
Assumption 2. There is no confounder in each single dataset, but we allow the changes of different causal modules being dependent.
Remark: If changes in several causal modules are dependent, can be regarded as special "confounders" that simultaneously affect these causal modules. As a consequence of such confounders, previous causal discovery algorithms designed for i..d. or stationary data may output erroncous edges. See section 3.1 for an illustration. Thus, causal discovery from multiple-domain data with distribution shifts (i.e. heterogeneous data) can be much more difficult than that from single-domal data.
2.2 SKELETON ESTIMATION ON AUGMENTED GRAPHS
With Assumptions I and 2, is natural to consider g(C") as an extra variable in order to remove any potential influence caused by these special confounders. We assume that there are L such confounders (gr (C), 91 (C)). The causal relation between each observable variable Vi and its parents PA' can be formalized with
Vi = f.
(1)
where (C) S {o(( is the set of confounders that influence Vi, (c(C)) the effective parameters in Vi's causal module that are also assumed to be functions of C' and are mutually independent for all variables.
Let Go be the underlying causal graph over V. We denote the graph resulting from adding arrows (C() + V; and (C) - Vi on Go for each Vú in V as Gang over V U {g(C)) U (0:(C)) We call Gaug an augmented graph (see Figure (dd an example), which satisfies the following assumption.
Assumption 3. The joint distribution over V U (9(C)) U (0,())) is Markov and faithful to Gaug-
To remove the potential influence from confounders and recover causal relations from multiple domains, one way is to perform causal discovery algorithms on the augmented graph. While (g) and {C)) not directly observed, we take C as a surrogate variable (Huang et al., 2019) for them because C always available as a domain index. Given Assumptions 1, 2 and 3, one can apply any score-based method over VU [U(C} to recover the causal relations among variables V as if (c)()) U (0,(C)) were known. For simplicity, we denote the graph over {C} as augmented graph as well. Since the domain index, P(C') follows a discrete uniform distribution. Correspondingly, the generating process of non- -tasttato data can be considered as follows: First we generate random values from P(C), and then we generate data points over V according to the SEM in Equation 1. Finally, generated data points are sorted in ascending order according to the values C' (i.c., data points having the same value of C are regarded as belonging to the same domain). In other words, we observe the distribution P(VIC), where P(V[C) may change across different values of C, resulting in non-stationary data. Note that if we do not include C into the system explicitly, samples of V not i..dd. However, after explicitly including the domain index C into the system, P(V, C' is fixed, and thus the pooled data are i.i.d. samples from distribution P(V,C),
Before stating our main result, first give the definitions of globally consistent scoring criterion and locally consistent scoring criterion, which will be used in the paper.
Definition 1 (Globally Consistent Scoring Criterion). Let D be a dataset consisting T i.i.d. samples from some distribution P(-). Let HH and be any DAGs. A scoring criterion S is globally consistent if the following two properties hold as T - Do:
1. If H contains P and G does not contain P, then S(H,D) > D()
2. If'H and G both contain P, and G contains fewer parameters than H. then S(H, D) < S(9, D). Definition 2 (Locally Consistent Scoring Criterion). Let D be a dataset consisting of i..dd samples from some distribution P..). Let G be any DAG, and let G' be the DAG that results from adding the edge V; - V, on G. A scoring criterion S(G, D) is locally consistent if the following two properties hold as - DO:
1. If IVVPAA then S(G', D) D).
2. LVV.LAPAA then S(g' S D).
'eere, larger score means the corresponding graph is closer to the equivalent class of the true DAG, while the MDS defined in Section 2.3 should be regarded a type of "loss function" which needs be minimized.
3

---

Under review as a conference paper at ICLR 2021
It has been shown that the BIC score and the GS score are both globally and locally consistent (Chickering, 2002: Huang et al., 2018).
The procedure for skeleton estimation augmented graphs is described in Algorithm 1. The predefined graph search algorithms will be discussed in Section 2.4. Apart from the recovered skeleton over V. the changing modules can be detected as well in Step 4 of Algorithm 1. It is important note that we allow causal relations to be either linear or nonlinear. If they are nonlinear, we apply GS as a score function. When they are linear, although we can also use GS, we use linear BIC instead because it is less likely be overfitting for linear data and is computationally more efficient.
Algorithm 1 Skeleton Search on Augmented Graph
Input: n. datasets, each has T observations, d variables and index C.
Output: skeleton S of Gaug's 's subgraph 11 over V, and variables Vc E V that are connected with C.
1: Pool all datasets with an extra surrogate variable to form data matrix X E R=Tx(d++1,
2: Use the predefined graph search algorithm with BIC or GS plus acyclicity constraints to recover the augmented graph. Eliminate any direction Vi + C' in the graph with the prior that any variable Vú does not affect domain index. This step leads to the recovered augmented graph Gaug.
3: Discard the index variable in Gaug to obtain the induced subgraph G1. Discard the directions in G1 and
output the skeleton S of G1-
4: Detect changing causal modules by inspecting Gang recovered in Step 2. and output Vo.
The validity of searching on augmented graph is guaranteed by Theorem 1.
Theorem 1. Let D be the pooling of all datasets, De be the augmented dataset with the domain index an extra random variable. Let Go be the underlying causal graph for the distribution of D over V, Ge be the underlying causal graph for the distribution De over V U C. we denote G'c: the graph after the following modifications on Gc: 1. adding any edges, 2. deleting any edges or 3. reversing any edges that changes the conditional dependence relation of Gc, then we have S(Ge, De) > S(Go,, De), where S is any globally consistent scoring criterion.
Proof of the theorem is given in Appendix A. Intuitively, this theorem means we will obtain an augmented graph that is in the same Markov equivalence class as the true augmented graph if we maximize the score.
2.3 CAUSAL DIRECTION DETERMINATION BY CILLTIPLE-DDMIN--SS
For each variable V E Vc. we prove that is possible determine the directions of edges that connect to VR. We denote Vi as any variable that connects to Vn. There are two possible cases:
1) Vi E Vo. In this case, C Vi= Vi forms an unshielded triple. It is intuitive to incorporate the prior that C + Vi (i.e. change of domain leads to the distribution shift of Vt). There are two possible patterns in this case: C' - Vv + Vi2 and C' + Vv + Ví, which we denote P and P' respectively. For P, the causal mechanisms Pofffect (caus)) is invariant when P(cause) changes. For P., we have the invariance of P(cause) when the causal mechanism P P(effect/cause) changes, which is complementary to the invariance of causal mechanisms. The direction between V₁, and Ví2 can be determined as long as a globally consistent score is used. To be specific, suppose P is the true causal pattern underlying the generative distribution, the score of P will be larger than that of P' if the score function used is globally consistent and decomposable, because compared with P. P' eliminates a conditional independence (C IL V/(V)) that actually holds in the generative distribution. This causal direction determination is not achievable for algorithms designed for stationary data from single domain (because domain index cannot be used as an additional variable in this case). To utilize this prior, we simply eliminate any direction V; - C' as described in Step 2 of Algorithm 1. Figure (()) is a graphical illustration for this case.
2) VI E Vc. In this case, both VE and V12 are connected to C, which is much more difficult than case 1). We propose novel multiple-domain score (MDS) that utilizes the property of independent changes of causal modules to determine causal directions based on the causal skeleton derived from Algorithm 1. To specify the idea, take the two-variable case
example. Here assume the true causal direction is Vi + Vg. Figure 1(b) stands for the case where 0, and 02 are independent. (We drop the notation of domain index C for simplicity). In other words, P PVV1; and P (V: 02) change independently. If the recovered
4

---

Under review as a conference paper at ICLR 2021
direction is reverse (see Figure I(c)), factorize the joint distribution
P (V1, V2; 02 0%) P (V2: (22) P [V2
(2)
where 8f and 8/2 are assumed be sufficient for P(V2) and P(V (V2) respectively. Since V/2 + 1/2 not the true direction, of and 622 are not independent, and they are determined jointly by 01 and 02. Based on this point, the causal direction can be determined by comparing the dependence between f1, Bg and the dependence between 8f and 22 and choose the direction with smaller dependence.
a
(a)
(b)
(c)
(d)
Figure 1: (a) In skeleton estimation on augmented graph, we force the direction be C-VV; if algorithm finds link between them. (b)(c) True causal graph for Vi and 1/2 and the graph with reversed causal direction in the two-variable example. (c) A two-variable example where confounder exists.
For linear systems, the dependence can be described with covariance. 02 and 02 can be easily obtained by regressing V1 and V22 on their parents respectively. We first perform the regression for each domain then calculate the covariance between 0, (C)) and (C)... When there are more than two variables that are connected to C, we denote such set of variables as Ve with cardinality m. For each variable Vi E Vc and its parents P.G S Vc, we calculate the sum of the dependence between parameters of V causal module and the parameters of the causal module of each variable in PAY. To incorporate the minimization of such dependence into score bose- method, propose MDS for linear systems:
(dIn(T)) - 2 2l)())))))+))))))) + & DAGs) +cM(A)++++++++++++ . MDS n
lcou(O P (3)
where n, d, T3 and L, represent the number of domains, the number of variables (here we assume this quantity is the same for all domains), sample size, and the maximized log likelihood for domain i, respectively. Å X2 and A3 are regularization coefficients. A3 fixed to 0.001 in our experiments, A and X2 are adjusted dynamically while training (see Zhu & Chen (2019) for how X and X2 are adjusted). A is the weighted adjacency matrix recovered by the algorithm. m is the number of nonstationary variables. cou(.) the covariance operator. The first term in Equation is the average of BIC on domains, the second and third terms are acyclicity constraints proposed in Zhu & Chen (2019) to narrow down the search space to DAGs. See Appendix A.2 for more details about the acyclicity constraints.
For nonlinear systems, cannot be calculated explicitly. In the two-variable case, the dependence between 0j and 0g can be characterized by the dependence between P(V)) and P(V/ V/1 with the assumption that 0₁ and 8g are sufficient for the corresponding distribution module. To calculate the dependence between P(V)) and P(V2 V/1,, Huang et al. (2019) proposes to first use kernel embeddings of distributions P(Vi) and P(V/ Vi), then measure their dependence with extended Hilbert Schmidt Independence Criterion (HSIC (Gretton et al., 2008)) Reproducing Kernel Hilbert Space (RKHS). When there are more than variables that are connected to C. for each variable Va E Ve and its parents P PAK E Vc. we calculate the dependence between P(PA and P(V/) PAK). We propose corresponding MDS for nonlinear systems by integrating such dependence with GS:
A
MDS ponlinear
Σ (GS.) + X I(G & DAGs) A2h(A)
(SIC (4)
where GS,, is the generalized score for domain i, PVvjPPA and ('PAA are the kemne embeddings of distribu- tions P(VVI and P(PAA) respectively. HSIC((-) is HSIC operator that measures the dependence of two random variables. See Appendix A.3, A.4 and A.5 for brief descriptions of GS, kernel embedding of distributions and HSIC respectively.
5

---

Under review as a conference paper at ICLR 2021
Degeneration issue. If we apply search strategies over the entire space of graphs over V to optimize MDS, the MDS penalty (corresponds to the fourth term in Equation 3 and 4) tends to eliminate any edges between each Ve and its parents because the dependence between Byg and empty set (i.e. no parents for VO) is O. We call this degeneration issue. To tackle this issue, we optimize the MDS score based on the skeleton of graph G1 from Algorithm 1. To be specific, we fix the skeleton S of G₁ and apply search strategies over the space defined by S to optimize the MDS score. In other words, MDS is optimized by only altering the direction of each edge in S. With the solution to degeneration issue, we claim that the proposed MDS can recover more correct directions compared with G1. This is supported by Theorem and Theorem 3. See Appendix A.6 and A.7 for proofs.
Let D be the pooling of n datasets with distribution shifts and Go be the DAG underlying the distribution of D. Let S be the skeleton of Go and G1 (recall G1 is in the same equivalent class as Go. which means they have the same skeleton). Let E, be the set of edges that exist in both Gg and G but have different directions in Go and G Let E2 E E₁ be the set of edges whose left node (or variable) and right node are both nonstationary. Let n and n2 be the cardinality of E₁ and E2.
MDSninear (G,D), let E E Eg be the set of edges Theorem 2. For linear systems. let G = arg min mig
whose directions correctly determined by G*, let n" denote the cardinality of Eğ. Given Eg is not empty and G1, G have the same skeleton S, then G is in the same equivalent class as Go and n n2.
Theorem 3. For nonlinear systems, let G = arg min MDSSonntnear(CC (G, let Eg E E2 be the set of G edges whose directions are correctly determined by G*, let n² denote the cardinality of Eiz. Given E2 is not empty and G1, G* have the same skeleton S, then G is in the same equivalent class as Go and E n n22
Theorem 2 and mainly state that under proper assumptions, the directions of some edges whose left node and right are both nonstationary can be correctly determined by the proposed method.
Confounding case. When the confounder g(C) exists (e.g. Figure (cd), the above approach still works the influence from the confounder is not very strong for the following reason: for the correct direction, the dependence of 0, and 02 would come from the confounder, while for the wrong direction, the dependence would come from the confounder as well as the wrong direction.
2.4 GRAPH SEARCH STRATEGIES
The proposed MDSS can be readily incorporated into off-the-shelff search strategies. In this paper, we adopt the policy-gradient based search strategy (Zhu & Chen, 2019) to search for the optimal causal structure. Compared with other search strategies such as greedy equivalence search (GES (Chickering, 2002)), max-min hill climbing (Tsamardinos et al., 2006), direct search by regarding the weighted graph adjacency matrix as parameters (Zheng et al., 2018: Yu et al., 2019; Lachapelle et al., 2019), the policy-gradient-ben gradient-based search by a reinforcement learning (RL) agent with stochastic policy can determine automatically where to search given the uncertainty information of the learned policy, which gets updated promptly by the stream of reward signals (Zhu & Chen, 2019). The graph search strategy using RL is proven to be better than other search strategies mentioned above empirically.
The idea of causal discovery with RL can be summarized as follows. The algorithm uses an encoder-decoder neural network model to generate directed graphs from the observed data, which are then used to compute rewards consisting of the predefined score function as well as some regularization terms for acyclicity. The encoder-decoder model can be regarded as an "actor" that learns to generate "actions" (i.e., graph adjacency matrices) in actor-critic algorithm, an algorithm commonly used in RL. The reward function can be regarded as the "environment" that evaluates how good the "action" is (i.c., how good the produced graph adjacency matrix fits the observed data). The weights of the encoder-decoder model is trained by policy gradient and stochastic optimization methods. The output of the algorithm is the graph that achieves the best reward during the training process.
To integrate MDS with the policy-gradiens gradient-basee search, we replace the predefined score function in the original paper (where BIC is used) with MDS. Apart from policy-gradient gra--rrrtdn-tbsee based search, we also experiment with greedy equivalence search, details of which can be found in Appendix A.8.
The complete search procedure is described in Algorithm 2.
6

---

Under review as a conference paper at ICLR 2021
Algorithm 2 Multiple-domain Score Search
Input datasets each has '' observations, d variables and index C.
Output causal graph G2 over V.
1: Execute Algorithm 1, input all the datasets and corresponding domain index, output skeleton S and nonstationary variables Vc.
2: Execute the predefined graph search algorithm with MDS in the space defined by S, output G2 over V.
3: Perform any pruning methods on G2 if needed.
3 EXPERIMENTS
In this section, we conduct empirical studies to show the effectiveness of our MDSS method combined with the MDS score. We compare MDSS to some well-known causal discovery algorithms that are designed for i..dd or stationary data from single domain (GES (Chickering, 2002), PC (Spirtes et al., 2000), LiNGAM (Shimizu et al., 2006), NO-TEARS (Zheng et al., 2018) and RL (Zhu & Chen, 2019)) as well as algorithms designed for heterogeneous data from multiple domains (CD-NOD (Huang et al., 2019), MC and IB (Ghassami et al., 2018)). The comparison is made on both synthetic and real data.
We evaluate the estimated graphs using three metrics: True Negative Rate (TNR), True Positive Rate (TPR), and Structural Hamming Distance (SHD, i.e. the smallest number of edge additions, deletions, and reversals to convert the estimated graph into the true DAG). A lower SHD indicates better estimate of the causal graph. For algorithms that output completed partially directed acyclic graph (CPDAG), we randomly choose direction for those undirected edges.
3.1 A TOY EXAMPLE
We use a synthetic toy example to illustrate the influence of confounders g(C) for algorithms (we use RL) designed for homogeneous data, and demonstrate that MDSS can avoid such influence and further identify more directions. See Appendix A.9 for this example.
3.2 SYNTHETIC DATA
In this section, we conduct extensive experiments with MDSS and other causal discovery algorithms on linear and nonlinear synthetic data. We denote as the number of datasets, each has d variables and T observations. We set n E {6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20), E {6, 7,8), 8). T 100 for both linear and nonlinear data. We repeat each setting 20 times with DAGs randomly generated by Erdōs Erdös-Rényi model (ER) with parameter 0.3. Each variable VÚ2 is chosen as nonstationary with probability 0.6. Similar to Section 3.1, linear data generated using linear SEM V, P.P by G,, fix w;; and e; across domains and vary b, if V is chosen as nonstationary. Nonlinear data is generated using nonlinear SEM V; = ft()A,,) + b, + G, fi(-) is randomly picked from (sin()), cos(-), sigmoid(()) and b, varies if VG is nonstationary, e; stays invariant.
We first consider the setting when n 6 and d 10. MDSS, MC, IB, and CD-NOD are tested on data from all domains. GES, PC, LiNGAM, NO-TEARS, and RL are tested on data from all domains as well as data from single domain (the domain randomly chosen). For GES, we use fast GES (FGES (Ramsey et al., 2017)), which is an improved version of the original GES. The mean and standard deviation are reported in Table and 2. As we can see, MDSS outperforms other algorithms on both linear and nonlinear data. The performances of PC, FGES, LiNGAM, NO-TEARS, and RL on pooling of all domain data are worse than that on single domain data. Specifically, despite the minor increase in TPR, their TNR decrease dramatically when data from all domains are used. This phenomenon further proves our proposition that distribution shifts will introduce spurious edges if not properly dealt with.
We further compare the performance of MDSS, IB, MC and CD-NOD by varying d and n. The results are reported in Figure 2. The black curves in figures of row 2 and row 3 are shorter than others because CD-NOD takes too much time to give any result when d > 6 and n > 15. According to the results, MDSS outperforms the others in most cases.
To demonstrate that the proposed MDS contributes to the performance, conduct some ablation studies. To be specific, we keep the directions in step 3 of Algorithm I and output G, (Algorithm 2, or MDS search, is not executed). We use the same experimental setting Table 1. The results (TPR, TNR and SHD) are
7

---

Under review as a conference paper at ICLR 2021
Table 1: Empirical results for MISS, MC, IB and CD-NOD linear and nonlinear data.
MDSS
MC
IB
CD-NOD
TPR
0.98++0.88
0.51+0.37
0,40=0.20
0.60+0.25
Lincar
TNR
0.90:0.11
0.65=0.20
0,55=0.17
0.92:0.09
SHD
1.25++1226
4.35+2.35
1.55+1156
1.70+1119
TPR
0.6540.17
0.18:0.18
0.23=0.19
0.6540.30
Nonlincar
TNR
0.80+0.15
0,78=0.11
0.82=0.09
0.64#0.14
SHD
2.3011139
5.20
4.65=1.31
2.6011.20
Table 2: Empirical results for PC, FGES, LINGAM, NO-TEARS and RL on linear and nonlinear data.
PC
FGES
LINGAM
NO-TEARS
RI.
single
pool
single
pool
single
pool
single
pool
single
pool
TPR 0.65=0.25
0.78+0.21
0.70+0229
0.85+0.17
0.23+0.21
0.20=0:22
0.68=0.18
0.73+0.13
0.1100111
0.87+0.10
Linear
TNR
0.48-00220
0.28++0.22
0.44-0.32
0,4740.22
0.69440.6
0.35--022
0.71+0.100
0.60+0.13
0.77+0.09
0.374-0.22
SHD
6.00=2.00
6.60+2.69
5.60+2.73
5.10#1.51
3.40+0.66
7.40=2.42
4.20=11.11
4.9val.45
1.50=1.16
7.50
TPR
0.25
0.53:10.11
0.3511.32
0.48:00288
0.1310013
0.2810.24
0.51:10.37
0.55
0.52:10.28
0.25:10.18
Nonlincar
TNR
0.43=0.28
0.18:20.6
0.4410.26
0.22#0.45
0.75:0013
0.48=0224
0.85=0.11
0.69:0.10
0.75=0.10
0.67:0.14
SHD
7.30=2245
7,80+3.25
7,30#3.10 10
8.10:33.27
4.30=0.9
7.50=2.58
2.55==1448
3.05+1.14
4.10
5.00++1.00
Results of TPR when d-6
Results of TNN when d=6
Results of SHD when dw6
0.4
1.5
F1
475
CDNCO
1.8
544
number domaloo
Rosults TPR
Results TNR
Results
6.7
-b- CENOC
6.1
CD9C0
34. numberr domaaiss
20
number
Furbbar dommalie
Results IP% die
Results of TNR when dwO
Results of SHD when d=
0.2
0.0
VDS
MES
MES
8
05
1.2
22
Figure 2: Comparison of MDSS, MC, IB and CD-NOD when 72 and d vary.
8

---

Under review as a conference paper at ICLR 2021
PRC
DG
PR.
DG
ERC
CA1
FRC
CA
PHC
Sub
PHc
Sub
(a) MDSS. TPR: 0.57, TNR: 0.88, SHD: 3
(b) MC, TPR: 0.29, TNR: 0.38, SHD; 9
PR
DG
PRC
DG
ERe
CA1
ERC
CAL
PHC
Sub
Pfl
Sub
(c) IB. TPR: 0.29, TNR: 0.38, SHD:
(d) CD-NOD. TPR: 0.43, TNR: 0.38, SHD:
Figure 3: Recovered graphs from MDSS, MC, and CD-NOD on hippocampus data.
0.58 + 0.21, 0.92 = 0.07 and 7.95 + 4.51 for linear data, 0.32 + 0.22, 0.58 + 0.22 and 6.20 + 0.26 for nonlinear data. When compared with the results of MDSS Table 1, it is obvious that MDS search identifies more directions and boost the performance.
3.3 REAL DATA
We apply MDSS to fMRI hippocampus dataset (Poldrack et al., 2015). This dataset records signals from 6 separate brain regions: perirhinal cortex (PRC), parahippocampal cortex (PHC), entorhinal cortex (ERC), subiculum (Sub), CAL, and CA3/Dentate Gyrus (DG) of a single person with resting states in 84 successive domain. We select 10 of them. The results from MDSS, MC. IB and CD-NOD are given in Figure 3. The directions of anatomical ground truth are: PHC - ERC, PRC + ERC, ERC + CA3/DG, CA3/DG + CAI, CAI > Sub, Sub + ERC and ERC - CAI. See Appendix A.10 for illustration of the ground true causal graph. As we can see, MDSS outperforms the others.
days. The records for each day can be regarded
4 CONCLUSIONS
This paper proposes a Multiple-Domale Score Search (MDSS) algorithm for causal discovery from heteroge- neous data. It first performs skeleton learning over the space of augmented graphs. Then a Multiple- Multiple-Domain Score (MDS) used determine causal directions based on the skeleton of the recovered graph. The MDS is proposed based on distribution shifts across domains and the assumption of independent change. Compared with previous methods, MDSS can remove the influence distribution shifts and further recover more causal directions. In future work, we aim to improve MDSS from the following two aspects: 1) Score calculation takes more time than training NNs for searching, so it is essential to optimize score computing to accelerate the entire framework. 2) The current framework of MDSS cannot deal with the more general case where causal directions also change, while this phenomenon does exist in some real- world circumstances.
9

---

Under review as a conference paper at ICLR 2021
REFERENCES
Ryan Prescott Adams and David JC MacKay. Bayesian online changepoint detection. arXiv preprint arXiv:0710.3742, 2007.
Marcus Bendtsen. Regime aware learnin. In Conference on Probabilistic Graphical Models, pp. 1-12, 2016.
Vince D Calhoun, Robyn Miller, Godfrey Pearlson, and Tulay Adal. The chronnectome: time-narying connectivity networks as the next frontier in finri data discovery. Neuron, 84(2):262-27 274, 2014.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):007-554 2002.
AmirEmad Ghassami, Negar Kiyavash, Biwei Huang, and Kun Zhang. Multi-domain causal structure learning in linear systems. In Advances neural information processing systems, pp. 6266-6276, 2018.
Arthur Gretton, Kenji Fukumizu, Choon H Teo, Le Song. Bernhard Schölkopf, and Alex J Smola. A kernel statistical test of independence. In Advances in neural information processing systems, pp. 5-5-592, 2008.
Kevin D Hoover. The logic of causal inference: Econometrics and the conditional analysis of causation. Economics & Philosophy, 6(2):207-234, 1990.
Patrik O Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In NIPS, pp. 689-696, 2009.
Biwei Huang, Kun Zhang, and Bernhard Schölkopf. Identification of time-dependent causal model: A gaussian process treatment. In Twenty-Fourth International Joint Conference on Artificial Intelligence, 2015.
Biwei Huang, Kun Zhang, Yizhu Lin, Bernhard Schölkopf. and Clark Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery Data Mining, pp. 1551-1560, 2018
Biwei Huang, Kun Zhang. Jiji Zhang, Joseph Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from helerngeneoos//oonsttiinnyy data. arXiv preprini 190002190..01672, 2019.
Erich Kummerfeld and David Danks. Tracking time-varying graphical structure. In Advances in neural information processing systems, pp. 1205-12133 2013.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-basce neural dag learning. arXiv preprint 17X.22690. 2019.
Young Ho Lee. Meta-analysis of genetic association studies. Annals of laboratory medicine, 35(3):283-287, 2015.
Guillemette Marot, Jean-Louis Foulley, Claus-Dieter Mayer, and Florence Jaffrézic. Moderated effect size and p-value combinations for microarray meta-analyses. Bioinformatics, 25(20):2692-2699 2009.
Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, New York, NY, USA, 2000. ISBN 0-521-77362-8.
Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Statistical Society: Series B (Statistical Methodooogyy,, 785): 2016.
Russell A Poldrack, Timothy Laumann, Oluwasanmi Koyejo, Brenda Gregory, Ashleigh Hover, Mei- Yen Chen, Krzysztof Gorgolewski, Jeffrey Luci, Sung Jun Joo, Ryan I. Boyd, et al. Long-term neural and physiological phenotyping of a single human. Nature communications, 6(1):1-15, 2015.
Joseph Ramsey, Madelyn Glymour, Ruben Monche and Clark Glymour. million variables and more: the fast greedy equivalence search algorithm for learning high- digensiimnii graphical causal models, with an application to functional magnetic resonance images. International journal of data science and analyties, 3(2):121-129, 2017.
10

---

Under review as a conference paper at ICLR 2021
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non- non-gaussian acyelic model for causal discovery. Journal of Machine Learning Research, 7(Oct):2003-200) 2030, 2006.
Le Song, Mladen Kolar, and Eric P Xing. Time-varying dynamic bayesian networks. In Advances in neural information processing systems, pp. 1732-1740, 2009.
Le Song, Kenji Fukumizu, and Arthur Gretton. Kerel embeddings of conditional distributions: unified kernel framework for nonparametrie inference in graphical models. IEEE Signal Processing Magazine, 30 (4):98-111, 2013.
Peter Spirtes, Clark N Glymour, Richard Scheines, David Heckerman, Christopher Mcek, Gregory Cooper, and Thomas Richardson. Causation, prediction, and search. MIT press, 2000.
Makram Talih and Nicolas Hengartner. Structural learning with time-varying components: tracking the cross-section of financial time series. Journal of the Royal Statistical Society: Series (Statistical Methodology), 67(3):321-341, 2005.
Bo Thiesson, Christopher Meek, David Maxwell Chickering, and David Heckerman. Learning mixtures of bayesian networks. Cooper & Moral, 1998.
J Tian and J Pearl. Causal discovery from changes: bayesian approach, ucla cognitive systems laboratory. Technical report, Technical Report, 2001.
Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
Eric Xing, Wenjie Fu, Le Song, al. A state-space mixed membership blockmodel for dynamic network tomography. The Annals of Applied Statistics, 4(2):535 -66, 2010.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: Dag structure learning with graph neural networks. arXiv preprint arXiv: 10098, 2019.
Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlinear causal model. In UAI, pp. 647-655. AUAI Press, 2009.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. In Advances Neural Information Processing Systems, pp. 9472-9483, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. Learning sparse nonparametric dags. arXiv preprint arXiv: 1909. 13189, 2019.
Shengyu Zhu and Zhitang Chen. Causal discovery with reinforcement learning. 1rXiv:1906.04477, 2019.
arXiv preprint
A APPENDIX
A.1 PROOF OF THEOREM 1
Proof. Let P(V) be the distribution of the data over V. (o(V. C) be the distribution of the augmented data over VUC. According to Markov and faithfulness conditions, Ge is the perfect map of PV,V,CO C)..
1. If G' is the graph after deleting any edges from Gc, then Go contains Po(V, C') while Go does not. According to Definition 1, we have S(GG, Dc) > Dc).
2. If Go is the graph after reversing any edges from Gc, that changes the conditional dependence of Gc, then obviously Ge contains Pc(V,C) while G¹ does not. According to Definition 1, we have S(G,, Dc) > (G),,,))
3. If G' is the graph after adding any edges from Go. Although both G' and Go contain Po(V. C), Go has more parameters. According Definition 1, we have S(Go, De) S(G', Dc).
11

---

Under review as a conference paper at ICLR 2021
A.2 ACYCLICITY CONSTRAINTS
Causal discovery from samples of a joint distribution is challenging combinatorial problem because of the intractable search space super-exponential in the number of graph nodes. Recently, Zheng et al. (2018) formulates the structure learning problem as a purely continuous optimization problem by a new characterization of acyclicity that not only smooth but also exact. Zheng et al. (2018) proposes to measure the "DAG-ness" of a graph by
(5) where A is weighted adjacency matrix and d is the number of node in the graph. Function h(-) satisfies the following properties:
h(A) = tr (e4oA)- -d,
(AA) = 0 if and only if A is acyclic.
The values of h quantify the "DAG-ness" of the graph.
à is smooth.
h and its derivatives are easy to compute.
Further, Zhu & Chen (2019) finds that h (A), which is non-negative, can be small for certain cyclic graphs and its minimum over all non-DAGs is not easy to compute. As consequence, it would require very large penalty weight for (((A) to obtain exact DAGs if only I((A) is used. To address the issue, Zhu & Chen (2019) proposes another acyclicity penalty term I[G & DAGs), which is the indicator function w.r.t. acyclicity to induce exact DAGs. The combination the above two acyclicity constraints can be written as XgI DAGs) + (A), which corresponds to the second and third terms in proposed MDS.
See their original papers Zheng et al. (2018) and Zhu & Chen (2019) for more details.
A.3 GENERALIZED SCORE
We use generalized score (GS (Huang et al., 2018)) as model selection criteria to measure how well the a graph fits the data. Here we give a brief introduction of the calculation of GS.
Assume X is random variable with domain X. and H.x is reproducing kernel Hilbert space (RKHS) on X with continuous feature mapping ox : X + Hx. Similarly we define variable Y, Z with domain ) Z, the corresponding RKHS Hy. Hz and feature mapping dy, Φz. Let 2 == (Y, Z). consider the following two regression functions in RKHS:
xx(()) F1( + U1 (x = F2(2) + U2.
(6)
where F. : 2 > H. and F2 E + Hx. If X I1 2| the following equation holds:
Bz[Varxie 1 X) ] | = x
(7)
which means that it is not useful to incorporate Y predictor of X given Z, the first model (i.e. the model with less complexity) in Equation 6 is preferred.
Cross-validee validated likelihood used express such preference. To perform cross validation, the whole data set D is split into a training set and a test set and repeat this procedure K times, i.e. K-fold cross validation. Let ()*) and ((A)) be the kth training and kth validation respectively. Let D and D(k)) be the data of X; and its parents in training set and validation set respectively. The GS of DAG Gi using cross-valin validatad
likelihood is calculated with
Sev (9h; D) = Σ Sev (X.,
1=1
(8)
(i ( 1
where PAG are parents of X, in graph Gr. F((0) is the regression function estimated from kth training data (R)
D(A) e [ D!!)) is the likelhhodoo evaluated on the kth validation set with the estimated regression function.
Another type of GS based marginal likelihood are also proposed, see Huang al. (2018) for more details.
12

---

Under review as a conference paper at ICLR 2021
A.4 KERNEL EMBEDDING OF DISTRIBUTIONS
According to Equation 4 in our paper, we need to calculate the dependence between distributions P(V and P(P.*). The dependence between random variables can be measured by Hilbert Schmidt Independence Criterion (HSIC), which will be discussed in section. To transform the distribution for data from different domains to a random variable in RKHS, we use kernel embeddings conditional distributions (Song et al., 2013). In the rest of this section, we denote P as X and V Y for simplicity.
Let X' be the domain of X, and (H, be reproducing kernel Hilbert space (RKHS) with a measurable kernel on X. Let (ax) E HH be a continuous feature mapping dx X + H. Similar notations are for variables Y and C. We define the cross-covariance operator Cyx H. - G as Cyx == E (Y)).
The kernel embedding of the conditional distribution P(X|C c,)) for data from given domain C' Co can be calculated as
Px/C=ea = =CxoCoo(cc,)
(9)
The empirical estimate for μ'|C==o, is
1 N
1 $
D + +ar)²
px|C=en
Dcn
N
(10)
= xx (Ko + +1I) ko,cn
where N is sample size, 4x (b1x))......) ...., (IN)), de == (ci)..... ..... & (cN)). Ke (c:; ce)
(6 (ce) @ (ce)), kc,ca = [k: (a), ...,..., (eN, cn)l The corresponding Gram matrix with Gaussian kernel with oz is
diag (Mx) IN + IN diag (Mx) -2M
M = exp
(1I)
20
where diag(-) sets all off-diagonal entries of the matrix zero, and 1N is N N matrix with all entries being 1. M2x is the Gram matrix with a linear kernel:
M = Ke (Ke + +1)= Kx (Ke + A1) Ke,
(12)
whose (c, entry can be calculated by
M'x (c, = p c=c xloze
I xx (Ke AI) = =k (Ke Al)
x1)
ke.
(13)
= = (Ke - Kx (Ke + AT) ke.e'
Similarly we can calculate the empirical kernel embedding of the conditional distribution P(Y C cn) and the corresponding Gram matrix, which we denote as (ayjX,C=e, and M respectively. For more details about kernel embeddings of distributions, see Song et al. (2013) and Huang et al. (2019).
A.5 EXTENDED HILBERT SCHMIDT INDEPENDENCE CRITERION
With the notations and results in the above section, calculate the dependence between P(X) and P(Y X) by extended Hilbert Schmidt Independence Criterion:
(N - 1)2 tr ( f
(14)
where H a matrix used center the features, whose entries Hij 8y - Huang et (2019) uses a normalized version of the estimated HSIC, which is invariant to the scale M and M M y|x
HSICP(X).P(1IX)
1 HSIC P(YXX)
N-1 tr (M N-, tr (MS. (
(15)
Ir (MH
Ir (M ir (M
13

---

Under review as a conference paper at ICLR 2021
A.6 PROOF OF THEOREM 2
Definition 3. Score Equivalence. Let D be dataset consisting Tf records that are i..d samples from some distribution P(.). A score function S is score equivalent for any two DAGs 9 and G' which are in the same Markov equivalence class, we have S D) = =S(G
Proof. Let G" = arg min M DDD (,, G" = arg min DDSnoottinart (G,D) and G""
G
arg argmin min MDS DDDSnntinerrrC D) be three DAGs with the same skeleton S as G1, let Ezz, E and Ex E Eg be the set of edges with cardinality n', n" and n" whose directions are correctly determined by G', G" and G" respectively.
If G' not in the same equivalent class as Go and n' n2- According to score consis- tency of BIC, the first term in M Ditnnear (G', is larger than that in MDStinear (G", and the last term in MDSiineal (G', D) and (Stieerrr"" (G*,D) is same because n' n2. then M DS (G', D) > M DStnneor (G*,D).
If G' is in the same equivalent class as Gn and n' < n2- According to score equivalence BIC, the first term in MDStinear (G', D) and MDStinear (G", D) is same and the last term in M DSinnear ( (G', DD is larger than that in M DD (G", D) because n' n2 then M DStineer (G', > DDSimeearG'.. (G*, D).
If G' is not in the same equivalent class as GD and n' < n2. then M MD (G', D) > M DS (G", D) clearly holds according to the above two cases.
All the above three cases contradict with the condition that these three graphs are DAGs that minimize M L.n.Snotinerrr(C (G, D).
A.7 PROOF OF THEOREM 3
The conclusion that 0 < n² & T2 holds clearly. The conclusion that G is in the same equivalent class as Go can be proved similar to the proof of Theorem 2. M D)S DSmontinear is not able to guarantee n* -22 mainly because GS not score equivalent.
A.8 EXPERIMENTS OF GREEDY EQUIVALENCE SEARCH
The proposed MDSS be readily incorporated into off-the-shelf search strategies. In the main paper, we adopt policy-gradient- -raddint--aeeeeedd search strategy (Zhu & Chen, 2019) to search the optimal causal structure. In this section, we demonstrate that greedy equivalence search (Chickering, 2002) can also be utilized the search strategy.
Similar the two-stage search in our main paper, we first perform greedy equivalence search on the augmented graphs (i.e. graphs with domain index as an additional node) to optimize the score (BIC for linear systems and GS for nonlinear systems). The output of this step is an equivalence class of the augmented graph. Then we utilize distribution shifts of the nonstationary variables to detect more edge directions. Consider the same setting as in Section 3.2 (i.e. n = 6 and d = 10). In linear case, TPR, TNR, SHD for MDSS with greedy equivalence search is 0.67 + 0.15, 0.75 + 0.16 and 4.75 + 2.05 respectively. In nonlinear case, TPR, TNR, SHD for MDSS with greedy equivalence search is 0.69 + 0.13, 0.63 + 0.21 and 5.80 = 2.67 respectively. Compared with the results in Table 1, 1) MDSS with greedy equivalence search outperforms MC, IB and CD-NOD in both linear and nonlinear cases, 2) although MDSS with greedy equivalence search is as good as MDSS with aaiic--gaadeettbasee search in linear case, it achieves comparable results as pplicy-gradient-beee one in nonlinear case.
A.9 A TOY EXAMPLE
The example is consisted of 10 linear datasets with 4 variables, whose underlying causal graph is given in Figure 4(a). We use linear SEM Vi P.A.A + b, + 4 to generate the data. For each variable V,, w; is fixed to and e; is also fixed as a standard Gaussian noise across all datasets. To introduce distribution shifts, we vary by across datasets, here by is chosen to be invariant (i.e. Vg is stationary).
14

---

Under review as a conference paper at ICLR 2021
We first RL single dataset (randomly chosen from 10 datasets) and the pooling of all datasets respectively, with results shown in Figure 4(b) and 4(c). RL misidentifies the direction between V and V/2 in both cases, which is reasonable because RL uses BIC plus acyclicity constraint score function, and BIC is score equivalent. Furthermore, when multiple datasets with distribution shifts are used, RI. outputs erroneous edges (i.e. edges between Vi, V4 and V2, V4) due confounders. Next pooling of all datasets, with result shown in Figure 4(d). MDSS correctly detects variables with changing causal modules. It recovers all the directions between those variables correctly as well.
MDSS on the
(b)
(c)
Figure 4: (a) True causal graph for the v-varialle toy example (b) Graph from RL running on single dataset. (c) Graph from RL running on multiple datasets. (d) Graph from MDSS.
A.10 GROUND TRUE CAUSAL GRAPH FOR REAL DATA
PR
DG
ERc
CA1
PHc
Sub
Figure 5: Ground true causal graph for real data
15


Reviews:
[{"review_id": "TO6pqnJEM6c", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper presents a score-based approach for causal graph discovery from heterogeneous data. The approach is claimed to be guaranteed to find the correct graph skeleton asymptotically, and can detect more causal directions than previous algorithms designed for single domain data.  \n\nPros:\n-The paper addresses an important problem of causal discovery from heterogeneous data.\n\nCons:\nI have doubts about the correctness of many conclusions in the paper. For example, it looks to me that for Theorem 1 to be true, the pooled data $D_C$ must be i.i.d. samples from distribution P(V,C). However, $C$ is defined as the domain index. The problem setup simply assumes we are given $n$ data sets from $n$ domains. Then it’s natural to assume the domain index “C” will stand for 1, 2,  …, n. Under this setting, the pooled data are not i.i.d. samples, Theorem 1 will not hold, and I believe some other conclusions in the paper will not hold either. For Theorem 1 (and other conclusions) to hold, I believe one has to assume the domain index C is a random variable, such that the data sets from different domains are drawn randomly based on some P(C) and P(V|C). However, this may not be a realistic setting, and is not how the problem is set up in the paper. I don’t think the real data sets in the experiments satisfy this setting (not clear how the synthetic data sets are generated).\n\nOverall, I think the paper has mistakes and vote for reject. \n", "extra_scores": {}, "extra_text": {}}, {"review_id": "FcdQcUHhHkZ", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary of the paper: The paper introduces a novel score-based approach for causal discovery from heterogeneous data. The novel score is called Multiple-Domain Score, and can be incorporated into existing search strategies. It is shown that the MDS is guaranteed to find the correct graph skeleton asymptotically.\n\nStrengths: The paper is rather well written. The related work is properly addressed. The proposed criteria are interesting. \n\nWeaknesses: The results are quire incremental. The paper takes a lot from Huang et al., 2019. Assumptions 1 and 2 are somewhat obvious. The paper discusses criteria based on covariance and independence tests, however, their relation to the causal inference is not so straightforward. Algorithm 2 and Section 2.4 need more details.  \n\nQuestions: The authors tell that the proposed criterion can identify more causal directions compared to the state-of-the art. I did not really find any confirmation of this statement.\n\nSection 2.3: Why do we need $\\theta$? I see that it is also the formulation adopted by Huang et al., 2019, but why the 'augmented' graph is needed? Why $\\theta$ are separate nodes if they are (eq. 2 and 3) just parameters of distributions? \n\nSection 2.3: \"the dependence can be described with covariance\": I agree but what about causality? \n\nThe criteria eq. 3 and 4 contain several terms. I wonder, what is the performance (and whether it is very different) of the criteria containing the last terms (the results of the covariance and HSIC) only. What is the impact of other terms?\n\nProbably I missed it: what is h (what operation applied to the adjacency matrix)?\n\nHow the policy-gradient-based search integrates the MDS? (Section 2.4) \n\nIn experiments: is a CPDAG a causal graph? Is it always so? \nThe same for the Bayesian networks (the hill-climbing mentioned in section 2.4). Can we consider that a BN is a causal network? We need at least some assumptions. \n\nI also doubt that a reinforcement learning algorithm without necessary assumptions can be considered as a causal discovery method.\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "7jGBt5Gh-4f", "reviewer": "AnonReviewer5", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes strategies for learning the structure of multiple sets of data observed over a common set of variables which may exhibit distribution shift. The authors address this problem by augmenting the dataset with an indicator variable which indicates membership to  dataset. After augmenting the dataset standard algorithms for structure learning are applied, with the additional restriction that the indicator variable may only be an ancestor. The authors provide theory that shows the procedure consistently estimates the local structures. The authors then show how the additional information obtained from the structure learned with the context variable can be used to disambiguate directions. Experimental results show the efficacy of the proposed approach. \n\nOverall, I think this is a sensible idea and contains some nice results. I do have a couple of questions: \n\n(1) In the paper the authors explicitly limit the algorithm to the case where all domains observe the same variables, however it seems like this need not be the case?\n(2) The authors should cite the literature on learning with mixtures of Bayesian networks. While the aims are different, the model representation is very similar. \n(3) Can the authors provide intuition on whether this degeneracy will also become an issue if there are errors in the skeleton detection algorithm?", "extra_scores": {}, "extra_text": {}}, {"review_id": "3njY0I1L9pk", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Strengths:\n\n1. This paper studies an interesting problem of causal discovery from heterogeneous data, and does a comprehensive survey of related literature.\n\n2. This paper proposes a novel Multiple-Domain Score Search algorithm for score-based causal discovery from different domains.\n\n3. This paper conducts experiments on both synthetic and real-world data with comparison to state-of-the-art baselines, and the results seem promising.\n\nWeaknesses:\n\n1. The difficulty of causal discovery in multiple domains is not clearly shown in the experiments.\n\n2. Some parts of the writing can be improved.\n\n3. There might still be some important aspects of the method left to evaluate.\n\nThis paper studies a problem of causal discovery from heterogeneous data, the problem is well-motivated, and the authors give a comprehensive survey of the related work. This paper proposes a novel algorithm Multiple-Domain Score Search (MDSS) to address the problem, which searches over the space of augmented graphs to avoid the spurious edges, characterizes distribution shift with the additional domain index, and uses a novel Multiple-Domain Score (MDS) to identify causal directions. The paper evaluates the proposed method on synthetic and real-world datasets, it achieves better performance than baselines.\n\nHere are some concerns about this paper:\n\n1. In the experiment, it is not very clear how difficult it is for causal discovery in multiple domains (e.g. the number of spurious edges, the difference of data distributions in multiple domains, especially in the real-world dataset), which brings challenges to convincingly evaluate the superiority of the proposed method.\n\n2. As for the writing, some parts of the paper are not very easy to follow for non-experts. It would be better if the terms (e.g., causal sufficiency, causal modules) can be briefly introduced when first mentioned.\n\n3. Also, I think there may be many interesting aspects to evaluate in the experiments, e.g., how robust the algorithm is under different levels of domain shift; how each component in the method contributes to the performance (some ablation studies may be considered).", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This submission tackles an important problem and presents interesting ideas. I am confident that the research will lead to good publications. However, in the particular situation here, AnonReviewer2 had serious concerns that are shared by me. The authors made a great effort to clarify the situation, but the current situation still leaves me uncertain about the presentation and correctness of everything. Because some issues were major, it is not easy to re-evaluate and take new conclusions in the short time of this process. I hope the authors do not take this too negatively, but given all the comments and discussions, it is best that another round of improvements and reviews be conducted.
Decision: Reject