Ordering-Based Causal Discovery with Reinforcement Learning

1 Introduction

Identifying causal structure from observational data is an important but also challenging task in many practical applications. This task can be formulated as that of finding a Directed Acyclic Graph (DAG) that minimizes a score function defined w.r.t. the observed data. However, searching over the space of DAGs for the best DAG is known to be NP-hard, even if each node has at most two parents [Chickering, 1996]. Consequently, traditional methods mostly rely on local heuristics to perform the search, including greedy hill-climbing and greedy equivalence search that explores the Markov equivalence classes [Chickering, 2002].

∗Work was done during an internship at Huawei Noah's Ark Lab.

†Corresponding author.

Along with various search strategies, existing methods have also cast causal structure learning problem as that of learning an optimal variable ordering, considering that the ordering space is significantly smaller than that of directed graphs and searching over the ordering space can avoid dealing with the acyclicity constraint [Teyssier and Koller, 2005]. Many methods, such as genetic algorithm [Larranaga et al., 1996], Markov chain Monte Carlo [Friedman and Koller, 2003] and greedy local hill-climbing [Teyssier and Koller, 2005], have been exploited as the search strategies to find desired orderings. In practice, however, these methods often cannot effectively find a globally optimal ordering for their heuristic nature.

Recently, with smooth score functions, several gradientbased methods have been proposed by exploiting a smooth characterization of acyclicity, including NOTEARS [Zheng et al., 2018] for linear causal models and several subsequent works, e.g., [Yu et al., 2019; Lachapelle et al., 2020; Ng et al., 2019b; Ng et al., 2019a; Zheng et al., 2020], which use neural networks for modelling non-linear causal relationships. As another attempt, [Zhu et al., 2020] utilize Reinforcement Learning (RL) to find the underlying DAG from the graph space without the need of smooth score functions. Unfortunately, [Zhu et al., 2020] achieved good performance only with up to 30 variables, for at least two reasons: 1) the action space, consisting of directed graphs, is tremendous for large scale problems and is hard to be explored efficiently; and 2) it has to compute scores for many non-DAGs generated during training but computing scores w.r.t. data is generally time-consuming. It appears that the RL-based approach may not be able to achieve a close performance to other gradientbased methods that directly optimize the same score function for large causal discovery problems, due to its search nature.

By taking advantage of the reduced space of variable orderings and the strong search ability of modern RL methods, we propose Causal discovery with Ordering-based Reinforcement Learning (CORL), which incorporates RL into the orderingbased paradigm and is shown to achieve a promising empirical performance. In particular, CORL outperforms NOTEARS, a state-of-the-art gradient-based method for linear data, even with 150-node graphs. Meanwhile, CORL is also competitive with a strong baseline, Causal Additive Model (CAM) method [Bühlmann et al., 2014], on non-linear data models.

Contributions. We make the following contributions in this work: 1) We formulate the ordering search problem as a multistep Markov Decision Process (MDP) and propose to implement the ordering generating process in an effective encoderdecoder architecture, followed by applying RL to optimizing the proposed model based on specifically designed reward mechanisms. We also incorporate a pretrained model into CORL to accelerate training. 2) We analyze the consistency and computational complexity of the proposed method. 3) We conduct comparative experiments on synthetic and real data sets to validate the performance of the proposed methods. 4) An implementation has been made available at https://github. com/huawei-noah/trustworthyAI/tree/master/gcastle. 1

2 RelatedWorks

Besides the aforementioned heuristic ordering search algorithms, [Schmidt et al., 2007] proposed L1OBS to conduct variable selection using ℓ1-regularization paths based on the method from [Teyssier and Koller, 2005]. [Scanagatta et al., 2015] further proposed an ordering exploration method on the basis of an approximated score function so as to scale to thousands of variables. The CAMmethod [Bühlmann et al., 2014] was specifically designed for non-linear additive models. Some recent ordering-based methods such as sparsest permutation [Raskutti and Uhler, 2018] and greedy sparsest permutation [Solus et al., 2017] can guarantee consistency of Markov equivalence class, relying on some conditional independence relations and certain assumptions like faithfulness. A variant of greedy sparest permutation was further proposed in [Bernstein et al., 2020] for the setting with latent variables. In the present work, we mainly work on identifiable cases which have different assumptions from theirs.

In addition, exact algorithms such as dynamic programming [Xiang and Kim, 2013] and integer or linear programming [Bartlett and Cussens, 2017] are also used for causal discovery problem. However, these algorithms usually can only work on small graphs [De Campos and Ji, 2011], and to handle larger problems with hundreds of variables, they usually need to incorporate heuristics search [Xiang and Kim, 2013] or limit the maximum number of parents of each node.

