Causal Discovery with Reinforcement Learning
Authors: Shengyu Zhu
Ignavier Ng
Zhitang Chen
Venue:  (2020)
Keywords: causal discovery
structure learning
reinforcement learning
directed acyclic graph

Abstract: Discovering causal structure among a set of variables is a fundamental problem in many empirical sciences. Traditional score-based casual discovery methods rely on various local heuristics to search for a Directed Acyclic Graph (DAG) according to a predefined score function. While these methods, e.g., greedy equivalence search, may have attractive results with infinite samples and certain model assumptions, they are less satisfactory in practice due to finite data and possible violation of assumptions. Motivated by recent advances in neural combinatorial optimization, we propose to use Reinforcement Learning (RL) to search for the DAG with the best scoring. Our encoder-decoder model takes observable data as input and generates graph adjacency matrices that are used to compute rewards. The reward incorporates both the predefined score function and two penalty terms for enforcing acyclicity. In contrast with typical RL applications where the goal is to learn a policy, we use RL as a search strategy and our final output would be the graph, among all graphs generated during training, that achieves the best reward. We conduct experiments on both synthetic and real datasets, and show that the proposed approach not only has an improved search ability but also allows for a flexible score function under the acyclicity constraint. 

Published as a conference paper at ICLR 2020
CAUSAL DISCOVERY WITH REMENTTCEMENT LEARNING
Shengyu Zhu Ignavier Zg Zhitang Chen
Huawei Noah's Ark Lab SUniversity Toronto (zhushengyu, chenzhitang? @huawei. com $ignaviernge ttr.:rrtooccs..oooooooooo
ABSTRACT
Discovering causal structure among set of variables is fundamental problem in many empirical sciences. Traditional score-based casual discovery methods rely on various local heuristics to search for a Directed Acyclic Graph (DAG) according to a predefined score funetion. While these methods, e.g., greedy equivalence search, may have attractive results with infinite samples and certain model assumptions, they are less satisfactory in practice due to finite data and possible violation assumptions. Motivated by recent advances in neural combinatorial optimization, we propose to use Reinforcement Learning (RL) to search for the DAG with the best scoring. Our encoder-decode model takes observable data input and generates graph adjacency matrices that are used to compute rewards. The reward incorporates both the predefined score function and two penalty terms for enforcing acyclicity. In contrast with typical RL applications where the goal is to learn a policy, we use RL among all graphs generated during training, that achieves the best reward. We conduct experiments on both synthetic and real datasets, and show that the proposed approach not only has an improved search ability but also allows flexible score function under the acyclicity constraint.
search strategy and our final output would be the graph,
INTRODUCTION
Discovering and understanding causal mechanisms underlying natural phenomena are important to many disciplines of sciences. An effective approach is to conduct controlled randomized experiments, which however expensive or even impossible in certain fields such as social sciences (Bollen 1989) and bioinformatics (Opgen-Rhar Rhein and Strimmer, (007). Causal discovery methods that infer causal relationships from passively observable data are hence attractive and have been an important research topic in the past decades (Pearl 2009 Spirtes et al. 2000 Peters et al. 2017).
major class of such causal discovery methods are score-based, which assign a score S(G), typically computed with the observed data, each directed graph G and then search over the space of all Directed Acyclic Graphs (DAGs) for the best scoring:
min S(g), subject to g E DAGs. g
While there have been well-defined score functions such as the Bayesian Information Criterion (BIC) or Minimum Description Length (MDL) score (Schwarz) 1978, Chickering, 2002) and the Bayesian Gaussian equivalent (BGe) score (Geiger and Heckerman 1994). Problem I is generally NP-hard to solve (Chickering 1996 Chickering et al. 2004). largely due to the combinatorial nature of its acyclicity constraint with the number of DAGs increasing super-exponentially in the number of graph nodes. To tackle this problem, most existing approaches rely on local heuristics to enforce the acyclicity. For example, Greedy Equivalence Search (GES) enforees acyclicity one edge at time, explicitly checking for the acyclicity constraint when an edge is added. GES is known to find global minimizer with infinite samples under suitable assumptions (Chickering) 2002 Nandy et al.) 2018). but this is not guaranteed in the finite sample regime. There are hybrid methods, e.g., the max-min hill climbing method (Tsamardinos al. 2006), which use constraint-based approaches reduce
Work was done during an internship Huawei Noah's Ark Lab.

---