Recently, RL has been used to tackle several combinatorial problems such as the maximum cut and traveling salesman problem [Bello et al., 2016; Khalil et al., 2017; Kool et al., 2019]. These works aim to learn a policy as a solver based on the particular type of combinatorial problems. However, causal discovery tasks generally have different relationships, data types, graph structures, etc., and moreover, are typically off-line with focus on a or a class of causal graph(s). As such, we use RL as a search strategy, similar to [Zoph and Le, 2017; Zhu et al., 2020]. Nevertheless, a pretrained model or policy can offer a good starting point to speed up training, as shown in our evaluation results (cf. Figure 3).

1 The extended version can be found at https://arxiv.org/abs/2105. 06631.

3 Background

3.1 Causal Structure Learning

Let G = (d, V, E) denotes a DAG, with d the number of nodes, V = {v1, · · · , vd} the set of nodes, and E = {(vi, vj)|i, j = 1, . . . , d} the set of directed edges from vi to vj. Each node vj is associated with a random variable Xj. The probability model associated with G factorizes as p(X1, · · · ,Xd) = ∏dj=1 p(Xj|Pa(Xj)), where p(Xj|Pa(Xj)) is the conditional probability distribution for Xj given its parents Pa(Xj) := {Xk|(vk, vj) ∈ E}. We assume that the observed data xj is obtained by the Structural Equation Model (SEM) with additive noises: Xj := fj(Pa(Xj)) + ϵj, j = 1, . . . , d, where fj represents the functional relationship between Xj and its parents, and ϵj's denote jointly independent additive noise variables. We assume causal minimality, which is equivalent to that each fj is not a constant for any Xk ∈ Pa(Xj) in this SEM [Peters et al., 2014].

Given a sample X = [x1, · · · , xd] ∈ Rm×d where xj is a vector of m observations for random variable Xj . The goal is to find a DAG G that optimizes the Bayesian Information Criterion (BIC) (or equivalently, minimum description length) score, defined as

where xkj is the k-th observation of Xj, θj is the parameter associated with each likelihood, and |θj | denotes the parameter dimension. For linear-Gaussian models, p(xkj |Pa(xkj ); θj) = N (xj |θT j Pa(xj), σ2 ) and σ2 can be estimated from the data.

The problem of finding a directed graph that satisfies the ayclicity constraint can be cast as that of finding a variable ordering [Teyssier and Koller, 2005; Schmidt et al., 2007]. Specifically, let Π denote an ordering of the nodes in V , where the length of the ordering |Π| = |V | and Π is indexed from 1. If node vj ∈ V lies in the p-th position, then Π(p) = vj. Notation Π≺vj denotes the set of nodes that precede node vj in Π. One can easily establish a canonical correspondence between an ordering Π and a fully-connected DAG G Π; an example is presented in Figure 1. A DAG G can be consistent with more than one orderings and the set of these orderings is denoted by

Φ(Π)={Π : fully-connected DAGG Π is a super-DAG of G}, where a super-DAG of G is a DAG whose edge set is a superset of that of G. The the search for the true DAG G ∗ can be decomposed to two phases: finding the correct ordering and performing variable selection; the latter is to find the optimal DAG that is consistent with the ordering found in the first step.

3.2 Reinforcement Learning

Standard RL is usually formulated as an MDP over the environment state s ∈ S and agent action a ∈ A, under an (unknown) environmental dynamics defined by a transition probability T (s′ |s, a). Let πφ(a|s) denote the policy, parameterized by φ, which outputs a distribution used to select an action from action space A based on state s. For episodic tasks, a trajectory τ = {st, at}Tt=0 , where T is the finite time horizon, can be collected by executing the policy repeatedly. In many cases, an immediate reward r(s, a) can be received when agent executes an action. The objective of RL is to learn a policy which can maximize the expected cumulative reward along a trajectory, i.e., J(φ) = Eπφ [R0] with R0 = ∑Tt=0 γ t rt(st, at) and γ ∈ (0, 1] being a discount factor. For some scenarios, the reward is only earned at the terminal time (also called episodic reward), and J(φ) = Eπφ [R(τ )] with R(τ ) = rT (sT , aT ).

4 Method

In this section, we first formulate the ordering search problem as an MDP and then describe the proposed approach. We also discuss the variable selection methods to obtain DAGs from variable orderings, as well as the consistency and computational complexity regarding the proposed method.