Published as a conference paper at ICLR 2020
the search space before applying score-based methods. However, this methodology generally lacks principled way of choosing a problem-spee -ppciific combination of score functions and search strategies.
Recently, Zheng et al. (2018) introduced smooth characterization for the acyclicity. With linear models, Problem ① was then formulated as a continuous optimization problem w.r.t the weighted graph adjacency matrix by picking proper loss function, e.g., the least squares loss. Subsequent works Yu et al. (219)) and Lachapelle et al. (2019) have also adopted the evidence lower bound and the negative log-likelihood as loss functions, respectively, and used Neural Networks (NNs) to model the causal relationships. Note that the loss functions in these methods must be carefully chosen in order to apply continuous optimization methods. Unfortunately, many effective score functions. e.g., the generalized score function proposed by Huang et al. (2018) and the independence based score function from Peters et al. (2014), either cannot be represented closed forms or have very complicated equivalent loss functions, and thus cannot be easily combined with this approach.
We propose to use Reinforcement Learning (RL) to search for the DAG with the best score according to predefined score function, as outlined in Figure [ The insight is that an RL agent with stochastic policy can determine automatically where to search given the uncertainty information of the learned policy, which can be updated promptly by the stream of reward signals. To apply RL to causal discovery, we use an encoder-decoder NN model to generate directed graphs from the observed data, which are then used to compute rewards consisting of the predefined score function well as two penalty terms to enforce acyclicity. We resort to policy gradient and stochastic optimization methods to train the weights of the NNs, and our output is the graph that achieves the best reward, among all graphs generated in the training process. Experiments on both synthetic and real datasets show that our approach has much improved search ability without sacrificing any flexibility choosing score functions. In particular, the proposed approach with BIC score outperforms GES with the same score function on Linear Gon-Gausiian Acyclic Model (LiNGAM) and Gauessan datasets, and also outperforms recent gradient based methods when the causal relationships are nonlinear.
Critic
Actor
rewards
Score function
data
encoder
decoder
encs,
graphs
Figure 1: Reinforcement learning for score-based causal discovery.
2
RELATED WORK
Constraint-basee causal discovery methods first use conditional independence tests to find causal skeleton and then determine the orientations of the edges up the Markov equivalence class, which usually contains DAGs that can be structurally diverse and may still have many unoriented edges. Examples include (un al. (2007): Zhang et al. (2012) that use kernel-based conditional independence criteria and the well-known PC algorithm (Spirtes et al. 2000). This class of methods involve multiple testing problem where the tests are usually conducted independently. The testing results may have conflicts and handling them is not easy, though there are certain works, e.g., [gyttnnnn jtt al. 20014), attempting to tackle this problem. These methods are also not robust as small errors in building the graph skeleton can result in large errors in the inferred Markov equivalence class.
Another class of causal discovery methods are based on properly defined functional causal models. Unlike constraint-basee based methods that assume faithfulness and identify only the Markov equivalence class, these methods are able to distinguish between different DAGs in the same equivalence class, thanks to the additional assumptions on data distribution and/or functional classes. Examples include
LiNGAM (Shimizu et al. 2006 2011). the nonlinear additive noise model (Hoyer et al., 2009 Peters
et al.
2014|20177 and the nost-nonlinear causal model (Zhang and Hyvärinen| (2009).
Besides Yu et al. (2019); Lachapelle et al. (2019), other recent NN based approaches to causal discovery include Goudet et al. (2018) that proposes causal generative NNs functional causal
2

---

Published as a conference paper at ICLR 2020
modeling with prior knowledge of initial skeleton of the causal graph and Kalainathan et al. (2018) that learns causal generative models in an adversarial way but does not guarantee acyelicity.
learning (Sutskever et al. 2014) have motivated the use of NNs for optimization in various domains (Vinyals et al. 2013 Zoph and Le 2017] Chen et al. 2017). A particular example is the traveling salesman problem that revisited in the work of pointer networks (Vinyals et al.) 2015). Authors proposed a recurrent NN with nonparametric softmaxes trained in supervised manner to predict the sequence of visited cities. [Bello et al.) (2016) further proposed to use the RL paradigm to tackle the combinatorial problems due to their relatively simple reward mechanisms, It was shown that an RL agent can have a better generalization even when the optimal solutions are used as labeled data in the previous supervised approach. Alternatively, the RL based approach in Dai et al. (017) considered combinatorial optimization problems on (undirected) graphs and achieved a promising performance by exploiting graph structures, in contrast with the general seneence-to-seueene modeling.
Recent advances in seuuence-to-squueee
There are many other successful RL applications in recent years, e.g., AlphaGo (Siverr et al. 2017 where the goal is to learn a policy for given task. As an exception, Zoph and Le (2017) applied RL to neural architecture search. While we use a similar idea as the RL paradigm can naturally include the search task, our work is different in the actor and reward designs: our actor is an encoder-decoder model that generates graph adjacency matrices (cf. Section 4 and the reward is tailored for causal discovery by incorporating a score function and the acyclicity constraint (cf. Section 3.Tp.
Zoop
MODEL. DEFINITION
We assume the following model for data generating procedure, as in Hoyer et al. (2009); Peters et al. (2014). Each variable x; is associated with a node i in d-node DAG G, and the observed value of a;; is obtained as a function of its parents in the graph plus an independent additive noise n2, i.c.,
== Ji (Xpa(()) ) + n,, i = 1,2...., d,
where Xpa(t) denotes the set of variables Iy so that there is an edge from xj to I; in the graph, and the noises n, are assumed to be jointly independent. We also assume causal minimality, which in this case reduces to that each function f: is not a constant in any of its arguments (Peters al.) 2014). Without further assumption on the forms of functions and/or noises, the above model can be identified only up Markov equivalence class under the usual Markov and faithful assumptions (Spirtes et al. 2000 Peters et al. 2014]; in our experiments we will consider synthetic datasets that are generated from fully identifiable models so that it practically meaningful evaluate the estimated graph w.r.t. the true DAG. If all the functions f; are linear and the noises Tı; are Gaussian distributed, the above model yields the class of standard [189): and Heckerman| (1994); (Spirtes et al.) (2000); (Pterss et al. (2017). When the functions are linear but the noises are non-Gaussian, one can obtain the LiNGAM described in Shimizu et al.) (2006) 2011) and the true DAG can be uniquely identified under favorable conditions.
Gaassiann model that has been studied in Bollen
In this paper, we consider that all the variables x; are scalars; extending to more complex cases is straightforward, provided with a properly defined score function. The observed data X, consisting of a number of vectors x == [x1, 22. x=]" E Rd, are then sampled independently according to the above model on unknown DAG, with fixed functions f. and fixed distributions for n., The objective of causal discovery is to use the observed data X, which gives the empirical version of the joint distribution of to infer the underlying causal DAG g.
4 NEURAL NETWORK ARCHITECTURE FOR GRAPH GENERATION
Given a dataset X = [x* 1m 1== where x denotes the k-th observed sample, we want to infer the causal graph that best describes the data generating procedure. We would like to use NNs to infer the causal graph from the observed data; specifically, we aim to design an NN based graph generator whose input is the observed data and the output is graph adjacency matrix. A naive choice would be using feed-forward NNs to output df scalars and then reshape them to adjacency matrix in dd However, this NN structure failed to produce promising results, possibly because the forward NNs could not provide sufficient interactions amongst variables capture the causal relations.
3

---

Published as a conference paper at ICLR 2020
Motivated by recent advances in neural combinatorial optimization, particularly the pointer networks (Bello et al. 2016 Vinyals et al. 2015). we draw nt random samples (with replacement) [=>= from X and reshape them as s = {x, )i=1 where x; E 12" is the vector concatenating all the i-th entries of the vectors [x' In an analogy to the traveling salesman problem, this represents a sequence of cities lying in an r-dim space. We are concerned with generating a binary adjacency matrix A E {0, 1dxa so that the corresponding graph is acyelic and achieves the best score. this work consider encoder-decode
models for graph generation:
Encoder We use the attention based encoder in the Transformer structure proposed by Vaswani et al. (017)) We believe that the self-attention scheme, together with structural DAG constraint, is capable of finding the causal relations amongst variables. Other attention based models such as graph attention network (Veličković al., (0118) may also be used, which will be considered in future work. Denote the outputs of the encoder by enci, i = 1,2,..., d, with dimension de.
Decoder Our decoder generates the graph adjacency matrix in element-wise manner, by building relationships between two encoder outputs enc; and enej. We consider the single layer decoder
Pij (W1: W2, u" tanh(W1 enc; W2 ency),
where Wi, W2 E Rddxxe, E Rdd. x1 are trainable parameters and d,, is the hidden dimension associated with the decoder. To generate a binary adjacency matrix A, pass each entry Piy into logistic sigmoid function (-)) and then sample according to a Bernoulli distribution with probability (gi)) that indicates the probability of existing an edge from z; to zy. To avoid self-loops, we simply mask the (i, i)-th entry in the adjacency matrix.
Other decoder choices include the neural tensor network model (Socher et al. 2013) and the bilinear model that build the pairwise relationships between encoder outputs. Another choice is the Trans- former decoder which can generate an adjacency matrix in a row-wise manner. Empirically, we find that the single layer decoder performs the best, possibly because it contains less parameters and is easier to train to find better DAGs while the self-attention based encoder has provided sufficient interactions amongst the variables for causal discovery. Appendix [A] provides more details regarding these decoders and their empirical results with mncar-Gaussian data models.
REINFORCEMENT LEARNING FOR SEARCH
In this section, we use RL as our search strategy to find the DAG with the best score, as outlined in Figure As one will see, the proposed method possesses an improved search ability over traditional score-based methods and also allows flexible score functions subject to the acyclicity constraint.
5.1 SCORE FUNCTION, ACYCLICITY, AND REWARD
Score Function In this work, we consider only existing score functions to construct the reward that will maximized by an RI agent. Often score-based methods assume a parametric model for causal relationships (e.g., a set of parameters Among all score functions that can be directly included here, we focus on the BIC score that is not only consistent (Haughton et al. 1988) but also locally consistent for its decomposability (Chickering 1996).
equations or multinomial distribution), which introduces
The BIC score for given directed graph 9 is
Sgi(()) = -oog p (X; + de logm,
where i is the maximum likelihood estimator and de denotes the dimensionality of the parameter f. We assume i.i.d. Gaussian additive noises throughout this paper. If we apply linear models to each causal relationship and let a be the corresponding estimate for = the i-th entry in the k-th observed sample, then we have the BIC score being (up to some additive constant)
SBic(G) Σ (m log (oSSR(mm + # #edges) logm, m,
where RSS, Ekn (xf 2²)² denotes the residual sum of squares for the i-th variable. The first term in Eq. 2 equivalent to the log-likelihood objective used by GraN-DAG (Lachapelle et al.
4

---

Published as a conference paper at ICLR 2020
2019) and the second term adds penalty on the number of edges in the graph G. Further assuming that the noise variances are equal (despite the fact that they may be different), we have
d Σ
SBIC (9) = -mddlo
RSS,
/(md) + (cdges) log
(3)
We notice that >, RSS; is the least squares loss used in NOTEARS (Zheng et al. 2018 Besides assuming linear models, other regression methods can also be used to estimate arf. In Section 6 will use quadratie regression and Gaussian Process Regression (GPR) to model causal relationships based on the observed data.
we
Acyelicity A remaining issue is the acyclicity constraint. Other than GES that explicitly checks for acyclicity each time an edge is added, we add penalty terms w.r.t. acyclicity to the score function to enforce acyelicity in implicit manner and allow the generated graph to change more than one edges at each iteration. In this work, we use a recent result from Zheng et al. (2018): a directed graph with binary adjacency matrix A is acyclic if and only if
((A) trace (c*) =0
(4)
where eA is the matrix exponential of A. We find that h(A), which is non-negative, can be small for cyclic graphs and the minimum over all non-DAGs is not easy to find. We would require very large penalty weight to guarantee acyclicity if only h(A) is used. We thus add another penalty term, the indicator function w.r.t. acyclicity, to induce exact DAGs. We remark that other functions (e.g., the total length of all cyclic paths in the graph), which compute some 'distance' from a directed graph to DAGs and need not be smooth, may also be used to construct the acyclicity penalty in approach.
Reward Our reward incorporates both the score function and the acyclicity constraint:
reward == - (S(9) Ig & DAGs) + )²h(AA),
(5)
where I(-) denotes the indicator function and ÅI, A2 2 0 are two penalty parameters. It is not hard to that the larger A and A2 are, the more likely generated graph with high reward is acyelic. We then aim to maximize the reward over all possible directed graphs, or equivalently, we have
min (S(g) + I & DAGs) + )gh(A)). G
(6)
An interesting question is whether this new formulation is equivalent to the original problem with hard acyclicity constraint. Fortunately, the following proposition guarantees that Problems ① and (6) are equivalent with properly chosen Å₁ and X2 which can be verified by showing that minimizer of one problem is also a solution to the other. A proof is provided in Appendix B for completeness.
Proposition 1. Let hmin > O be the minimum of h(A) over all directed cyclic graphs. i.e., hmin min GÆDAGs (A).. Let S* denote the optimal score achieved by some DAG in Problem (7). Assume that St. E R is a lower bound of the score function over all possible directed graphs, iee, St. < min g SG)), and Sy E R is an upper bound on the optimal score with S I Sy. Then Problems and (6) are equivalent if
A1 + Xzhmin 2 Su -S-
For practical use, we need to find respective quantities in order to choose proper penalty parameters. An upper bound Sy can be easily found by drawing some random DAGS or using the results from other methods like NOTEARS, lower bound St depends on the particular score function. With BIC score, we can fit each variable r; against all the rest variables, and use only the RSS; terms but ignore the additive penalty on the number of edges. With the independence based score function proposed by Peters et al. (2014), we may simply set St. = 0. The minimum term Amin as previously mentioned, may not easy to find. Fortunately, with A = Sy SS,, Proposition guaarntees the equivalence of Problems (1) and (6) for any A 2 O. However, simply setting X2 0 could only get good performance with very small graphs (see a discussion in Appendix (). We will pick a relatively small value for A2 which helps to generate directed graphs that become closer to DAGs.
Empirically, we find that the initial penalty weights are set too large, the score function would have little effect on the reward, which then limits the exploration of the RL agent and usually results in DAGs with high scores. Similar to Lagrangian methods, we can start with small penalty weights and gradually increase them so that the condition in Proposition T is satisfied. Meanwhile, notice that
5

---

Published as a conference paper at ICLR 2020
Algorithm 1 The proposed RL approach to score-based causal discovery
Require: score parameters: St. St. and So: penalty parameters: Å1, At, 22 22. and Ag; iteration number for parameter update: tg-
1: for t 1,2..... do
2: Run actor-critic algorithm, with score adjustment by S + So(S Sr.)(Su - Sr.)
3: if (mod (o) = then
4:
if the maximum reward corresponds to a DAG with score Smin then update Sy + min(Su, Smin)
5:
6:
end if
7: update A t min(A1 + A Sy) and A2 + min(A2222 A2)
8: update recorded rewards according to new X1 and X2
9: end if
end for
different score functions may have different ranges while the acyclicity penalty terms independent of the particular range of the score function. We hence also adjust the predefined scores to a certain range by using So(S- St) (SS -S,) ) for some So > 0 and the optimal score will lie in [0, Sa] Our Our
St,)
algorithm is summarized in Algorithm[] [ where A and A₂ are the updating parameters associated with A and Ag, respectively, and to denotes the updating frequency. The weight A2 is updated in similar manner to the updating rule on the Lagrange multiplier used by NOTEARS and we set A as an upper bound X2. as previously discussed. In all our experiments that BIC as score function, St is obtained from a complete directed graph and St: is from an empty graph. Since St with the empty graph can be very high for large graphs. we also update it by keeping track of the lowest score achieved by DAGs generated during training. Other parameter choices in this work are Sn = 5. to 1, 000, A = O, A I, A2 10 [4/33/3], A2 = 10 and A2 = 0.01. We comment that these parameter choices may be further tuned for specific applications, and the inferred causal graph would be the one that is acyelic and achieves the best score, among all the final outputs (cf. Section 5.3) of the RL approach with different parameter choices.
5.2 CRROICCCIIIC ALGORITHM
We believe that the exploitation and exploration scheme in the RL paradigm provides appropriate way to guide the search. Let (- s) and y denote the policy and NN parameters for graph generation, respectively. Our training objective is the expected reward defined
J( s) = E E x) (S(g) + I & DAGs) + X2h(A))}
During training, the input s is constructed by randomly drawing samples from the observed dataset X, as described in Section 4
We resort to policy gradient methods and stochastic methods to optimize the parameters vb The gradient V y (ib s) can be obtained by the well-known REINFORCE algorithm (Williams) 1992. Sutton et al.) 2000). We draw B samples $1, ....... Sg as batch to estimate the gradrent which is then used to train the NNs through stochastic optimization methods like Adam (Kingma and Ba 2014). Using parametric baseline to estimate the reward can also help training (Konda and Tsitsiklis
2000 For the present work, our critic is a simple 2-layer feed-forward NN with ReI.U units, with the The critic is trained with Adam between its predictions and the true rewards. An entropy regularization term (Williams and Peng al. 2016) is also added encourage exploration of the RL agent. Although policy gradient methods only guarantee local convergence under proper conditions (Sutton et al., (2000). we remark that the inferred graphs from the actor-crine critic algorithm are all DAGs in our experiments.
input being the encoder outputs {enc;
mean squared error
1991, Mnih Mnih
Training an RL agent typically requires many iterations. In the present work, we find that computing the rewards for generated graphs is much more time-consuming than training NNs. Therefore, we record the computed rewards corresponding to different graph structures. Moreover, the BIC score can be decomposed according to single causal relationships and we also record the corresponding RSS; to avoid repeated computations.
When Su St 0, then we have obtained the solution if we know the graph that achieves Sy or St, otherwise we may simply pick new upper bound Su + 1.

---

Published as a conference paper at ICLR 2020
5.3
FINAL OUTPUT
Since we are concerned with finding a DAG with the best score rather than policy, we record all the graphs generated during the training process and output the one with the best reward. In practice, the graph may contain spurious edges and further processing is needed.
To this end, we can prune the estimated edges in a greedy way, according either the regression performance or the score function. For an inferred causal relationship, we remove a parental variable and calculate the performance of the resulting graph, with all other causal relationships unchanged. If the performance does not degrade or degrade within predefined tolerance, we accept pruning and continue this process with the pruned causal relationship. For linear models, pruning can be simply done by thresholding the estimated coefficients.
Related to the above pruning process is to add to the score function an increased penalty weight on the number of edges of graph. However, this weight is not easy to choose, incur missing edges. In this work, we stick to the penalty weight log m that is included in the BIC score and then apply pruning to the inferred graph in order to reduce false discoveries.
large weight may
EXPERIMENTAL RESULTS
We report empirical results on synthetic and real datasets to compare our approach against both traditional and recent gradient based approaches, including GES (with BIC score) (Chickering 2002
Ramsey et al. 2017), ). the PC algorithm (with Fisher-z test and p-value 0.01) (Spirtes et al. 2000).
ICA-LINGAM (Shimizu et al. 2006), the Causal Additive Model (CAM) based algorithm proposed by Bühlmann et al. (2014), NOTEARS (Zheng et al. 2018). DAG-GNN (Yu et al. 2019). and GraN- DAG (Lachapelle al.) 2019), among others. All these algorithms have available implementationss and we give a brief description on these algorithms and their implementations in Appendix D Default hyper-parameters of these implementations are used unless otherwise stated. For pruning, we use the same thresholding method for ICA-LiNGAM, NOTEARS, and DAG-GNN. Since the authors of CAM and GraN-DAG propose apply significance testing of covariates based on generalized additive models and then declare significance if the reported p-values are lower than or equal to 0.001, we stick to the same pruning method for CAM and GraN-DAG.
The proposed RL based approach is implemented based on existing Tensorflow (Abadi et al. 2016) implementation of neural combinatorial optimizer (see Appendix OppendixD]forr more details). The decoder is modified as described in Section Jand the RL algorithm related hyper-parameters are left unchanged. We pick B- = 64 as batch size at each iteration and dh = 16 as the hidden dimension with the single layer decoder. Our approach is combined with the BIC scores under Gaussianity assumption given in Eqs. (2) and 3 and are denoted as RL-BIC and RL-BIC2, respectively.
We evaluate the estimated graphs using three metrics: False Discovery Rate (FDR), True Positive Rate (TPR), and Structural Hamming Distance (SHD) which is the smallest number of edge additions, deletions, and reversals to convert the estimated graph into the true DAG. The SHD takes into account both false positives and false negatives and a lower SHD indicates a better estimate of the causal graph. Since GES and PC may output unoriented edges, follow Zheng et al. (2018) to treat GES and PC favorably by regarding undirected edges as true positives as long as the true graph has a directed edge in place of the undirected edge.
6.1 LINEAR MODEL WITH GAUSSIAN AND CON-GAUSSIAN NOISE
Given number of variables d, we generate d x d upper triangular matrix as the graph binary adjacency matrix, in which the upper entries are sampled independently from Bern(0.5). We assign edge weights independently from Unif (l-2, 0.5) U [0.5, 2)) to obtain a weight matrix W e R and then sample x = W x + n E Idd from both Gaussian and non-Gaussian noise models. The non-Gaussian noise is the same as the one used for ICA-LiNGAM (Shimizu et al. 2006), which generates samples from Gaussian distribution and passes them through power nonlinearity to make them non -ou-Gassiinn We pick unit noise variances in both models and generate m = 5, 000 samples as our datasets. A random permutation of variables is then performed. This data generating procedure is similar to that used by NOTEARS and DAG-GI GNN and the true causal graphs in both cases are known be identifiable (Shimizu al. 2006
Peters and Bühlmann 2013).
7