4.1 Ordering Search as Markov Decision Process

To incorporate RL into the ordering-based paradigm, we formulate the variable ordering search problem as a multi-step decision process with a variable as an action at each decision step, and the order of the selected actions (or variables) is treated as the searched ordering. The decision-making process is Markovian, and its elements are described as follows.

State. One can directly take the sample data xj as the state. However, preliminary experiments (see Appendix A.1) show that it is difficult for feed-forward neural network models to capture the underlying causal relationships directly using observed data as states, and that the data pre-processed by an encoder module is helpful to find better orderings. The encoder module embeds each xj to state sj and all the embedded states constitute the state space S := {s1, · · · , sd}. In our case, we also need an initial state, denoted by s0 (detailed choice is given in Section 4.2), to select the first action. The complete state space would be ˆ S := S ∪ {s0}. We will use ˆ st to denote the actual state encountered at the t-th decision step when generating a variable ordering.

Action. We select an action (variable) from the action space A := {v1, · · · , vd} consisting of all the variables at each decision step, and the action space size is equal to the number of variables, i.e., |A| = d. Compared to the previous RL-based method that searches over the graph space with size O(2d×d ) [Zhu et al., 2020], the resulting action space becomes much smaller.

State transition. The state transition is related to the action selected at the current decision step. If the selected variable is vj at the t-th decision step, then the state is transferred to the state sj ∈ S which corresponds to xj embedded by the encoder, i.e., ˆ st+1 = sj .

Reward. In ordering-based methods, only the variables selected in previous decision steps can be the potential parents of the currently selected variable. Hence, we design the rewards in the following cases: episodic reward and dense reward. In the former case, we calculate the score for a variable ordering Π with d variables as the episodic reward, i.e.,

where T = d - 1 and SBIC has been defined in Equation (1), with Pa(Xj) replaced by the potential parent variable set U(Xj); here U(Xj) denotes the variables associated with the nodes in Π≺vj . If the score function is decomposable (e.g., the BIC score), we can calculate an immediate reward by exploiting the decomposability for the current decision step. That is, for vj selected at time step t, the immediate reward is

This belongs the second case with dense rewards. Here we keep -|θj |/2 log m to make Equation (3) consistent with the form of BIC score.

4.2 Implementation and Optimization with Reinforcement Learning

We briefly describe the neural network architectures implemented in our method, as shown in Figure 2. More details can be found in Appendix A.

Encoder. f enc φe : ˜X ↦→ S is used to map the observed data to the embedding space S = {s1, · · · , sd}. Similar to [Zhu et al., 2020], we adopt mini-batch training and randomly draw n samples from m samples of the data set X to construct ˜X ∈ Rn×d at each episode. We also set the embedding sj to be in the same dimension, i.e., sj ∈ Rn . For encoder choice, we conduct an empirical comparison among several representative structures such as MLP, LSTM and the self-attention based encoder [Vaswani et al., 2017]. Empirically, we validate that the self-attention based encoder in the Transformer structure performs the best (see Appendix A.1).

Decoder. f dec φd : ˆ S ↦→ A maps the state space ˆ S to the action space A. Among several decoder choices (see also Appendix A.1 for an empirical comparison), we pick an LSTM based structure that proves effective in our experiments. Although the initial state is generated randomly in many applications, we pick it as s0 = 1d ∑di=1 si, considering that the source node is fixed in a correct ordering. We restrict each node only be selected once by masking the selected nodes, in order to generate a valid ordering [Vinyals et al., 2015].

Optimization. The optimization objective is to learn a policy maximizing J(φ), where φ = {φe, φd} with φe and φd being parameters associated with encoder f enc and decoder f dec , respectively. Based on the above definition, policy gradient [Sutton and Barto, 2018] is used to optimize the ordering generation model parameters. For the episodic reward case, we have the following policy gradient ∇J(φ) = Eπφ [R(τ ) ∑Tt=0 ∇φ log πφ (at|ˆ st)], and the algorithm in this case is denoted as CORL-1. For the dense reward case, policy gradient can be calculated as ∇J(φ) = Eπφ [∑Tt=0 Rt∇φ log πφ (at|ˆ st)], where Rt = ∑T-t l=0 γ l rt+l denotes the return at time step t. We denote the algorithm in this case as CORL-2. Using a parametric baseline to estimate the expected score typically improves learning [Sutton and Barto, 2018]. Therefore, we introduce a critic network Vφv (ˆ st) parameterized by φv, which learns the expected return given state ˆ st and is trained with stochastic gradient descent using Adam optimizer on a mean squared error objective between its predicted value and the actual return. More details about the critic network are described in Appendix A.2.