---

Published as a conference paper at ICLR 2020
average reward per batch
maximum reward batch
maximum reward
megalte swopae
2200
5600
75l 1930 125ce 19630 17d00 20302
Beration
101
16a
16t
2300 3000
7500 10330 12500 18030 21500 20303 Beration
2500
3000
7500
30000
17500 20000 32500
15933
iteraioo
(a) penalty weights
(b) negative reward
Figure 2: Learning process of the proposed method RL-BIC2 on
dataset.
Table 1: Empirical results on LiNGAM and
data models with 12-node graphs
RL-BIC
RL-BIC2
PC
GES
ICA-LINGAM
CAM
NOTEARS
DAG-GNN
GraN-DAG
FDR 0.28 0.11
0.1:0
0.06 1: 0.04
0.621
01.00
0.67 1:
0.04..0.00
0.11 1.00
0.10
LiNGAM
TPR
0.11
14:0
0.254:0.0
0.04
14:0
0.09
0.95
0.94::0.04
0.37
SHD
7,54
0=0
31.8
3298
0=0
40.4
5.00 2,4
1.4
36.0
FDR 0.38 0.11
0:::
0.52
0.06
0.00
0.70
0.02
0.10
0.70
Linear-
TPR
0.12
1:0
0.31 +0.0
0.24
0.73
0.44
0.90
0.95
0.17
Gaussian
SHD 22.2
0=0
2.06
3.42
2.72
40.8
4,40 0.0
2.06
6.68
We first consider graphs with d 12 nodes. We use n = 64 for constructing the input sample and set the maximum number of iterations to 20, 000. We use threshold 0.3. same as NOTEARS and DAG-GNN with this data model, to prune the estimated edges. Figure 2 shows the learning process of the proposed method RL-BIC2 on 683, 784 different graphs during training, much lower than the total number (around 5.22 x 1026) of DAGs. The pruned DAG turns out to be exactly the same as the underlying causal graph.
dataset. In this example, RL-BIC2 generates
We report the empirical results on LiNGAM and PC and GES perform poorly, possibly because we consider relatively dense graphs for our data generating procedure. CAM does not perform well either, as it assumes nonlinear causal relationships. ICA-LiNGAM recovers all the true causal graphs for LiNGAM data but performs poorly on linear- Gaussian data. This is not surprising because ICA-LiNGAM works for non- non-Gaussian noise and does not provide guarantee for
data models in Table D Both
Gaussian datasets. Both NOTEARS and DAG-GNN have good causal discovery results whereas GraN-DAG performs much worse. We believe that it is because GraN-DAG uses 2-layer feed-forward NNs to model the causal relationships, which may not be able to learn good linear relationship in this experiment. Modifying the feed-forward NNs to linear functions reduces to NOTEARS with negative log-likelihood as loss function, which yields similar performance on these datasets (see Appendix H.T for detailed results). As our proposed methods, we observe that RL-BIC2 recovers all the true causal graphs on both data models this experiment while RL-BIC has worse performance. One may wonder whether this observation is due to the same noise variances that are used in our data models; we conduct additional experiments where the noise variances are randomly sampled and RL-BIC2 still outperforms RL-BIC by a large margin (see also Appendix I.I.p Nevertheless, with the same BIC score, RL-BIC performs much better than GES on both datasets, indicating that the RL approach brings in a greatly improved search ability.
Finally, we test the proposed method on larger graphs with d = 30 nodes, where the upper entries are sampled independently from (0.2...2). This edge probability choice corresponds to the fact that large graphs usually have low edge degrees in practice; see, e.g., the experiment settings of Zheng et al. (2018); Yu et al. 2019);| Lachapelle et al. (201)). To incorporate this prior information in our approach, we add to each gıj common bias term initialized to - 10 (see Appendix E. for details). Considering the much increased search space, we also choose a larger number of observed samples,
128, to construct the input for graph generator and increase the training iterations 40, 000. On LiNGAM datasets, RL-BIC2 has FDR, TPR, and SHD being 0.14 + 0.15, 0.94 = 0.07, and 19.8 + 23.0, respectively, comparable to NOTEARS with O. + 0.09, 0.94 + 0.04, and 17.2 + 12.
8