Inspired by the benefits from pretrained models [Bello et al., 2016], we also consider to incorporate pretraining to our method to accelerate training. In practice, one can usually obtain some observed data with known causal graphs or correct orderings, e.g., by simulation or real data with labeled graphs. Hence, we can pretrain a policy model with such data in a supervised way and use the pretrained model as initialization for new tasks. Meanwhile, a sufficient generalization ability is desired and we hence include diverse data sets with different numbers of nodes, noise types, causal relationships, etc.

4.3 Variable Selection

One can obtain the causal graph from an ordering by conducting variable selection methods, such as sparse candidate [Teyssier and Koller, 2005], significance testing of covariates [Bühlmann et al., 2014], and group Lasso [Schmidt et al., 2007]. In this work, for linear data models, we apply linear regression to the obtained fully-connected DAG and then use thresholding to prune edges with small weights, as similarly used by [Zheng et al., 2018]. For the non-linear model, we adopt the CAM pruning used by [Lachapelle et al., 2020]. For each variable Xj, one can fit a generalized additive model against the current parents of Xj and then apply significance testing of covariates, declaring significance if the reported p-values are lower that or equal to 0.001. The overall method is summarized in Algorithm 1.

Algorithm 1 Causal discovery with Ordering-based RL.

Require: observed data X, initial parameters φe, φd and φv, two empty buffers D and Dscore, initial value (negative

• infinite) BestScore and a random ordering BestOrdeing. 1: while not terminated do 2: draw a batch of samples fromX, encode them to S and calculate the initial state ˆ s0 3: for t = 0, 1, . . . , T do 4: collect a batch of data 〈ˆ st, at, rt〉 with πφ: D = D ∪ {〈ˆ st, at, rt〉} 5: if 〈vt,Π≺vt , rt〉 is not in Dscore then 6: store 〈vt, Π≺vt , rt〉 in Dscore 7: end if 8: end for 9: update φe, φd, and φv as described in Section 4.2 10: if ∑Tt=0 rt > BestScore then 11: update the BestScore and BestOrdering 12: end if 13: end while
• 14: get the final DAG by pruning the BestOrdering

4.4 Consistency Analysis

So far we have presented CORL in a general manner without specifying explicitly the distribution family for calculating the scores or rewards. In principle, any distribution family could be employed as long as its log-likelihood can be computed. However, whether the maximization of the accumulated reward recovers the correct ordering, i.e., whether consistency of the score function holds, depends on both the modelling choice of reward and the underlying SEM. If the SEM is identifiable, then the following proposition shows that it is possible to find the correct ordering with high probability in the large sample limit.

Proposition 1. Suppose that an identifiable SEM with true causal DAG G ∗ on X = {Xj}dj=1 induces distribution P(X). Let G Π be the fully-connected DAG that corresponds to an ordering Π. If there is an SEM with G Π inducing the same distribution P(X), then G Π must be a super-graph of G ∗ , i.e., every edge in G ∗ is covered in G Π.

Proof. The SEM with G Π may not be causally minimal but can be reduced to an SEM satisfying the causal minimality condition [Peters et al., 2014]. Let ˜ G Π denotes the causal graph in the reduced SEM with the same distribution P(X). Since we have assumed that original SEM is identifiable, i.e., the distribution P(X) corresponds to a unique true graph, ˜ G Π is then identical to G ∗ . The proof is complete by noticing that G Π is a super-graph of ˜ G Π.

Thus, if the causal relationships fall into the chosen model functions and a right distribution family is assumed, then given infinite samples the optimal accumulated reward (e.g., the optimal BIC score) must be achieved by a super-DAG of the underlying graph. However, finding the optimal accumulated reward may be hard, because policy gradient methods only guarantee local convergence [Sutton and Barto, 2018], and we can only apply approximate model functions and also need to assume a certain distribution family for calculating the reward. Nevertheless, the experimental results in Section 5 show that the proposed method can achieve a better performance than those with consistency guarantee in the finite sample regime, thanks to the improved search ability of modern RL methods.

4.5 Computational Complexity

In contrast with typical RL applications, we treat RL here as a search strategy, aiming to find an ordering that achieves the best score. CORL requires the evaluation of the rewards at each episode with O(dm2 + d3 ) computational cost if linear functions are adopted to model the causal relations, which is same to RL-BIC2 [Zhu et al., 2020]. Fortunately, CORL does not need to compute the matrix exponential term with O(d3 ) cost due to the use of ordering search. We observe that CORL performs fewer episodes than RL-BIC2 before the episode reward converges (see Appendix C). The evaluation of Transformer encoder and LSTM decoder in CORL take O(nd2 ) and O(dn2 ), respectively. However, we find that computing rewards is dominating in the total running time (e.g., around 95% and 87% for 30- and 100-node linear data models). Thus, we record the decomposed scores for each variable vj with different parental sets Π≺vj to avoid repeated computations.

5 Experiments

In this section, we conduct experiments on synthetic data sets with linear and non-linear causal relationships as well as a real data set. The baselines are ICA-LiNGAM [Shimizu et al., 2006], three ordering-based approaches L1OBS [Schmidt et al., 2007], CAM [Bühlmann et al., 2014] and A* Lasso [Xiang and Kim, 2013], some recent gradient-based approaches NOTEARS [Zheng et al., 2018], DAG-GNN [Yu et al., 2019] and GraN-DAG [Lachapelle et al., 2020], and the RL-based approach RL-BIC2 [Zhu et al., 2020]. We use the original implementations (see Appendix B.1 for details) and pick the recommended hyper-parameters unless otherwise stated.

We generate different types of synthetic data sets which vary along: level of edge sparsity, graph type, number of nodes, causal functions and sample size. Two types of graph sampling schemes, Erdös-Rényi (ER) and Scale-free (SF), are considered here. We denote d-node ER and SF graphs with on average hd edges as ERh and SFh, respectively. Two common metrics are considered: True Positive Rate (TPR) and Structural Hamming Distance (SHD). The former indicates the probability of correctly finding the positive edges among the discoveries. Hence, it can be used to measure the quality of an ordering, and the higher the better. The latter counts the total number of missing, falsely detected or reversed edges, and the smaller the better.

5.1 Linear Models with Gaussian and Non-Gaussian Noise

We evaluate the proposed methods on Linear Gaussian (LG) with equal variance Gaussian noise and LiNGAM data models, and the true DAGs in both cases are known to be identifiable [Peters and Bühlmann, 2014; Shimizu et al., 2006]. We set h ∈

{2, 5} and d ∈ {30, 50, 100} to generate observed data (see Appendix B.2 for details). For variable selection, we set the thresholding as 0.3 and apply it to the estimated coefficients, as similarly used by [Zheng et al., 2018; Zhu et al., 2020].

Table 1 presents the results for 30- and 100-node LG data models; the conclusions do not change with 50-node graphs, which are given in Appendix D. The performances of ICALiNGAM, GraN-DAG and CAM is also given in Appendix D, and they are almost never on par with the best methods presented in this section. CORL-1 and CORL-2 achieve consistently good results on LiNGAM data sets which are reported in Appendix E due to the space limit.

We now examine Table 1 (the values in parentheses represent the standard deviation across data sets per task). Across all settings, CORL-1 and CORL-2 are the best performing methods in terms of both TPR and SHD, while NOTEARS and DAG-GNN are not too far behind. In Figure 3, we further show the training reward curves of CORL-1 and CORL-2 on 100-node LG data sets, where CORL-2 converges faster to a better ordering than CORL-1. We conjecture that this is because dense rewards can provide more guidance information for the training process than episodic rewards, which is beneficial to the learning of RL model and improves the training performance. Hence, CORL-2 is preferred in practice if the score function is decomposable for each variable. As discussed previously, RL-BIC2 only achieves satisfactory results on graphs with 30 nodes. The TPR of L1OBS is lower than that of A* Lasso, which indicates that L1OBS using greedy hill-climbing with tabu lists may not find a good ordering. Note that the SHD of L1OBS and A* Lasso reported here are the results after applying the introduced pruning method. We observe that the SHDs are greatly improved after pruning. For example, the SHDs of L1OBS decrease from 171.6 (29.5), 588.0 (66.2) and 1964.5 (136.6) to 85.2 (23.8), 215.4 (26.3) and 481.2 (49.9) for ER2 graphs with 30, 50 and 100 nodes, respectively, while the TPRs almost keep the same.

We have also evaluated our method on 150-node LG data models on ER2 graphs. CORL-1 has TPR and SHD being 0.95 (0.01) and 63.7 (9.1), while CORL-2 has 0.97 (0.01) and 38.3 (14.3), respectively. CORL-2 outperforms NOTEARS that achieves 0.94 (0.02) and 50.8 (21.8).