---

Published as a conference paper at ICLR 2020
6.2 NONLINEAR MODEL WITH QUADRATIC FUNCTIONS
We now consider nonlinear causal relationships with quadratic functions. We generate an upper triangular matrix in similar way to the first experiment. For causal relationship with parents Xpa(t) = [K "'g....]' at the i-th node, we expand Xpa(s) to contain both first- and second-order features. The coefficient for each term is then either 0 or sampled from Unif -0.5, U [0.5, 1]), with equal probability. If a parent variable does not appear in any feature term with a non-zero coefficient, then we remove the corresponding edge in the causal graph. The rest follows the same as in first experiment and here we use the non-Gaussian noise model with 10-node graphs and 5, 000 samples. The true causal graph is identifiable according to PPeters et al. (0144). For this quadratic model, there may exist very large variable values which cause computation issues for quadratic regression. We treat these samples as outliers and detailed processing is given in Appendix E.2
We use quadratic regression for a given causal relationship and calculate the BIC score (assuming equal noise variances) in Eq. (3). For pruning, we simply apply thresholding, with threshold as 0.3. to the estimated coefficients of both first- and second-order terms. If the coefficient of second-order term, e.g., 2i, Vig, is non-zero after thresholding, then we have two directed edges that are from zi, to I; and from Tig to F,, respectively. We do not consider PC and GES in this experiment due to their poor performance in the first experiment. Our results with 10-node graphs are reported in Table 2 which shows that RL-BIC2 achieves the best performance.
Table 2: Empirical results on nonlinear models with quadratic functions.
RI.BBIC2
NOTEARS
NOTEARS-2
NOTEARS.3
ICA-LINGAM
CAM
DAG-GNN
GraN-DAG
FDR
0.04
0,35 0.06
0.10
0
0.32 0.07
0.39
0.10
TPR
0.98
0.04
0,71 = 0.16
1.115
0.79 0:20
0,76 276:=000 0.09
0.78 2: 0.05
0.73 0.55 0,11
SHDD
0.6 1.20
148 3.37
3.82
5.19
20.4 5.00
1.112
28.
3.86
For fair comparison, we apply the same quadratic regression based pruning method to the outputs of NOTEARS, denoted NOTEARS-2. We see that this pruning further reduces FDR, i.e., removes spurious edges, with little effect on TPR. Since pruning does not help discover additional positive edges or increase TPR, we will not apply this pruning method to other methods as their TPRs are much lower than that of RL-BIC2. Finally, with prior knowledge that the function form is quadratic, we can modify NOTEARS apply quadratic functions to modeling the causal relationships, with an equivalent weighted adjacency matrix constructed using the coefficients of the first- and second-order terms, similar to the idea used by GraN-DAG (detailed derivations are given Appendix [F.2). The problem then becomes nonconvex optimization problem with (d - 1 1)d² /2 parameters (which are the coefficients of both first- and second-ordes features), compared to the original NOTEARS with df parameters. This method corresponds to NOTEARS-3 in Table 2 Despite the fact that NOTEARS-3 did not achieve a better overall performance than RL-BIC2, we comment that it discovered almost correct causal graphs (with SHD < 2) more than half of the datasets, but performed poorly on the rest datasets. We believe that it is due to the increased number of optimization parameters and the more complicated equivalent adjacency matrix which make the optimization problem harder solve. Meanwhile, we do not exclude that NOTEARS-3 can achieve a better causal discovery performance with other optimization algorithms.
6.3 NONLINEAR MODEL WITH GAUSSIAN PROCESSES
Given a randomly generated causal graph, we consider another nonlinear model where each causal relationship f; a function sampled from Gaussian process with RBF kernel of bandwidth one. The additive noise n, is normally distributed with variance sampled uniformly. This setting is known to be identifiable according to Peters et al. (2014). We use a setup that also considered by GraN-DAG (Lachapelle al. 2019): 10-node and 40-edge graphs with 1, 000 generated samples.
The empirical results are reported in Table 3 One can see that ICA-LiNGAM, NOTEARS, and model this type of causal relationship. More importantly, these methods operate on a notion of weighted adjacency matrix, which is not obvious here. For our method, we apply Gaussian Process Regression (GPR) with RBF kernel to model the causal relationships. Notice that even though the observed data from function sampled from Gaussian process, it is not guaranteed that GPR
DAG-GNN perform poorly on this data model. A possible reason is that they may not be able to

---

Published as a conference paper at ICLR 2020
with the same kernel can achieve good performance. Indeed, using a fixed kernel bandwidth would lead to severe overfitting that incurs many spurious edges and the graph with the highest reward is usually not DAG. To proceed, normalize the observed data and apply median heuristics for kernel bandwidth. Both our methods perform reasonably well, with RI--BIC outperforming all the other methods.
Table 3: Empirical results on nonlinear models with Gaussian processes.
RL-BIC
RL-BIC2
ICA-LINGAM
NOTEARS
DAG-GNN
GraN-DAG
CAM
FDR
0.14 1: 0.03
0.17 0.12
0.48 0.004
0.48 1 0.19
0.36 0.11
0.12 0.08
0.15 1 0.07
TPR
0.96 0.03
0.09
0.63
0.18 +
0.07 0.03
0.11 0.05
0.02
SHD
6.2 1 1.33
12.0 1: 5.18
30.4 2.550
33.8 1 2.56
34.6 1 1.36
10,2 2.39
10.2 2.93
6.4 REAL DATA
We consider a real dataset to discover a protein signaling network based on expression levels of proteins and phospholipids (Sachs et al. 2005). This dataset is common benchmark in graphical models, with experimental annotations well accepted by the biological community. Both observational and interventional data are contained in this dataset. Since we interested in using observational data to infer causal mechanisms, we only consider the observational data with m = 853 samples. The ground truth causal graph given by Sachs al. (2005) has 11 nodes and 17 edges.
Notice that the true graph is indeed sparse and an empty graph can have an SHD as low as 17. Therefore, we report more detailed results regarding the estimated graph: number of total edges, number of correct edges, and the SHD. Both PC and GES output too many unoriented edges, and we will not report their results here. We apply GPR with RBF kernel to modeling the causal relationships, with the same data normalization and median heuristics for kernel bandwidth as in Section 6.3 We also use CAM pruning on the inferred graph from the training process. The empirical results are given in Table 4 Both RL-BIC and RL-BIC2 achieve promising results, compared with other methods.
Table 4: Empirical results Sachs dataset.
RL-BIC
RL-BIC2 ICA-LINGAM
CAM
NOTEARS
DAG-GNN GraN-DAG
Total Edges
10
10
8
10
20
15
10
Correct Edges
7
4
6
6
6
5
SHD
11
19
13
7 CONCLUDING REMARKS AND FUTURE WORKS
We have proposed to use RL to search for the DAG with the optimal score. Our reward is designed incorporate predefined score function and two penalty terms to enforce acyelicity. We the actor-critic algorithm as our RL algorithm, where the actor is constructed based on recently developed encoder-decoder models. Experiments are conducted on both synthetic and real datasets to show the advantages of our method over other causal discovery methods.
We have also shown the effectiveness of the proposed method with 30-node graphs, yet dealing with large graphs (with more than 50 nodes) is still challenging. Nevertheless, many real applications, like Sachs dataset (Sachs et al. 2005), have a relatively small number of variables. Furthermore, it is possible to decompose large causal discovery problems into smaller ones; see, e.g.. Ma et al. (2008). Prior knowledge or constraint-basee methods is also applicable to reduce the search space.
There are several future directions from the present work. In our current implementation, computing scores is much more time consuming than training NNs. We believe that developing a more efficient and effective score function will further improve the proposed approach. Other powerful RL algo- rithms may also be used. For example, the asynchronous ad advantage actor-critic algorithm has been shown to be effective in many applications (Mmih al.) 2016:) Zoph and Le 2017). 7). In addition, observe that the total iteration numbers used in our experiments are usually more than needed (see, e.g., Figure 2(b).. A proper early stopping criterion will be favored.
10

---

Published as a conference paper at ICLR 2020
ACKNOWLEDGMENTS
The authors are grateful to the anonymous reviewers for valuable comments and suggestions. The authors would also like to thank Prof. Jiji Zhang from Lingnan University, Dr. Yue Liu from Huawei Noah's Ark Lab, and Zhuangyan Fang from Peking University for many helpful discussions.
REFERENCES
M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, Dean, M. Devin, S. Ghemawat, G. Irving. M. Isard, et al. Tensorflow: A system for large-scale scale machine learning. In 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI), 2016.
Bello, H. Pham, Q. V. Le, M. Norouzi, and S. Bengio. Neural combinatorial optimization with reinforcement learning. arXiv preprint arXiv:1611.09940 2016.
K. A. Bollen. Structural Equations with Latent Variables. Wiley, 1989.
Bühlmann, J. Peters, J. Ernest, et al. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 22(6):256-2256 2014.
Y. Chen, M. W. Hoffman, S. G. Colmenarejo, M. Denil, T. P. Lillicrap, M. Botvinick, and N. de Freitas. Learning learn without gradient descent by gradient descent. In International Conference on Machine Learning. 2017.
D. M. Chickering. Learning Bayesian networks is NP-complete. In Learning from Data, pages 121-130. Springer, 1996.
D. M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507-555) 2002.
D. M. Chickering, D. Heckerman, and C. Meek. Large-sample learning of Bayesian networks is NP-hard. Journal of Machine Learning Research, ((ct)):287-13300 2004.
H. Dai, E. Khalil, Y. Zhang, B. Dilkina, and L. Song. Learning combinatorial optimization algorithms over graphs. In Advances Neural Information Processing Systems, 2017.
D. Geiger and D. Heckerman. Learning Gaussian networks. In Conference on Uncertainty in Artificial Intelligence, 1994.
O. Goudet, D. Kalainathan, P. Caillou, I. Guyon, D. Lopez-Paz, and M. Sebag. Learning functional causal models with generative neural networks. In Explainable and Interpretable Models in Computer Vision and Machine Learning, pages 39-80. Springer, 2018
W. Han, Chen, M.-S. Cheon, and H. Zhong. Estimation directed acyclic graphs through two- stage adaptive LASSO for gene network inference. Journal of the American Statistical Association, 111(515):1004
2016.
D. M. Haughton et al. On the choice of a model to fit data from an exponential family. The Annals of Statistics, 16(1):342-355, 1988.
P. O. Hoyer. D. Janzing, M. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems 21, 2009.
B. Huang, K. Zhang, Y. Lin, B. Schölkopf, and C. Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. 2018.
A. Hyttinen, F. Eberhardt, and M. Järvisalo. Constraint-based causal discovery: Conflict resolution with answer set programming. In Conference on Uncertainty in Artificial Intelligence, 2014.
D. Kalainathan, O. Goudet, I. Guyon, D. Lopez-Paz, and M. Sebag. Structural agnostic modeling: Adversarial learning of causal graphs. arXiv preprint arXiv: (t0.://9330499 2018.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. ICI.R, 2014.
11

---

Published as a conference paper at ICLR 2020
V. R. Konda and J. N. Tsitsiklis. Actor-critic algorithms. In Advances in Neural Information Processing Systems, 2000.
S. Lachapelle, P. Brouillard, T. Deleu, and S. Lacoste-Julien. Gradient-based neural DAG learning. arXiv preprini (rXiv:9906002220 2019.
Z. Ma, X. Xie, and Z. Geng. Structural learning of chain graphs via decomposition. Journal of Machine Learning Research, (Dec)):2447-2880, 2008.
V. Mnih, A. P. Badia, M. Mirza, Graves, T. Lillicrap. T. Harley, D. Silver, and K. Kavukcuoglu. Asynchronous methods for deep reinforcement learning. In ICMI, 2016.
Nandy, A. Hauser, M. Maathuis, et al. High-dimensional consistency in score-based and hybrid structure learning. The Annals of Statistics, 4666A):3151-3118 2018.
R. Opgen- Rhein and K. Strimmer. From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data. BMC Systems Biology, 1(1)37, 2007.
J. Pearl. Causality, Cambridge University Press, 2009.
J. Peters and P. Bühlmann. Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1):219-22) 2013.
J. Peters, J. M. Mooij, D. Janzing. and B. Schölkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, (5(1):2009-2053 2014.
J. Peters, D. Janzing, and B. Schölkopf. Elements of Causal Inference Foundations and Learning Algorithms. The MIT Press, Cambridge, MA, 2017.
J. Ramsey, M. Glymour, R. Sanchez-Romero, and C. Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics, 3(2):121-129 2017.
K. Sachs, O. Perez, D. Pe'er, A. Lauffenburger, and G. P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721):523-529, 2005.
G. Schwarz. Estimating the dimension model. The Annals of Statistics, 6(2):461 464, 1978.
S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. A linear non-Gaussian acyelic model for causal discovery. Journal of Machine Learning Research, 7(0ct):2003-2030, 2006.
S. Shimizu, T. Inazumi, Y. Sogawa, A. Hyvärinen, Y. Kawahara, T. Washio, P. O. Hoyer, and K. Bollen. Directlingam: direct method for learning a linear non non-Gaussian structural equation model. Journal of Machine Learning Research, 12(Apr):125512258 2011.
D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, et al. Mastering the game of go without human knowledge. Nature, 550(7676): 354, 2017.
R. Socher, Chen, C. D. Manning, and A. Ng. Reasoning with neural tensor networks for knowledge base completion. In Advances in Neural Information Processing Systems, pages 226-934, 2013.
Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT press, Cambridge, MA, USA, 2nd edition, 2000.
X. Sun, D. Janzing, B. Schölkopf, and K. Fukumizu. A kernel-based causal learning algorithm. In International Conference Machine Learning, 2007.
Sutskever, O. Vinyals, and Q. V. Le. Sequence to sequence learning with neural networks. In Advances in Neural Information Processing Systems, 2014.
12

---

Published as a conference paper at ICLR 2020
R. S. Sutton, D. A. McAllester, S. Singh, and Y. Mansour. Policy gradient methods for reinforcement learning with function approximation. In Advances in Neural Information Processing Systems, 2000.
I. Tsamardinos, L. E. Brown, and C. F. Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, k.. Kaiser, and I. Polosukhin. Attention all you need. Advances in Neural Information Processing Systems, 2017.
Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio. Graph attention networks. International Conference Learning Representations, 2018.
O. Vinyals, M. Fortunato, and N. Jaitly. Pointer networks. In Advances in Neural Information Processing Systems, 2015.
R. J. Williams. Simple statistical tradient-following algorithms for connectionist reinforcement learning. Machine Learming, 8(3-4):229-256, 1992.
R. Williams and J. Peng. Function optimization using connectionist reinforcement learning algorithms. Connection Science, 3(3):241-26,, 1991.
Y. Yu, J. Chen, T. Gao, and M. Yu. DAG-GNN: DAG structure learning with graph neural networks. In ICML, 2019.
K. Zhang and A. Hyvärinen. On the identifiability of the post-nonlinear causal model. In Conference on Uncertainty in Artificial Intelligence, 2009.
K. Zhang, J. Peters, D. Janzing, and B. Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Conference Uncertainty in Artificial Intelligence, 2012.
X. Zheng. B. Aragam, P. Ravikumar, and E. P. Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, 2018.
B. Zoph and Q V. Le. Neural architecture search with reinforcement learning. In ICLR, 2017.
13

---

Published as a conference paper at ICLR 2020
APPENDIX
A MORE DETAILS ABOUT DECODERS
We briefly describe the NN based decoders for generating binary adjacency matrices:
Single layer decoder:
9ig 1 W2, u) = u² tanh(Wi enc; + W2 ency),
where Wi, W2 Rd, Rd are trainable parameters and da denotes the hidden dimension associated with the decoder.
Bilinear decoder:
Bey (W) ence Wennenn Wennj
where W& E Rd a trainable parameter.
Neural Tensor Network (NTN) decoder (Socher al. 2013):
Dij (W((k) V,b) u tanh (enc,
Wil: (enc]TW(l:k ene, + , + +b).
where W[ E Tidde x is a tensor and the bilinear tensor product enc? wyli:klenc results a vector with each entry being ency W Wiklene; for 1,2,..., K, V Kx2d. u E Kx1 and K
Transformer decoder uses the multi-head attention module obtain the decoder outputs {dec, followed by ferd-forward NN whose weights are shared across all dec; (Vaswani et al. 2017). The output consists of vectors in R2 which are treated as the row vectors of d x matrix. We then pass each element of this matrix into sigmoid functions and sample a binary adjacency matrix accordingly.
Table provides the empirical results on linear-Gaussian Gaussian data models with 1 12-node graphs and unit variances (see Section 6.1 for more details on this data generating procedure). In our implementation, pick de 64 as the dimension of the encoder output, d,. 16 for the single layer decoder, and K for the NTN decoder. We find that single layer decoder performs the best, possibly because has less parameters and easier train find better DAGs, while the Transformer encoder has provided sufficient interactions amongst variables.
Table 5: Empirical results of different decoders on Ginssr-Gaussian data models with 12-node graphs.
Single layer
Bilinear
NTN
Transformer
FDR
0=0
0.07
0.12
0.07
Linear- Gaussian
TPR
1
0.09 0.95 + 0.08
0.11
SHD
0+:0
4.0 7.04
7.52
38.2
B EQUIVALENCE OF PROBLEMS AND
66
Proof of Proposition Let G be solution to Problem 1 Then we have S = S(G) and & must be a DAG due to the hard acyclicity constraint. Assume that G is not solution to Problem (6). which indicates that there exists directed graph G' (with binary adjacency matrix A') so that
S S(g)) I( & DAGs) + >2h(A').
(8)
Clearly, 9' cannot be DAG, for otherwise we would have a DAG that achieves lower score than the minimum S*. By our assumption, it follows that
r.h.s. of Eq. Sr. + Ai + Agtimin 2 Su,
which contradicts the fact that Sy is upper bound on S*.
14

---

Published as a conference paper at ICLR 2020
For the other direction, let g be solution to Problem (6) but not a solution to Problem T.. This indicates that either 9 is not a DAG or G is DAG but has higher score than the minimum score, i.c. (()) > S* The latter case clearly contradicts the definition of the minimum score. For the former case, assume that some DAG g' achieves the minimum score. Then plugging g' into the negative reward, we can get the same inequality in Eq. (8) since both penalty terms zeros for a DAG. This then contradicts the assumption that G minimizes the negative reward.
C PENALTY WEIGHT CHOICE
Although setting > = 0, or equivalently using only the indicator function w.r.t. acyclicity, can still make Problem (6) equivalent the original problem with hard acyclicity constraint, we remark that this choice usually does not result in good performance of the RL approach, largely due to that the reward with only the indicator term is likely to fail guide the RL agent to generate DAGS.
To see why it is the case, consider two cyclic directed graphs, one with all the possible directed edges in place and the other with only two edges (i.e., 4; + 2j and zy + 24 for some i * j). The latter is much 'closer' to acyclicity in many senses, such as h(.A) given in Eq. (4) and number of edge deletion, addition, and reversal make a directed graph acyclic. Assume linear data model that has relatively dense graph. Then the former graph will have lower BIC score when using linear regression for fitting causal relations, yet the penalty terms of acyclicity are the same with only the indicator function. The former graph then has higher reward, which does not help the agent to tend to generate DAGs, Also notice that the graphs in our approach are generated according to Bernoulli distributions determined by NN parameters that are randomly initialized. Without loss of generality, consider that each edge is drawn independently according to (ern(0.5). For small graphs (with less than equal to nodes so), few hundreds of samples of directed graphs are very likely to contain a DAG. Yet for large graphs, the probability of sampling a DAG is much lower. If no DAG generated during training, the RL agent can hardly learn to generate DAGs. The above facts indeed motivate us to choose small value of A2 so that the agent can be trained to produce graphs closer to acyclicity and finally to generate exact DAGs.
A question is then what if the RL approach starts with a DAG, e.g., by initializing the probability of generating each edge to be nearly zero. This setting did not lead good performance, either. The generated directed graphs at early iterations can be very different from the true graphs in that many true edges do not exist, and the resulting score is much higher than the minimum under the DAG constraint. With small penalty weights of the acyclicity terms, the agent could be trained to produce cyclic graphs with better scores, similar to the case with randomly initialized NN parameters. On the other hand, large initial penalty weights, as we have discussed in the paper, limit exploration of the RL agent and usually result in DAGs whose scores are far from optimum.
D IMPLEMENTATION DETAILS
We use existing implementations of causal discovery algorithms in comparison, listed below:
ICA-LINGAM (Shimizu et al. 2006) assumes linear non-Gaussian additive model for data generating procedure and applies Independent Component Analysis (ICA) to recover the weighted adjacency matrix, followed by thresholding on the weights before outputting the inferred graph. A Python implementation is available at the first author's website https://ss http///stessseee google. com/ site/eete /esee//sssiiiiiim wwt.ssswwi...........oo. zu0 / nngam
GES and PC: we use the fast greedy search implementation of GES (Ramsey et al. 2017) which has been reported to outperform other techniques such as max-m min hill climbing (Han al. 2016 Zheng et al. 2018). Implementations of both methods are available through the py-sausal package at wtthtpss: / thub. com/ wwtpcccmccc///eeeeecc http/////cccccccuuu/eeepaa causal written in
highly optimized Java codes.
CAM (Peters et al. 2014) decouples the causal order search among the variables from feature or edge selection in DAG. CAM also assumes additive noise as in work, with an additional condition that each function is nonlinear. Codes are available through the CRAN R package repository at httttt:://rrrte r-project org/ wttp///bbpaaaaa weeppackages//Cr
15

---

Published as a conference paper at ICLR 2020
NOTEARS (Zheng et al. 2018) recovers the causal graph by estimating the weighted adjacency matrix with the least squares loss and the smooth characterization for acyclicity constraint, followed by thresholding on the estimated weights. Codes are available at the first author's github repository https:/ / / t.hb..... com/ xunzheng loinng/nottaaa /notears We also re-implement the augmented Lagrangian method following the same updating rule the Lagrange multiplier and the penalty parameter in Tensorflow, so that the augmented Lagrangian each iteration can be readily minimized without the need of obtaining closed- form gradients. We use this implementation in Sections 6.2 and 6.3 when the objective function and/or the acyclicity constraint are modified.
DAG-GNN (Yu al., 2019) formulates causal discovery in the framework of variational autoencoder, where the encoder and decoder are two shallow graph NNs. With modified smooth characterization on acyclicity. DAG-GNN optimizes a weighted adjacency matrix with the evidence lower bound as loss function. Python codes are available at the first http://gtttt.u.g.o comm shmoon 234 2341124--0A2200
author's github repository
GraN-DAG (Lachapelle et al. and chooses the sum of all product paths between variables I; and Φy as the (i,j)-th element of an equivalent weighted adjacency matrix. GraN-DAG uses the same smooth constraint from Zheng et al. (2018) to find DAG that maximizes the oggllikelhhod of the observed samples. Codes are available at the first author's github repository https DAG
2019) uses feed-foward NNs to model each causal relationship
7/github cww http://koowwsaannaae
Our implementation is based on an existing Tensorflow implementation of neural com- binatorial optimizer that is available https: ponozall-ombinttorllll--ooo weiiiiiiiiiiniai wow.oo------eesoooooooooooo We add an entropy regularization term, and modify the reward and decoder as described in Sections and 5.1 respectively.
E FURTHER EXPERIMENT DETAILS AND RESULTS
E.1 EXPERIMENT IN SECTION 6.1
We replace the feed-forward NNs with linear functions in GraN-DAG and obtain similar experimental results NOTEARS (FDR, TPR, SHD): 0.05 + 0.04, 0.93 = 0.06, 3.2 + 2.93 and 0.05 + 0.04, 0.95 + 0.03, 2.40 + 1.85 for LiNGAM and
data models, respectively.
We conduct additional experiments with linear models where the noise variances uniformly sampled according to Unif([0.5, 2)). Results are given in Table
Table 6: Empirical results on LiNGAM and Gaussian data models with 12-node graphs and
different noise variances.
RL-BIC
RL-BIC2
PC
GES
ICA-LINGAM
CAM
NOTEARS
DAG-GNN
GraN-DAG
FDR 1:
0.09 0.06
0.57 0.10
0.13
0::0
0.70
0.8 10.00
0.14 0.07 0.71 0.110
LINGAM TPR 0.77 0.15
0.94
0,28=006
0.11
1=0
0.12
0.94
0.04
0.25 0.09
SHD 14.4 = 7.17
2.0
3.14
3200
0:0
41.6
3.99
6.6 1.02
3.87
Linear- FDR 4:
0.10:1.077
0.54 0.10
0.61 40.14
0.67
0.05
0.60
0.77
0.09
0.12
0.112
TPR 0.68 = 0.09
0.93
0.29
0.21
0,75
0.11
0.90
0.94 0.04
0.21 0.08
Gaussian
SHD 18.8 - 3.43
3.07
30.0 2.83 3:42
49.0
37.8
3.0
5:4 2.06
39.6
Knowing a sparse true causal graph a priori is also helpful. To incorporate this information in our experiment with 30-node graphs, we add additional biased term E IR to each decoder output: for the single layer decoder, we have
9ry W2, u) = nT tanh (W1 enc, + W2 ency) +
where we let è be trainable and other parameters have been defined in Appendix A In our experiments, è is initialized to be -00; this choice aims to set a good starting point for generating graph adjacency matrices, motivated by the fact that a good starting point is usually helpful to locally convergent algorithms.
16

---

Published as a conference paper at ICLR 2020
E.2 EXPERIMENT IN SECTION 6.2
To remove 'outlier' samples with large values that may cause computation issues for quadratic regression, we sort the samples in ascending order according their &2-norms and then pick the first 3, 000 samples for causal discovery.
We use similar idea from GraN-DAG to build an equivalent weighted adjacency matrix for NOTEARS. Take the first variable for example. We first expand the rest variables ......, contain both first- and second-order features: I2,,... Bdy 22 ... ".. ,d with i,j = ...... d and < j. There are in total d(d - 1)/2 terms and we use to denote the vector that concatenates these feature terms. Correspondingly. we use c; and Cty denote the coeffi- cients associated with these features and с₁ to denote the concatenating vector of the coefficients. Notice that the variable vt, + affects #1 only through the terms 2t, 2;30 with 1, and æyzz with > Therefore, an equivalent weighted adjacency matrix W lying in Radd can constructed with the (1, 1)-th entry Wn = |ai + Ei=2 |cal + Σj=141 lctl: in this way, Wi = 0 implies that Ij has no effect on Z1. The least squares term, corresponding to variable Tı, in the loss function will become Σ (ur
of xt)² where is the total number of samples. In summary, we have the
following optimization problem
min Cg.0o
c c/ x
k=1 subject to trace (e (Wol)) d 0,
where o denotes the element- wise product and the constraint enforces acyclicity w.r.t weighted adjacency matrix (Zheng al. al.:22002001 2018). The problem has (d 1)d² parameters to optimize, while the original NOTEARS optimizes d parameters. We solve this problem with augmented Lagrangian method where at each iteration the augmented Lagrangian approximately minimized by Adam (Kingma and Ba 2014) with Tensorflow. The Lagrange multiplier and the penalty parameter are updated the same fashion as in the original NOTEARS.
17


Reviews:
[{"review_id": "SJguv9wTtB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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 work addresses the task of causal discovery. The proposed contribution is to apply prior work which uses reinforcement learning for combinatorial optimization to structure learning. Specifically, the proposed optimization problem seeks to maximize a penalized score criterion subject to the acyclicity constraint proposed by Zheng, et al. Empirical results show the proposed method performing favorably in contrast to prior art. \n\nOverall I think this is a sensible idea, and the authors do a nice job of exposition, and empirical evaluation. \n\nMy concerns are as follows:\n\n* The novelty is somewhat limited, since the paper is combining two previously proposed ideas (combinatorial search and the acyclicity constraint) for structure learning.\n\n* The paper is loose with technical points. Specifically, the authors claim to use the additive noise model, but then make no restrictions on f(). In this setting, it is fairly well known that we can only hope to learn up to the Markov equivalence class (not the fully directed graph), but there is no mention of this in the paper. \n\nWith all of this said, I think overall the paper is an interesting addition to the causal discovery literature. ", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HklEzhj6FB", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "Update: after the revision, I have decided to increase my score to 8.\n\nOriginal comments:\n\nIn this paper, the authors proposed a new reinforcement learning based algorithm to learn causal graphical models. Simulations on real and synthetic data also shows promise.\n\nPros\n\n1. It's great to see the authors has done a comprehensive comparison with the other methods, especially under different simulation scenarios.\n\n2. The novel idea of applying reinforcement learning to DAG search sounds intriguing. Reinforcement learning offers a powerful tool for policy evaluation and decision making. It’s good to see that the author can successfully extend such toolbox to the field of causal structure learning. To the best of the author’s knowledge, such idea has never been considered by previous work in causal graphical models.\n\nCons. \n\n1. In the introduction section, the authors claimed that “GES is not guaranteed in the finite sample regime”. This seems to be incorrect. For example, the Nandy et al. paper tackles exactly the finite sample problem. \n\nIn conclusion, overall this is a sensible idea, although some of the preliminaries still remain to be polished.", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "N/A"}, "extra_text": {}}, {"review_id": "H1eabJol9B", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "In this paper, the authors propose an RL-based structure searching method for causal discovery. The authors reformulate the score-based causal discovery problem into an RL-format, which includes the reward function re-design, hyper-parameter choose, and graph generation. To my knowledge, it’s the first time that the RL algorithm is applied to causal discovery area for structure searching.\n \nThe authors’ contributions are:\n(1) re-design the reword function which concludes the traditional score function and the acyclic constraint\n\n(2) Theoretically prove that the maximizing the reward function is equivalent to maximizing the original score function under some choices of the hyper-parameters.\n\n(3) Apply the reinforce gradient estimator to search the parameters related to adjacency matrix generation.  \n\n(4) In the experiment, the authors conduct experiment on datasets which includes both linear/non-linear model with Gaussian/Non-gaussian noise.\n\n(5) The authors public their code for reproducibility.\n \nOverall, the idea of this paper is novel, and the experiment is comprehensive. I have the following concerns.\n \n(1) In page 4 Encoder paragraph, the authors mention that the self-attention scheme is capable of finding the causal relationships. Why? In my opinion, the attention scheme only reflects the correlation relationship. The authors should give more clarifications to convince me about their beliefs.\n \n(2) The authors first introduce the h(A) constraint in eqn. (4), and mentioned that only have that constraint would result in a large penalty weight. To solve this, the authors introduce the indicator function constraint. What if we only use the indicator function constraint? In this case, the equivalence is still satisfied, so I am confused about the motivation of imposing the h(A) constraint.\n \n(3) In the last paragraph of page 5, why the authors adjust the predefined scores to a certain range?\n \n(4) Whether the acyclic can be guaranteed after minimizing the negative reward function (the eqn.(6))? I.e., After the training process, whether the graph with the best reward can be theoretically guaranteed to be acyclic?\n \n(5) In section 5.3, the authors mention that the generated graph may contain spurious edges? Whether the edges that in the cyclic are spurious? Whether the last pruning step contains pruning the cyclic path?\n \n \n(6) In the experiment, the authors adopt three metrics. For better comparison, the author should clarify that: the smaller the FDR/SHD is, the better the performance, and the larger the TPR is, the better the performance.\n\n(7) From the experimental results, the proposed method seems more superiors under the non-linear model case. Why? Could the authors give a few sentences about the guidance of the model selection in the real-world? i.e., when to select the proposed RL-based method? And under which case to choose RL-BIC, and which case to selection RL-BIC2?\n \n(8) What’s training time, and how many samples are needed in the training process?\n \n \nMinor:\n1. In the page 4 decoder section, the notation of enc_i and enc_j is not clarified.\n\n2. On page 5, the \\Delta_1 and \\Delta_2 are not explained.\n\n3. For better reading experience, in table 1,2,3,4, the authors should bold value that has the best performance.\n\n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an RL-based structure search method for causal discovery. The reviewers and AC think that the idea of applying reinforcement learning to causal structure discovery is novel and intriguing. While there were initially some concerns regarding presentation of the results, these have been taken care of during the discussion period. The reviewers agree that this is a very good submission, which merits acceptance to ICLR-2020.
Decision: Accept (Talk)