Pretraining. We show the training reward curve of CORL-2pretrain in Figure 3, where the model parameters are pretrained in a supervised manner. The data sets used for pretraining contain 30-node ER2 and SF2 graphs with different causal relationships. Note that the data sets used for evaluation are different from those used for pretraining. Compared to that of CORL-2 using random initialization, a pretrained model can accelerate the model learning process. Although pretraining requires additional time, it is only carried out once and when finished, the pretrained model can be used for multiple causal discovery tasks. Similar conclusion can be drawn in terms of CORL-1, which is shown in Appendix G.

Running time. We also report the running time of all the methods on 30- and 100-node linear data models: CORL1, CORL-2, GraN-DAG and DAG-GNN ≈ 15 minutes for 30-node graphs; CORL-1 and CORL-2 ≈ 7 hours against GraN-DAG and DAG-GNN ≈ 4 hours for 100-node graphs; CAM ≈ 15 minutes for both 30- and 100-node graphs, while L1OBS and A* Lasso ≈ 2 minutes for that tasks; NOTEARS ≈ 5 minutes and ≈ 1 hour for the two tasks respectively; RLBIC2 ≈ 3 hours for 30-node graphs. We set the maximal running time up to 15 hours, but RL-BIC2 did not converge on 100-node graphs, hence we did not report its results. Note that the running time can be significantly reduced by paralleling the evaluation of reward. The neural network based learning methods generally take longer time, and the proposed method achieves the best performance among these methods.

5.2 Non-Linear Model with Gaussian Process

In this experiment, we consider causal relationships with fj being a function sampled from a Gaussian Process (GP) with radial basis function kernel of bandwidth one. The additive noise follows standard Gaussian distribution, which is known to be identifiable [Peters et al., 2014]. We consider ER1 and ER4 graphs with different sample numbers (see Appendix B.2 for the generation of data sets), and we only report the results with m = 500 samples due to the space limit (the remaining results are given in Appendix F). For comparison, only the methods that have been shown competitive for this non-linear data model in existing works [Zhu et al., 2020; Lachapelle et al., 2020] are included. For a given ordering, we follow [Zhu et al., 2020] to use GP regression to fit the causal relationships. We also set a maximum time limit of 15 hours for all the methods for fair comparison and only graphs with up to 30 nodes are considered here, as using GP regression to calculate the scores is time-consuming. The variable selection method used here is the CAM pruning from [Bühlmann et al., 2014].

The results on 10- and 30-node data sets with ER1 and ER4 graphs are shown in Figure 4. Overall, both GraN-DAG and DAG-GNN perform worse than CAM. We conjecture that this is because the number of samples are not sufficient for GraNDAG and DAG-GNN to fit neural networks well, as also shown by [Lachapelle et al., 2020]. CAM, CORL-1, and CORL-2 have similar results, with CORL-2 performing the best on 10-node graphs and being slightly worse than CAM on 30node graphs. All of these methods have better results on ER1 graphs than on ER4 graphs, especially with 30 nodes. We also notice that CORL-2 only runs about 700 iterations on 30-node graphs and about 5000 iterations on 10-node graphs within the time limit, due to the increased time fromGP regression. Nonetheless, the proposed method achieves a much improved performance compared with the existing RL-based method.

5.3 Real Data

The Sachs data set [Sachs et al., 2005], with 11-node and 17-edge true graph, is widely used for research on graphical models. The expression levels of protein and phospholipid in the data set can be used to discover the implicit protein signal network. The observational data set has m = 853 samples and is used to discover the causal structure. We similary use Gaussian Process regression to model the causal relationships in calculating the score. In this experiment, CORL-1, CORL-2 and RL-BIC2 achieve the best SHD 11. CAM, GraN-DAG, and ICA-LiNGAM achieve SHDs 12, 13 and 14, respectively. Particularly, DAG-GNN and NOTEARS result in SHDs 16 and 19, respectively, whereas an empty graph has SHD 17.

6 Conclusion

In this work, we have incorporated RL into the ordering-based paradigm for causal discovery, where a generated ordering can be pruned by variable selection to obtain the causal DAG. Two methods are developed based on the MDP formulation and an encoder-decoder framework. We further analyze the consistency and computational complexity for the proposed approach. Empirical results validate the improved performance over existing RL-based causal discovery approach.

Acknowledgments

Xiaoqiang Wang and Liangjun Ke were supported by the National Natural Science Foundation of China under Grant 61973244.