Gradient-Based Neural DAG Learning
Authors: Sébastien Lachapelle
Philippe Brouillard
Tristan Deleu
Simon Lacoste-Julien
Venue:  (2020)
Keywords: Structure Learning
Causality
Density estimation

Abstract: We propose a novel score-based approach to learning a directed acyclic graph (DAG) from observational data. We adapt a recently proposed continuous constrained optimization formulation to allow for nonlinear relationships between variables using neural networks. This extension allows to model complex interactions while avoiding the combinatorial nature of the problem. In addition to comparing our method to existing continuous optimization methods, we provide missing empirical comparisons to nonlinear greedy search methods. On both synthetic and real-world data sets, this new method outperforms current continuous methods on most tasks while being competitive with existing greedy search methods on important metrics for causal inference.

Published as a conference paper at ICLR 2020
GRADIENT-BASED NEURAL DAG LEARNING
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu & Simon Lacoste-Julien'
Mila & DIRO
Université de Montréal
ABSTRACT
We propose novel score-based approach learning directed acyclic graph (DAG) from observational data. We adapt a recently proposed continuous con- strained optimization formulation to allow for nonlinear relationships between variables using neural networks. This extension allows to model complex in- teractions while avoiding the combinatorial nature of the problem. In addition to comparing our method to existing continuous optimization methods, we pro- vide missing empirical comparisons to nonlinear greedy search methods. On both synthetic and real-world data sets, this new method outperforms current continu- ous methods on most tasks, while being competitive with existing greedy search methods important metries for causal inference.
INTRODUCTION
Structure learning and causal inference have many important applications in different areas of sci- ence such as genetics (Koller& & Friedman, 2009; Peters et al., 2017), biology (Sachs et al., 2005) and cconomics (Pearl, 2009). Bayesian networks (BN), which encode conditional independencies using directed acyclic graphs (DAG), are powerful models which are both interpretable and computation- ally tractable. Causal graphical models (CGM) (Peters et al., 2017) are BNs which support inter- ventional queries like: What will happen if someone external to the system intervenes on variable X? Recent work suggests that causality could partially solve challenges faced by current machine learn- ing systems such as robustness to out-of-distribution samples, adaptability and explainability (Pearl, 2019; Magliacane et al., 2018). However, structure and causal learning are daunting tasks due to both the combinatorial nature of the space of structures (the number of DAGs grows super expo- nentially with the number of nodes) and the question of structure identifiability (see Section 2.2). Nevertheless, these graphical models known qualities and promises of improvement for machine intelligence renders the quest for structure/causal learning appealing.
The typical motivation for learning causal graphical model is predict the effect of various in- terventions. A CGM can be best estimated when given interventional data, but interventions are often costly impossible data and rely on different assumptions which make the graph identifiable from the distribution (see Section 2.2). This the approach employed in this paper.
obtained. As an alternative, one can use exclusively observational
We propose score-based method (Koller & Friedman, 2009) for structure learning named GraN- DAG which makes use of a recent reformulation of the original combinatorial problem of finding an optimal DAG into a continuous constrained optimization problem. In the original method named NOTEARS (Zheng et al., 2018), the directed graph is encoded as a weighted adjacency matrix which represents coefficients in linear siruciural equation model (SEM) (Pearl, 2009) (see Section 2.3) and enforces acyclicity using a constraint which is both efficiently computable and easily differen- tiable, thus allowing the use of numerical solvers. This continuous approach improved upon popular methods while avoiding the design of greedy algorithms based on heuristics.
Our first contribution is to extend the framework of Zheng et al. (2018) to deal with nonlinear relationships between variables using neural networks (NN) (Goodfellow et al., 2016). To adapt the acyclicity constraint to our nonlinear model, we argument similar what is used in
*Caaaaa CIFAR AI Chair Correspondence sebast lachape Techapesle@uuoore reau..ntreel...

---

Published as a conference paper at ICLR 2020
Zheng et al. (2018) and apply it first the level of neural network paths and then at the level of graph paths. Although GraN-DAG is general enough to deal with a large variety of parametric fam- ilies of conditional probability distributions, Gaussian additive noise models since, under specific assumptions, provides appealing theoret- ical guarantees easing the comparison to other graph search procedures (see Section 2.2 & 3.3). On both synthetic and real-world tasks, we show GraN-DAG often outperforms other approaches which leverage the continuous paradigm, including DAG-GNN (Yu al., 2019), a recent nonlinear extension of Zheng et al. (2018) which uses an evidence lower bound as score.
experiments focus the special case of nonlinear
Our second contribution is provide a missing empirical comparison to existing methods that support nonlinear relationships but tackle the optimization problem in its discrete form using greedy search procedures, namely CAM (Bühlmann al., 2014) and GSF (Huang et al., 2018). We show that GraN-DAG is competitive on the wide range of tasks we considered, while using pre- and post-processing steps similar to CAM.
We provide implementation of GraN-DAG here.
2
BACKGROUND
Before presenting GraN-DAG, we review concepts relevant to structure and causal learning.
2.1 CAUSAL GRAPHICAL MODELS
We suppose the natural phenomenon of interest can described by a random vector X E Rd entailed by an underlying CGM (Px, g) where Px is probability distribution over X and G (V, E) is a DAG (Peters et al., 2017). Each node j E V corresponds to exactly one variable in the system. Let n. denote the set of parents of node j 9 and let X_e denote the random vector containing the variables corresponding to the parents of j in G. Throughout the paper, we assume there are no hidden variables. In a CGM, the distribution Px is said be Markov G, i.e. we can write the probability density function (pdf) of Px as p(x) = II (, where P5 pz(zy (1y) is the conditional pdf of variable Xj given X-g. A CGM can be thought of as BN in which directed edges are given a causal meaning, allowing it to answer queries regarding interventional distributions (Koller & Friedman, 2009).
2.2 STRUCTURE IDENTIFIABILITY
In general, is impossible to recover 9 given only samples from Px, i.e. without interventional data. is, however, customary rely on a set of assumptions render the structure fully or partially identifiable.
Definition Given a set of assumptions A on a CGM M = (Px.9). its graph G is said to be identifiable from Px if there exists no other CGM M = (Px,9) satisfying all assumptions in A such that & 4 G and Px- = Px-
There are many examples of graph identifiability results for continuous variables (Peters et al., 2014; Peters & Bühlman, 2014; Shimizu al., 2006; Zhang & Hyvärinen, 2009) as well for discrete variables (Peters et al., 2011) 1). These results are obtained by assuming that the conditional densities belong a specific parametric family. For example, if one assumes that the distribution Px is entailed by a structural equation model of the form
x= fX,,)) + N₂ with N, o3) VjeV
where fy is a nonlinear function satisfying some mild regularity conditions and the noises N, are mutually independent, then is identifiable from Px (see Peters et al. (2014) for the complete theorem and its proof). This a particular instance of additive noise models (ANM). We will make of this result in our experiments in Section 4.
One can consider weaker assumptions such as faithfulness (Peters et al., 2017). This assumption allows one to identify, not & itself, but the Markov equivalence class which it belongs (Spirtes
2

---

Published as a conference paper at ICLR 2020
et al., 2000). A Markov equivalence class is set of DAGs which encode exactly the same set of conditional independence statements and can be characterized by a graphical object named a completed partially directed acyclic graph (CPDAG) (Koller & Friedman, 2009; Peters al., 2017). Some al algorithms we use as baselines in Section output only CPDAG.
2.3 NOTEARS: CONTINUOUS OPTIMIZATION FOR STRUCTURE LEARNING
Structure learning is the problem of learning & using a data set of n samples f (x)) from Px. Score-based approaches cast this problem an optimization problem, i.c. ف = arg maXGEDAG S(G) where (()) is regularized maximum likelihood under graph g. Since the number of DAGs is super exponential in the number of nodes, most methods rely on various heuristic greedy search procedures to approximately solve the problem (see Section for review). We now present the work of Zheng et al. (2018) which proposes to cast this combinatorial optimiza- tion problem into continuous constrained one.
To do so. the authors propose to encode the graph & on d nodes as weighted adjacency matrix U [u,] |ua] E paxd which represents (possibly negative) coefficients linear SEM of the form X; : ; X + N₂ Vj where Nj is a noise variable. Let Gu be the directed graph associated with the SEM and let Ay be the (binary) adjacency matrix associated with Gy. One can see that the following equivalence holds;
(Av); = U,
(2)
To make sure Gy is acyclic, the authors propose the following constraint on U:
UU
(3)
where e. 2 E
M*
is the matrix exponential and O is the Hadamard product. ET
To why this constraint characterizes acyclicity, first note that (Au of length passing through node j in graph Gu- Clearly, for Gy to be acyelic, we must have
(Au*)a ")i is the number of cycles
Tr Au 0 for k = 2.., xc. By equivalence (2), this is true when Tr( @ U) = 0 for k 1, ..., 30 From there, one can simply apply the definition of the matrix exponential to why constraint (3) characterizes acyclicity (see Zheng et al. (2018) for the full development).
The authors propose to regularized negative least square score (maximum likelihood for Gaussian noise model). The resulting continuous constrained problem is
1 | XUI 2n
max max SUU, &
A||U/|| s.t. Tre -d =d=0 0
(4)
where X E Rn is the design matrix containing all n samples. The nature of the problem has been drastically changed: went from combinatorial to continuous problem. The difficulties of combinatorial optimization have been replaced by those of non-convex optimization, since the feasible set is non-convex. Nevertheless, standard numerical solver for constrained optimization such has an augmented Lagrangian method (Bertsekas, 1999) can be applied to get an approximate solution, hence there is need design a greedy search procedure. Moreover, this approach is more global than greedy methods in the sense that the whole matrix U is updated at each itera- tion. Continuous approaches to combinatorial optimization have sometimes demonstrated improved performance over discrete approaches in the literature (see for example Alayrac et al. (201 8, $5.2) where they solve the multiple sequence alignment problem with continuous optimization method).
GRAN-DAG: GRADIENT-BASED NEURAL DAG LEARNING
We propose new nonlinear extension to the framework presented in Section 2.3. For each variable Xj. we Iearn a fully connected neural network with hidden layers parametrized by Φ6) A (1)
(i) ..... w((++1) where ()) is the 4th weight matrix of the jth NN (biases are omitted for clarity). Each NN takes as input X_,
v E Rd,, ie. the vector X with the jth component masked to zero, and outputs (g) E Rm the m-dimennsiona parameter vector of the desired distribution family
3

---

Published as a conference paper at ICLR 2020
for variable Xy.' The fully connected NNs have the following form
B & w (5) g.....(W(B 9(W () x-)))...) Vj
(5) where 9 is nonlinearity applied element-wise. Note that the evaluation of all NNs be parallelized GPU. Distribution families need not be the same for each variable. Let ø - (d())...,, ((4)) } represents all parameters of all d NNs. Without any constraint on its param- eter ()), neural network j models the conditional pdf p)(z)l=-z; Note that the product I (= pose according DAG. We now show how one can constrain e to make sure the product of all conditionals outputted by the NNs is joint pdf. The idea is to define a new weighted adjacency ma- trix Ag similar to the one encountered in Section 2.3, which can be directly used inside the constraint of Equation 3 to enforce acyclicity.
does not integrate to one (i.e. it is not joint pdf), since it does not decom-
3.1 NEURAL NETWORK CONNECTIVITY
Before defining the weighted adjacency matrix Ag. we need to focus on how one can make some NN outputs unaffected by some inputs. Since we will discuss properties of single NN, we drop the NN subscript (j) to improve readability.
We will use the term neural network path to refer to computation path in a NN. For example, in w-(3)) h, hah. khz is NN path from input i to output k. We say that a NN path is inactive if at least one weight along the path is zero. We can loosely interpret the path product [W h, (1 ||
NN with two hidden layers. the sequence of weights (W()) w(2)
khg 2 0 as the strength of the NN path, where a path product is equal to zero if and only if the path is inactive. Note that if all NN paths from input i output are inactive (i.e. the sum of their path products is zero), then output & does not depend on input i anymore since the information in input will never reach output k. The sum of all path products from input to output for all input and output it can be easily computed by taking the following matrix product.
[[
hyh,
C A ( (W (I E R 20
(6)
where WW [ is the dementtwwsee absolute value of W. Let us name C the neural network connectivity matrix. It can be verified that Csa is the sum of all NN path products from input i to output k. This means is sufficient to have Ck = 0 to render output independent of input i.
Remember that each NN in our model outputs parameter vector 0 for conditional distribution and that we want the product of all conditionals to be a valid joint pdf, i.c. we want its corresponding directed graph to be acyclic. With this in mind, we see that it could be useful to make a certain parameter 0 not dependent on certain inputs of the NN. To have 0 independent of variable X,, is sufficient to have / Cki O.
3.2 A WEIGHTED ADIACENCY MATRIX
We now define a weighted adjacency matrix As that can be used in constraint of Equation 3.
(Ao) A { Em (Cis) ke
if +
0,
otherwise
where Ciyy denotes the connectivity matrix of the NN associated with variable X5.
As the notation suggests, Ag E R depends on all weights of all NNs. Moreover, it can effectively be interpreted as a weighted adjacency matrix similarly to what we presented in Section 2.3, since we have that
(Ao)iz 0 = (G) does not depend on variable X;
(8)
We note G6 to be the directed graph entailed by parameter Φb We can now write our adapted acyclic- ity constraint:
h(0) Tre d =0
Not all parameter vectors need to have the same dimensionality, but to simplify the notation, suppose Vj
my
4

---

Published as a conference paper at ICLR 2020
Note that can compute the gradient of h()) w.r.t ø (except at points of non-differeniaabiiiy arising from the absolute value function, similar to standard neural networks with ReLU activa- tions (Glorot et al., 2011); these points did not appear problematic in our experiments using SGD).
3.3 A DIFFERENTIABLE SCORE AND ITS OPTIMIZATION
We propose solving the maximum likelihood optimization problem
max Ex~Px Px
log (X,lx 23:00) Σ j=1
s.t. T
(10)
where r denotes the set of parents of node j in graph 9g. Note that Σ3= log (X,1X (()) is
a valid log-likelihood function when constraint (9) is satisfied.
As suggested Zheng al. (2018), we apply an augmented Lagrangian approach to get an approx- imate solution to program (10). Augmented Lagrangian methods consist of optimizing a sequence of subproblems for which the exact solutions are known to converge stationary point of the con- strained problem under some regularity conditions (Bertsekas, 1999). In our case, each subproblem
max [(o, (ilo,,,,HA Ae, Ex~Px (X Ach(o) H M(O)2
(1I)
where À and je the Lagrangian and penalty coefficients of the tth subproblem, respectively. These coefficients are updated after each subproblem is solved. Since GraN-DAG rests on neural networks, we propose to approximately solve each subproblem using well-known stochastic gra- dient algorithm popular for NN in part for its implicit regularizing effect (Poggio et al., 2018), See Appendix A.I for details regarding the optimization procedure.
In the current section, we presented GraN-DAG in a general manner without specifying explicitly which distribution family parameterized by 0g)) In principle, any distribution family could be em- ployed as long as its ligelikeldoood can be computed and differentiated with respect to its parameter However, is not always clear whether the exact solution of problem (10) recovers the ground truth graph G. It will depend on both the modelling choice of GraN-DAG and the underlying CGM (Px.9).
Proposition 1 Let 6* and Gg- be the optimal solution to (10) and its corresponding graph, respec- tively. Let M(A) be the set of CGM (P', ,') for which the assumptions in A are satisfied and let C be the set of CGM (P', g') which can be represented by the model (e.g. NN outputting a Gaussian distribution). If the underlying CGM (Px, 9) E C and C = MM(A for a specific set of assumptions A such that G is identifiable from Px, then G g.
Proof: Let Po be the joint distribution entailed by parameter Φ. Note that the population log- likelihood Ex-Px log (XX) is maximal iff Po Px. We know this maximum can be achieved by specific parameter ø* since by hypothesis (Px,9) E C. Since G is identifiable from Px, we know there exists no other CGM (P) E C such that ② A G and Px = Px Hence Gø- has to be equal to G.
In Section 4.1, we empirically explore the identifiable setting of nonlinear Gaussian ANMs intro- duced in Section 2.2. practice, one should keep mind that solving (10) exactly is hard since the problem is non-convex (the augmented Lagrangian converges only to stationary point) and more- over we only have access to the empirical log-likelihood (Proposition holds for the population case).
3.4 THRESHOLDING
The solution outputted by the augmented Lagrangian will satisfy the constraint only up to numerical precision, thus several entries of As might not be exactly zero and require thresholding. To do so, we mask the inputs of each NN j using binary matrix Mg) E {0. 1] dxd initialized to have (Mo))a = Vi + and zeros everywhere else. Having (M = means the input of NN
5

---

Published as a conference paper at ICLR 2020
j has been thresholded. This mask is integrated in the product of Equation 6 by doing Cu) - ( W
( MM) without changing the interpretation of C (MG) can be seen simply as an extra layer in the NN). During optimization, if the entry (Ag) is smaller than the threshold c= =004, the corresponding mask entry (Mg)) is set to zero, permanently. The masks My) never updated via gradient descent. We also add an iterative thresholding step at the end to ensure the estimated graph Go is acyelic (described in Appendix A.2).
(5)
3.5 OVERFITTING
In practice, we maximize an empirical estimate of the objective of problem (10). It is well known that this maximum likelihood score is prone to overfitting in the sense that adding edges can never reduce the maximal likelihood (Koller & Friedman, 2009). GraN-DAG gets around this issue in four ways. First, as we optimize a subproblem, we evaluate its objective on held-out data set and declare convergence once it has stopped improving. This approach is known early stop- ping (Prechelt, 1997). Second, to optimize (11), we use a stochastic gradient algorithm variant which is now known to have an implicit regularizing effect (Poggio al., 2018). Third, once we have thresholded our graph estimate to be a DAG, we apply a final pruning step identical to what is done in CAM (Bühlmann et al., 2014) remove spurious edges. This step performs a regression of each node against its parents and uses a significance test to decide which parents should be kept or not. Fourth, for graphs of 50 nodes or more, we apply a preliminary neighbors selection (PNS) be- fore running the optimization procedure as was also recommended in Bühlmann et al. (2014). This procedure selects a set of potential parents for each variables. See Appendix A.3 for details on PNS and pruning. Many score-based approaches control overfitting by penalizing the number of edges in their score. For example, NOTEARS includes the LI norm of its weighted adjacency matrix U in its objective. GraN-DAG regularizes using PNS and pruning for ease of comparision to CAM, the most competitive approach in our experiments. The importance of PNS and pruning and their ability to reduce overfitting is illustrated in an ablation study presented in Appendix A.3. The study shows that PNS and pruning are both very important for the performance of GraN-DAG in terms of SHD. but do not have a significant effect in terms of SID. In these experiments, we also present NOTEARS and DAG-GNY with PNS and pruning, without noting a significant improvement.
3.6 COMPUTATIONAL COMPLEXITY
To learn a graph, GraN-DAG relies on the proper training of neural networks on which it is built. We thus propose using a stochastic gradient method which is standard choice when it comes to NN training because it scales well with both the sample size and the number of parameters and it implic- itly regularizes learning. Similarly to NOTEARS, GraN-DAG requires the evaluation of the matrix exponential of As at each iteration costing Odd)). NOTEARS justifies the use of a batch proximal quasi-Newton algorithm by the low number of (dd) iterations required to converge. Since GraN- DAG uses a stochastic gradient method, one would expect it will require more iterations to converge. However, in practice we observe that GraN-DAG performs fewer iterations than NOTEARS before the augmented Lagrangian converges (see Table 4 of Appendix A. We hypothesize this is due to early stopping which avoids having to wait until the full convergence of the subproblems hence lim- iting the total number of iterations. Moreover, for the graph sizes considered in this paper (d < 100), the evaluation of (()) in GraN-DAG, which includes the matrix exponentiation, does not dominate the cost of each iteration (= 4% for 20 nodes and = 13% for 100 nodes graphs). Evaluating the ap- proximate gradient of the log- ligelikeld (costing (dd)) assuming a fixed minibatch size, NN depth
and width) appears to be of greater importance for d < 100.
4 EXPERIMENTS
In this section, we compare GraN-DAG to various basclines in the continuous paradigm, namely DAG-GNN (Yu et al.. 2019) and NOTEARS (Zheng et al., 2018), and also in the combinatorial paradigm, namely CAM (Bühlmann et al., 2014), GSF (Huang et al., 2018), GES (Chickering, 2003) and PC (Spirtes al., 2000). These methods are discussed in Section 5. In all experiments, each NN learned by GraN-DAG outputs the mean of a Gaussian distribution AG), i.c. Og) == (() and X2 1X,g N(Po), 2
Vj. The parameters ( are learned as well, but do not depend on

---

Published as a conference paper at ICLR 2020
the parent variables X (unless otherwise stated). Note that this modelling choice matches the nonlinear Gaussian ANM introduced in Section 2.2.
We report the performance of random graphs sampled using the Erlös-Rényí (ER) scheme de- scribed in Appendix A.5 (denoted by RANDOM). For each approach, we evaluate the estimated graph on two metries: the structural hamming distance (SHD) and the structural interventional dis- tance (SID) (Peters & Bühlmann, 2013). The former simply counts the number of missing, falsely detected reversed edges. The latter is especially well suited for causal inference since it counts the number of couples (i,j) such that the interventional distribution p(æyldo(X; = F)) would be miscalculated if we were use the estimated graph form the parent adjustement GSF, GES and PC output only a CPDAG, hence the need to report lower and an upper bound on the SID. See Appendix A.7 for more details on SHD and SID. All experiments were ran with publicly available code from the authors website. See Appendix A.8 for the details of their hyperparameters. In Appendix A.9, we explain how one could use a held-out data set to select the hyperparameters of score-based approaches and report the results of such procedure on almost every settings discussed the present section.
Note that
4.1 SYNTHETIC DATA
We have generated different data sel lypes which vary along four dimensions: data generating pro- cess, number of nodes, level of edge sparsity and graph type. We consider two graph sampling schemes: Erdös-Rényí (ER) and scale-free (SF) (see Appendix A.5 for details). For each data set type, we sampled 10 data sets of 1000 examples as follows: First, a ground truth DAG G is ran- domly sampled following either the ER or the SF scheme. Then, the data is generated according to specific sampling scheme.
The first data generating process we consider is the nonlinear Gaussian ANM (Gauss-ANM) in- troduced in Section 2.2 in which data is sampled following Xj == (x(x²g) + N₂ with mutually independent noises N₂ -~(00, of) Vj where the functions fy are independently sampled from Gaussian process with unit bandwidth RBF kernel and with ož sampled uniformly. As mentioned in Section 2.2, we know G to be identifiable from the distribution. Proposition indicates that the modelling choice of GraN-DAG together with this synthetic data ensure that solving (10) optimal- ity would recover the correct graph. Note that NOTEARS and CAM also make the correct Gaussian noise assumption, but do not have enough capacity to represent the Sj functions properly.
We considered graphs of 10, 20, 50 and 100 nodes. Tables 1 & 2 present results only for 10 and 50 nodes since the conclusions do not change with graphs of 20 or 100 nodes (see Appendix A.6 for these additional experiments). We consider graphs of d and 4d edges (respectively denoted by ERI and ER4 in the case of ER graphs). We report the performance of the popular GES and PC in Appendix A.6 since they are almost never on par with the best methods presented in this section.
Table I: Results for ER and SF graphs nodes with Gauss- data
ERI
ER4
SFI
SF4
SHD
SID
SHD
SID
SHD
SID
SHD
SID
GraN-DAG
1.7:-2,5
1.7+331
8.3
21,8
1.2
4.1:16.1
9,9
16,4
DAG-GNN
11.413.11
14.4
35.
81.914.7
9.9
29.7.115.8
20.8.1.9
48.4115.66
NOTEARS
12.2+2.9
36.6=13.1
32.6-3.2
79.0+4.1
10.7+2.2
32.0++1..3
20.8+2.7
49.8+15.6
CAM
1.1-1.11
1.1+2.44
12.2-2.7
33.9
1.4-1.6
5.4
9.8
19.3+7.55
GSF
6.5.12.6
[6.2.110.8
21.7.8.4
137.2.1192
1.8
[2.0.5.11
8.5.4..2
[13.2.16.8
17.7=12.3]
62.7+14.9] 76.6+77.0
6.9++6.2]
20.6+12.11
RANDOM
26.3++9.8
25.8- 10.4
31.8+5.0
25.1+10.2
10.5++00.5
28.5++..0
47.2+12.2
Table 2: Results for ER and SF graphs of 50 nodes with Gauss-ANM data
ERI
ER4
SFI
SF4
SHD
SHD
SID
SHD
SID
GraN-DAG
5.1+2.8
22.4
102.6+21.2
1060.1:109.4
25.5++6.2
90.0+18.9
111.3+12,3
271.2=6654
DAG-GNN
49.2:779
304.4:105.1
1919215.2
2146.2164
182.814229
144.9113.3
540.81151.11
NOTEARS
62.8:992
327.3=119.9
202.3 14.3
2149.1 1:76.3
57.74:3.5
195.74554.9
153.7+11..
558.4=153.5
CAM
4.3:1.9
22.01:1.9
98.8=20.7
1197.2:125.9
24.11:6.2
85.7131.9
111.2:13,3
320.7:152.6
25.611..1
121.1:123.1
$1.8
[906662214.7
3166:66.7
[85.8:1:29.9
[284.7880.2
79.2433.51
1030.2::122.61
147.3:1:49.9
379.9119831
RANDOM
535.7=401.2
272.3+
708.44-234.4
1921.344-03.5
514.0+:360.0
380.34-190.3
660.6+194.99
1198.943044.

---

Published as a conference paper at ICLR 2020
We now examine Tables 1 & 2 (the errors bars represent the standard deviation across datasets per task). We can see that, across all settings, GraN-DAG and CAM are the best performing methods, both in terms of SHD and SID, while GSF is not too far behind. The poor performance of NOTEARS can be explained by its inability to model nonlinear functions. In terms of SHD, DAG-GNN per- forms rarely better than NOTEARS while in terms of SID, it performs similarly to RANDOM in almost all cases except in scale-free networks of 50 nodes or more. Its poor performance might be due to its incorrect modelling assumptions and because its architecture uses strong form of parame- ter sharing between the f functions, which is not justified in a setup like ours. GSF performs always better than DAG-GNN and NOTEARS but performs as good as CAM and GraN-DAG only about half the time. Among the continuous approaches considered, GraN-DAG is the best performing on these synthetic tasks.
Even though CAM (wrongly) assumes that the functions fj are additive, i.e. f (1(029) = LiEng (1y(2)) Vj. it manages to compete with GraN-DAG which does not make this incorrect modelling assumption". This might partly be explained by a bias-variance trade-off. CAM is biased but has a lower variance than GraN-DAG due to its restricted capacity, resulting in both methods performing similarly. In Appendix A.4, we present an experiment showing that GraN-DAG can outperform CAM in higher sample size settings, suggesting this explanation is reasonable.
Having confirmed that GraN-DAG is competitive on the ideal Gauss-ANM data, we experimented with settings better adjusted to other models to see whether GraN-DAG remains competitive. We considered linear Gaussian data (better adjusted NOTEARS) and nonlinear Gaussian data with additive functions (better adjusted to CAM) named LIN and ADD-FUNC, respectively. See Ap- pendix A.5 for the details of their generation. We report the results of GraN-DAG and other base- lines in Table 12& & 13 of the appendix. On linear Gaussian data, most methods score poorly in terms of SID which is probably due to the unidentifiability of the linear Gaussian model (when the noise variances are unequal). GraN-DAG and CAM perform similarly to NOTEARS in terms of SHD. On ADD-FUNC, CAM dominates all methods on most graph types considered (GraN-DAG is on par only for the 10 nodes ERI graph). However, GraN-DAG outperforms all other methods which can be explained by the fact that the conditions of Proposition I are satisfied (supposing the functions Eieng fiz(X, (X;) can be represented by the NNs).
We also considered synthetic data sets which do not satisfy the additive Gaussian noise assumption present in GraN-DAG, NOTEARS and CAM. We considered two kinds of post nonlinear causal models (Zhang & Hyvärinen, 2009), PNL-GP and PNL-MULT (see Appendix A.5 for details about their generation). A post nonlinear model has the form X, == (X( x + N2) where N, is noise variable. Note that GraN-DAG (with the current modelling choice) and CAM do not have the representational power express these conditional distributions, hence violating an assumption of Proposition 1. However, these data sets differ from the previous additive noise setup only by the nonlinearity 9j. hence offering a case of mild model misspecification. The results are reported in Table 14 of the appendix. GraN-DAG and CAM are outperforming DAG-GNN and NOTEARS except in SID for certain data sets where all methods score similarly to RANDOM. GraN-DAG and CAM have similar performance on all data sets except one where CAM is better. GSF performs worst than GraN-DAG (in both SHD and SID) on PNI.-GP but not on PNI-MULT where it performs better in SID.
4.2 REAL AND PSEUDO-READ DATA
We have tested all methods considered so far on a well known data set which measures the expression level of different proteins and phospholipids in human cells (Sachs et al., 2005). We trained only on the n = 853 observational samples. This datasct and its ground truth graph proposed in Sachs et al. (2005) (11 nodes and 17 edges) are often used in the probabilistic graphical model literature (Koller & Friedman, 2009). We also consider pseudo-real data sets sampled from the SynTReN generator (Van den Buicke, 2006). This generator designed to create synthetic transcriptional regulatory networks and produces simulated gene expression data that approximates experimental data. See Appendix A.5 for details of the generation.
2 Althughh it is true that GraN-DAG does not wrongly assume that the functions f, are additive, it is not clear whether its neural networks can exactly represent functions sampled from the Gaussian process.
8

---

Published as a conference paper at ICLR 2020
In applications, it is not clear whether the conditions of Proposition 1 hold since we do not know whether (Px. E C. This departure from identifiable settings is an occasion to explore a different modelling choice for GraN-DAG, In addition to the model presented at the beginning of this section, we consider an alternative, denoted GraN-DAG++, which allows the variance parameters o% to depend the parent variables X through the NN, i.e. 0()) == (ي(r)) log (.))) Note that this is violating the additive noise assumption (in ANMs, the noise independent of the parent variables).
In addition to metrics used in Section 4. we also report SHD-C. To compute the SHD-C between two DAGs, we first map each of them to their corresponding CPDAG and measure the SHD between the two. This metric is useful to compare algorithms which only outputs a CPDAG like GSF, GES and PC to other methods which outputs a DAG. Results are reported Table
Table 3: Results real and pseudo-real data
Protein signaling data set
SynTReN (20 nodes) SHD
SHD
SHD-C
STD
SHD-C
SID
GraN-DAG
13
47
3..0
36.4
161.7:153.4
GraN-DAG++
13
10
48
33.7+3.77
39.4+:499
127.5+55.8
DAG-GNN
16
44
93.6+9.2
97.6++00.3
157.5+44.6
NOTEARS
21
44
151.8-28.2
156.14228.7
110.7::66.7
CAM
12
9
55
4..5
41.411.11
152.3:1:48
GSF
18
10
[44, 61]
61.8+9.6
63.3+1144
[76.7=5111, 109.9+39.9]]
GES
26
[34, 45]
82.6+99.3
85.6+10
[157.2+44.3, 168.8+47.8]
PC
17
[47. 62]
41225.11
42.414.66
[154.8=47.6, 179.3455.6]
RANDOM
21
20
60
84.71:53.8
86.7.1.55.8
175.8:164.7
First, all methods perform worse than what was reported for graphs of similar size in Section 4.1, both in terms of SID and SHD. This might be due to the lack of identifiability guarantees we face in applications. On the protein data set, GraN-DAG outperforms CAM in terms of SID (which differs from the general trend of Section 4.1) and arrive almost par in terms of SHD and SHD-C. On this data set, DAG-GNN has reasonable performance, beating GraN-DAG SID, but not in SHD. On SynTReN, GraN-DAr obtains the best SHD but not the best SID. Overall, GraN-DAG is always competitive with the best methods of each task.
5 RELATED WORK
Most methods for structure learning from observational data make use of some identifiability results similar to the ones raised in Section 2.2. Roughly speaking, there are two classes of methods: independence-based and score-based methods. GraN-DAG falls into the second class.
basre-based methods (Koller Friedman, 2009; Peters et al., 2017) cast the problem of structure learning as optimization problem over the space of structures (DAGs or CPDAGs). Many popular algorithms tackle the combinatorial nature of the problem by performing a form of greedy search. GES (Chickering, 2003) is a popular example. It usually assumes a linear parametric model with Gaussian noise and greedily search the space of CPDAGs in order to optimize the Bayesian infor- mation criterion. GSF (Huang et al., 2018), is based on the same search algorithm as GES, but uses a generalized score function which can model nonlinear relationships. Other greedy approaches rely on parametric assumptions which render & fully identifiable. For example, Peters & Bühlman (2014) relies on a linear Gaussian model with equal variances to render the DAG identifiable. RE- SIT (Peters et al., 2014), assumes nonlinear relationships with additive Gaussian noise and greedily maximizes an independence-basee based score. However, RESIT does not scale well to graph of more than 20 nodes. CAM (Bühlmann al., 2014) decouples the search for the optimal node ordering from the parents selection for each node and assumes an additive noise model (ANM) (Peters et al., 2017) in which the nonlinear functions are additive. As mentioned Section 2.3, NOTEARS, proposed in Zheng et al. (2018), tackles the problem of finding an optimal DAG as a continuous constrained optimization program. This a drastic departure from previous combinatorial approaches which en- ables the application of well studied numerical solvers for continuous optimizations. Recently, Yu et al. (2019) proposed DAG-GNN, graph neural network architecture (GNN) which can be used to learn DAGs via the maximization of evidence lower bound. By design, GNN makes use of parameter sharing which we hypothesize is not well suited for most DAG learning tasks. To the best of our knowledge, DAG-GNN is the first approach extending the NOTEARS algorithm for structure
9

---

Published as a conference paper at ICLR 2020
learning to support nonlinear relationships. Although Yu et al. (2019) provides empirical compar- isons to linear approaches, namely NOTEARS and FGS (a faster extension of GES) (Ramsey et al., 2017), comparisons to greedy approaches supporting nonlinear relationships such as CAM and GSF are missing. Moreover, GraN-DAG significantly outperforms DAG-GNN on our benchmarks. There exists certain score-basee based approaches which uses integer linear programming (IIP)) (Jaakkola al., 2010; Cussens, 2011) I) which internally solve continuous linear relaxations. Connections between such methods and the continuous constrained approaches are yet be explored.
When used with the additive Gaussian noise assumption, the theoretical guarantee of GraN-DAG rests on the identifiability of nonlinear Gaussian ANMs. Analogously to CAM and NOTEARS, this guarantee holds only if the correct identifiability assumptions hold in the data and if the score maximization problem is solved exactly (which is not the case in all three algorithms). DAG-GNN provides no theoretical justification for its approach. NOTEARS and CAM are designed to handle what is sometimes called the high--imensional setting in which the number of samples is signifi- cantly smaller than the number of nodes. Bühhmaan et al. (2014) provides consistency results for CAM in this setting. GraN-DAG and DAG-GNN were not designed with this setting in mind and would most likely fail if confronted to it. Solutions for fitting neural network on less data points than input dimensions are not common in the NN literature.
Methods for causal discovery using NNs have already been proposed. SAM (Kalainathan et al., 2018) learns conditional NN generators using adversarial losses but does not enforce acyclicity. CGNN (Goudet al., 2018), when used for multivariate data, requires initial skeleton to learn the different functional relationships.
GraN-DAG has strong connections with MADE (Germain al., 2015), method used learn distributions using a masked NN which enforces the so-called autoregressive property. The autore- gressive property and acyclicity are in fact equivalent. MADE does not learn the weight masking, it fixes it at the beginning of the procedure. GraN-DAG could be used with a unique NN taking as input all variables and outputting parameters for all conditional distributions. In this case, would be similar to MADE, except the variable ordering would be learned from data instead of fixed priori.
CONCLUSION
The continuous constrained approach to structure learning has the advantage of being more global than other approximate greedy methods (since updates all edges at each step based on the gradient of the score but also the acyclicity constraint) and allows to replace task-specific greedy algorithms by appropriate off-the-shelf numerical solvers. In this work, we have introduced GraN-DAG, a novel basee approach for structure learning supporting nonlinear relationships while leverag- ing a continuous optimization paradigm. The method rests on a novel characterization of acyclicity for NNs based on the work of Zheng et al. (2018). We showed GraN-DAG outperforms other gradient-based approaches, namely NOTEARS and its recent nonlinear extension DAG-GNN, on the synthetic data sets considered in Section 4.1 while being competitive on real and pseudo-real data sets of Section 4.2. Compared to greedy approaches, GraN-DAG is competitive across all datasets considered. To the best of our knowledge, GraN-DAG the first approach leveraging the continuous paradigm introduced in Zheng et al. (2018) which has been shown to be competitive with state of the art methods supporting nonlinear relationships.
ACKNOWLEDGMENTS
This research was partially supported by the Canada CIFAR AI Chair Program and by a Google Focused Research award. The authors would like to thank Rémi Le Priol, Tatjana Chavdarova, well as
Charles Guille-Escuref Nicolas Gagné and Yoshua Bengio for insightful discussions Alexandre Drouin and Florian Bordes for technical support. The experiments were in part enabled by computational resources provided by Calcul Québec, Compute Canada and Element AI.
REFERENCES
J. Alayrac, P. Bojanowski, N. Agrawal, J. Sivic, I. Laptev, and S. Lacoste- Julien. Learning from narrated instruction videos. TPAMI, 4009):214-2200 Sep. 2018.
10

---

Published as a conference paper at ICLR 2020
A.-L. Barabási. Scale-free networks: decade and beyond. Science, 2009.
A.-L. Barabási and R. Albert. Emergence of scaling random networks. Science, 1999.
James Bergstra and Yoshua Bengio. Random search for hyper-paramedte optimization. Journal of Machine Learning Research, 2012.
D.P. Bertsekas. Nonlinear Programming. Athena Scientific, 1999.
P. Bühlmann, J. Peters, and J. Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. Annals Statistics, 2014.
D.M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 2003.
A. Clauset, C.R. Shalizi, and M.E.J. Newman. Power-law distributions in empirical data. SIAM review, 2009.
J. Cussens. Bayesian network learning with cutting planes. In Proceedings of the 27th Conference on Uncertainty in Artificial Intelligence, 2011.
M. Germain, K. Gregor, I. Murray, and H. Larochelle. Made: Masked autoencoder for distribution estimation. In Proceedings of the 32nd International Conference on Machine Learning, 2015.
P. Geurts, D. Ernst, and L. Wehenkel. Extremely randomized trees. Machine learming, 2006.
X. Glorot and Y. Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, 2010.
X. Glorot, A. Bordes, and Y. Bengio. Deep sparse rectifier neural networks. In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, 2011.
I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016.
O. Goudet, D. Kalainathan, P. Caillou, D. Lopez-Paz, Guyon, and M. Sebag. Learning Functional Causal Models with Generative Neural Networks. Explainable and Interpretable Models in Computer Vision and Machine Learning. Springer International Publishing, 2018.
Biwei Huang, Kun Zhang, Yizhu Lin, Bernhard Schölkopf. and Clark Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Confer- ence on Knowledge Discovery &#38; Data Mining. 2018
Jaakkola, D. Sontag, A. Globerson, and M. Meila. Learning Bayesian Network Structure using LP Relaxations. In Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, 2010.
D. Kalainathan, O. Goudet, I. Guyon, D. Lopez-Paz, and M. Sebag. Sam: Structural agnostic model, causal discovery and penalized adversarial learning. arXiv preprint arXiv:1803 04999, 2018.
D. Koller and N. Friedman. Probabilistic Graphical Models: Principles and Techniques Adaptive Computation and Machine Learning. MIT Press, 2009.
S. Magliacane, van Ommen, T. Claassen, S. Bongers, Versteeg, and J.M. Mooij. Domain adaptation by using causal inference predict invariant conditional distributions. In Advances in Neural Information Processing Systems 31, 2018.
J. Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, 2nd edition, 2009.
J. Pearl. The seven tools of causal inference, with reflections on machine learning. Commun. ACM, 2019.
Peters and Bühlman. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 2014.

---

Published as a conference paper at ICLR 2020
J. Peters and P. Bühlmann. Structural intervention distance (sid) for evaluating causal graphs. Neural Computation, 2013.
J. Peters and P. Bühlmann. Structural intervention distance (SID) for evaluating causal graphs. Neural Computation, 2015.
Peters, D. Janzing, and B. Schölkopf. Causal inference on discrete data using additive noise models. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011.
J. Peters, J. M. Mooij, D. Janzing, and B. Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 2014.
J. Peters, D. Janzing, and B. Schölkopf. Elements of Causal Inference - Foundations and Learning Algorithms. MIT Press, 2017.
A. Poggio, K. Kawaguchi, Q. Liao, B. Miranda, Rosasco, X. Boix, J. Hidary, and H. Mhaskar. Theory of deep learning III: explaining the non-overfitting puzzle. arXiv preprint arXiv: 1801.00110 2018.
L. Prechelt. Early stopping but when? Neural Networks: Tricks of the Trade, volume 1524 of LNCS, chapter 2, 1997.
J. Ramsey, M. Glymour, R. Sanchez-Romero, and C. Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimennional graphical causal models, with an application to functional magnetic resonance images. J. Data Science and Analytics, 2017.
K. Sachs, O. Perez, D. Pe'er, D.A. Lauffenburger, and G.P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 2005.
Shimizu, P.O. Hoyer, A. Hyvärimen, and A. Kerminen. linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 2006.
Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT press, 2000.
Tieleman and Hinton. Lecture 35-Msmrrppp Divide the gradient by a running average of its recent magnitude. COURSERA: Neural Networks for Machine Learning, 2012.
Tim Van den Bulcke, et al. SynTReN: generator of synthetic gene expression data for design and analysis of structure learning algorithms. 2006.
Y. Yu, J. Chen, T. Gao, and M. Yu. DAG-GNN: DAG structure learning with graph neural networks. In Proceedings of the 36th International Conference Machine Learning, 2019.
K. Zhang and Hyvärinen. On the identifiability of the post-nonlinear causal model. Proceed- ings of the 25th Conference on Uncertainty in Artificial Intelligence, 2009.
Kun Zhang, Zhikun Wang, Jiji Zhang, and Bernhard Schölkopf. On estimation of functional causal models: General results and application to the post-nonlinear causal model. ACM Trans. Intell. Syst. Technol., 2015.
X. Zheng, B. Aragam, P.K. Ravikumar, and E.P. Xing. Dags with tears: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems 31, 2018.
12

---

Published as a conference paper at ICLR 2020
A APPENDIX
A.1 OPTIMIZATION
Let us recall the augmented Lagrangian:
max L(d, Aos (b)) - ExxPx [xix (3 - Ap((6) - Ht h()) 2
(12)
where D and j4 are the Lagrangian and penalty coefficients of the ith subproblem, respectively. In all our experiments, mately solve each non-convex subproblem using RMSprop (Tieleman & Hinton, 2012), a stochastic gradient descent variant popular for NNs. We use the following gradient estimate:
initialize those coefficients using λo and μl = 10 We approxi-
V Vol(o, da, (4) VoLB(o As,
d
(13)
with A,
log logni(mll Ach(d) - 14 h(0)² 2
B
ZEB
where B is minibatch sampled from the data set and B is the minibatch size. The gradient estimate V VoLa de: be computed using standard deep learning libraries. We consider subproblem has converged when Ln(0), At, (4. evaluated on a held-out data set H stops increasing. Let of be the approximate solution to subproblem I. Then, d and μe are updated according to the following rule:
Aet1 < de Meh (ot)
Ma41 + {
{ MH if h 7h (ii-1)
(14)
He, otherwise
with 77 = 10 and y = 0.9. Each subproblem is initialized using the previous subproblem solution 1- The augmented Lagrangian method stops when h(6) < 10-8
Total number of iterations before augmented Lagrangian converges: In GraN-DAG and NOTEARS, every subproblem is approximately solved using iterative algorithm. Let T be the number of subproblems solved before the convergence of the augmented Lagrangian. For given subproblem f, let K, be the number of iterations executed to approximately solve Let I --= K, be the total number of iterations before the augmented I Lagrangian converges. Ta- ble 4 reports the total number of iterations for GraN-DAG and NOTEARS, averaged over ten data sets. Note that the matrix exponential is evaluated once per itcration. Even though GraN-DAG uses stochastic gradient algorithm. it requires less iterations than NOTEARS which uses batch proximal quasi-Newtor method. We hypothesize early stopping avoids having to wait until full convergence before moving the next subproblem, hence reducing the total number of iterations. Note that GraN-DAG total run time is still larger than NOTEARS due to its gradient requiring more computation evaluate (total runtime 10 minutes against R 1 minute for 20 nodes graphs and s 4 hours against z hour for 100 nodes graphs). GraN-DAG runtime on 100 nodes graphs can be roughly halved when executed on GPU.
Table 4: Total number of iterations (x (x103) before augmented Lagrangian converges Gauss-ANM data.
20 nodes ERI
20 nodes ER4
100 nodes ER1
100 nodes ER4
GraN-DAG
27.3 + 3.6
30.4 + 4.2
23.1 + 0.7
23.1 + 0.8
NOTEARS
67.1 = 35.3
72.3 + 24.3
243.6 12.3
232.4 + 12.9
A.2 THRESHOLDING TO ENSURE ACYCLICITY
The augmented Lagrangian outputs of where T is the number of subproblems solved before declar- ing convergence. Note that the weighted adjacency matrix Ag; will most likely not represent an acyclic graph, even if threshold as learn, as explained in Section 3.4. We need remove
13

---

Published as a conference paper at ICLR 2020
additional edges to obtain DAG (edges are removed using the mask presented in Section 3.4). One option would be to remove edges one by one until a DAG is obtained, starting from the edge (i.j) with the lowest (Asz. )e up to the edge with the highest (Agz )ay. This amounts to gradually increas- ing the threshold until Aaz is acyelic. However, this approach has the following flaw: It is possible to have (Aa;)s Jay significantly higher than zero while having (()) almost completely independent of variable X,. This can happen for least two reasons. First, the NN paths from input i to output k might end up cancelling each others, rendering the input inactive. Second, some neurons of the NNs might always be saturated for the observed range of inputs, rendering some NN paths effec- rively inactive without being inactive in the sense described in Section 3.1. Those two observations illustrate the fact that having (Asr, )az 0 is only a sufficient condition have ( independent of variable X; and not necessary one.
To avoid this issue, we consider the following alternative. Consider the function L : Rdd I2d which maps all d variables to their respective conditional likelihoods, i.e. (i(X) - px(:; X 8,
o5) Vi.
We consider the following expected Jacobian matrix
0L
T 2
Ex~Px
(15)
0x
where 1x is the Jacobian matrix of L evaluated X, in absolute value (element-wise). Similarly to (Ao; lij. the entry Ti can be loosely interpreted as the strength of edge (i,j). We propose removing edges starting from the lowest J.y to the highest, stopping as soon as acyclicity is achieved. We believe I is better than Age, at capturing which NN inputs are effectively inactive since it takes into account NN paths cancelling each others and saturated neurons. Empirically, we found that using J instead of As; yields better results, and thus we report the results with J in this paper.
A.3 PRELIMINARY NEIGHBORHOOD SELECTION AND DAG PRUNING
For graphs of 50 nodes more, GraN-DAG performs preliminary neighborhood selec- tion (PNS) similar to what has been proposed in Bühlmann et al. (2014). This procedure applies a variable selection method to get a set of possible parents for each node. This is done by fit- ting an extremely randomized trees (Geurts et al., 2006) (using Ext rerrrrreeeeeeeeeooo from for each variable against all the other variables. For each node feature impor- tance score based on the gain of purity is calculated. Only nodes that have a feature importance score higher than 0.75 mean are kept potential parent, where mean is the mean of the feature impor- tance scores of all nodes. Although the use of PNS in CAM was motivated by gains in computation time, GraN-DAG uses it to avoid overfitting, without reducing the computation time.
PNS:
Once the thresholding is performed and a DAG is obtained as described in A.2, GraN- DAG performs a pruning step identical to CAM (Bühlmann al., 2014) in order to remove spurious edges. We use the implementation of Bühlmann et al. (2014) based on the R function gamboost from the mboost package. For each variable X,, generalized additive model is fitted against the current parents of X, and a significance test of covariates is applied. Parents with a p-value higher than 0.001 are removed from the parent set. Similarly to what Bühlmann al. (2014) observed, this pruning phase generally has the effect of greatly reducing the SHD without considerably changing the SID.
Pruning:
Ablation study: In Table 5, we present an ablation study which shows the effect of adding PNS and pruning to GraN-DAG on different performance metrics and on the negative log-likelihood (NLL) of the training and validation set. Note that, before computing both NLL, we reset all parameters of GraN-DAG except the mask and retrained the model on the training set without any acyclicity constraint (acyclicity is already ensure by the masks at this point). This retraining procedure is important since the pruning removes edges (i.e. some additional NN inputs are masked) and affects the likelihood of the model (hence the need retrain).
14

---

Published as a conference paper at ICLR 2020
Table 5: PNS and pruning ablation study for GraN-DAG (averaged 10 datasets from ER1 with 50 nodes)
PNS
Pruning
SHD
SID
NLL (train)
NLL (validation)
False
False
1086.84-48.88
31.6++23.6
0.36
1.44
True
False
540.4=70.3
17.4=16.7
0.52
1.16:0.17
False
True
11.815.0
39.7=25.5
0.78+0.112
0.84+0.12
True
True
6.1
29.3
0.78+00.1
0.83+0.12
first observation is that adding PNS and pruning improve the NLL on the validation set while de- teriorating the NLL on the training set, showing that those two steps indeed reducing overfitting. Secondly, the effect on SHD is really important while the effect on SID is almost nonexistent. This be explained by the fact that SID has more do with the ordering of the nodes than with false positive edges. For instance, if we have a complete DAG with node ordering coherent with the ground truth graph, the SID zero, but the SHD is not due all the false positive edges. Without the regularizing effect of PNS and pruning, GraN-DAG manages to find a DAG with a good ordering but with many spurious edges (explaining the poor SHD, the good SID and the big gap between the NLL of the training set and validation set). PNS and pruning helps reducing the number of spurious edges, hence improving SHD.
We also implemented PNS and pruning for NOTEARS and DAG-GNN to see whether their perfor- mance could also be improved. Table reports an ablation study for DAG-GNN and NOTEARS. First, the SHD improvement is not as important significant. The improved SHD does not come close to performance of GraN-DAG, Second, PNS and pruning do not have significant effect of SID, improvement for those methods is probably due to the fact that they are not overfitting like GraN- DAG, as the training and validation (unregularized) scores shows. NOTEARS captures only linear relationships, thus it will have hard time overfitting nonlinear data and DAG-GNN uses a strong form of parameter sharing between its conditional densities which possibly cause underfitting in setup where all the parameters of the conditionals are sampled independently.
for GraN-DAG and is almost not statistically
the for GraN-DAG. The lack of
Moreover, DAG-GNN and NOTEARS threshold aggressively their respective weighted adjacency matrix at the end of training (with the default parameters used in the code), which also acts as form of heavy regularization, and allow them to remove many spurious edges. GraN-DAG without PNS and pruning does not threshold as strongly by default which explains the high SHD of Table 5. To test this explanation, we removed all edges (i, ji for which (Aa)iz < 0.3* for GraN-DAG and obtained an SHD of 29, 15.99 and an SID of 85.6.145.7, showing a significant improvement over NOTEARS and
even without PNS and pruning.
Table 6: PNS and pruning ablation study for DAG-GNN and NOTEARS (averaged over 10 datasets from ERI with 50 nodes)
Algorithm
PNS
Pruning
SHD
SID
Score (train)
Score (validation)
DAG-GNN
False
False
56.8=11.1
322.9+103.8
-2.8
-2.2+1.6
True
False
55.5+11.2
314.5+107.6
-2.111.6
-2.11177
False
True
49.4-17.8
325.1+103.7
-1.8+1.1
-1.8+1.2
True
True
47.7+7.3
316.5+105.6
-1.9+1.6
-1.9+1.6
NOTEARS
False
False
64.2499.5
327.1+110.9
-3.1111.8
-3.22-2.11
True
False
54.1+10.9
321.5+
.25.2+2.7
-25.4+2.8
False
True
49.5+1.88
327.7+111.3
-26.742.0
-26.822.1
True
True
49.017.6
326.41106.9
-26.2312.2
-26.412.4
Thii the default value of thresholding used in NOTEARS and DAG-GNN.
15

---

Published as a conference paper at ICLR 2020
A.4 LARGE SAMPLE SIZE EXPERIMENT
In this section, we test the bias-variance hypothesis which attempts to explain why CAM is on par with GraN-DAG on Gauss-ANM data even if its model wrongly assumes that the ʃj functions are additive. Table 7 reports the performance of GraN-DAG and CAM for different sample sizes. We can see that, as the sample size grows, GraN-DAG ends up outperforming CAM in terms of SID while staying on par in terms of SHD. We explain this observation by the fact that a larger sample size reduces variance for GraN-DAG thus allowing to leverage its greater capacity against CAM which is stuck with its modelling bias. Both algorithms were run with their respective default hyperparameter combination.
This experiment suggests GraN-DAG could be an appealing option in settings where the sample size substantial. The present paper focuses on sample sizes typically encountered in the struc- ture/causal learning litterature and leave this question for future work.
Table Effect of sample size
Gauss-ANM 50 nodes ER4 (averaged over 10 datasets)
Sample size
Method
SHD
SID
500
CAM
123.5 + 13.9
1181.2 + 160.8
GraN-DAG
130.2 = 14.4
1246,4 + 126.1
1000
CAM
103.7 = 15.2
1074.7 125,8
GraN-DAG
104.4 15.3
942.1 + 69.8
5000
CAM
74.1
845.0 159.8
GraN-DAG
7
554.1 117.9
10000
CAM
66.3 16.0
808. + 142.9
GraN-DAG
65.9 + 19.8
453.4 + 171.7
A.5 DETAILS ON DATA SETS GENERATION
Synthetic data sets: For each data set type, 10 data sets are sampled with 1000 examples each. As the synthetic data introduced Section 4.1, for each data set, a ground truth DAG g is randomly sampled following the ER scheme and then the data is generated. Unless otherwise stated, all root variables are sampled from 11.
Gauss-ANM is generated following X, == [x(x,)) + N2 Vj with mutually independent noises Ny ~ (N(, o3) Vj where the functions fJ are independently sampled from a Gaus- sian process with unit bandwidth RBF kernel and of ~ 440044 0.8). Source nodes are Gaussian with zero mean and variance sampled from [1,1,
LIN is generated following Xy XylX,3 w/ X, + 0.2 (O o3) Vj where of ~ u[1,
and wy is vector of |T/] coefficients each sampled from 4(10, 1].
ADD-FUNC is generated following X,lX,ş ~ Exery f5.((,,) + 0.2 -N of) Vj
where of ~ U[1, 2] and the functions f3, are independently sampled from a Gaussian process with bandwidth one. This model is adapted from Bühlmann al. (2014).
PNL-GP is generated following X,IX,9 (X5g(),)) Laplace(0, 4)) Vj with the functions f, independently sampled from a Gaussian process with bandwidth one and 1j 00 1]. In the two-variable case, this model is identifiable following the Corollary 9 from Zhang & Hyvärinen (2009). To get identifiability according to this corollary, is important
non non-Gausiian noise, explaining our design choices.
PNL-MULT is generated following X,|X.g exp(log(Eieeg) X,) + LV((,o3))) Vj where of u10 1]. Root variables are sampled from u(0, 2). This model is adapted from Zhang et al. (2015).
SynTReN:
Ten datasets have been generated using the SynTReN generator (http: //bioiifooratiicomm intec. ugent. be/ /ma rcha htttk//aaaaaaaaneee//ea///aaa/// innerReNN/indee html) using the
16

---

Published as a conference paper at ICLR 2020
software default parameters except for the probability for complex 2-regulator interactions that was set to 1 and the random seeds used were 0 to 9. Each dataset contains 500 samples and comes from 20 nodes graph.
Graph types: Erdös-Rényi (ER) graphs are generated by randomly sampling a topological order and by adding directed edges were it is allowed independently with probability p = 2²-d were is the expected number of edges in the resulting DAG. Scale-free (SF) graphs were generated using the Barabási-A Albert model (Barabási & Albert, 1999) which is based on preferential attachment. Nodes Between the
added one by to d or 4d) will be added. An existing node i have the probability p(k,) Σ, to be chosen, where k; represents the degree of the node i. While ER graphs have degree distribution following Poisson distribution, SF graphs have degree distribution following power law: few nodes, often called hubs, have a high degree, Barabási (2009) have stated that these types of graphs have similar properties to ceal-world networks which can be found in many different fields, although these claims remain controversial (Clauset et al., 2009).
node and the existing nodes, edges (where is equal
ks
A.6 SUPPLEMENTARY EXPERIMENTS
Gauss-ANM: The results for 20 and 100 nodes are presented in Table 8 and using the same Gauss- ANM data set types introduced in Section 4.1. The conclusions drawn remains similar to the 10 and 50 nodes experiments. For GES and PC, the SHD and SID are respectively presented in Table 10 and 11. Their performances do not compare favorably to the GraN-DAG nor CAM. Figure 1 shows the entries of the weighted adjacency matrix Ag as training proceeds in a typical run for 10 nodes.
Experiments with LIN and ADD-FUNC data is reported in Table 12 & 13. The details of their generation are given in Appendix A.5.
LIN & ADD-FUNC:
PNL-GP & PNL-MULT: Table 14 contains the performance of GraN-DAG and other baselines on post nonlinear data discussed in Section 4.1.
Table 8: Results for ER and SF graphs of 20 nodes with Gauss-ANM data
ERI
ER4
SFI
SF4
SHD
SID
SHD
SID
SHD
SID
SHD
SID
GmN-DAG
4.0 +34
19.5
45.24-10.7
165.1+21.00
7.6+
28.84-10.4
36.8-5.11
62.5 18.88
DAG-GNN
25.6=7.5
109.1=53.1
75.0477.7
14.
19.5+11.8
60.1+:1.88
49.5
115.24433.3
NOTEARS
30,3=7.8
107.3247.6
79.018.0
346.6213.2
23.9::3.5
69.4119.7
$2.0:14.5
120.5::32.5
CAM
2.7+11.8
10.6=8.66
41.0+1119
157.9::41.2
5.7+-2.6
23.3+-18.0
59.1+1888
GSF
12.3=4.6
[15:0::19.9
4188+3.8
[153.1149.4
7,4:3
[5.7=7.1
38.6::3.6
15,,9=14,4
45.6::22.9
201.6:37.91
27.3:113.2]
86.7-:24.1
RANDOM
103.04399.6
94.3+53.0
117.5:25.9
298.5428.7
105.2=48.8
81.1+:5.44
121.5=28.5
204.8:38.5
Table 9: Results for ER and SF graphs of 100 nodes with Gauss-ANM data
ERI
ERI
SEI
SF4
SID
4207.3::49.7
59.2-7.77
DAG-ENN
110.--10.5
883.0+3320.9
8096.1+656
97.6-1.5
438.6-112.7
316.0+14.3
1394.6-166
NOTEARS
125
913.1+334..8
387 84:25
81247457714
111.7=54
484 334.3-138.44
527. 241..
1=42.:=210.1
CAM
17.3:45
124.9::65.0
186.44:28.8
46011.92482.7
$2.7=:9.3
230.3::36.
288.6::21.7
845.8=161.35
GSF
66.8273
[104,7159.5
12 hours"
714:1112
[212.7171.1
275.9121.0
17919-15225
238.6.159.31
193.4..149.1
RANDOM
1561.6511334
1175.35557.9
2380.9:1458.0
10260
1164.24:93.3
2485.0+1403.9
42066.4s
'Note that GSF results are missing for two data set types in Tables and 14. This is because the search algorithm could not finish within 12 hours, even when the maximal in-degree was limited to 5. All other methods could run in less than 6 hours.
17

---

Published as a conference paper at ICLR 2020
Table 10: SHD for GES and PC (against GraN-DAG for reference) with Gauss-ANM data
10 nodes
20 nodes
50 nodes
100 nodes
ERI
ER4
ER1
BR4
ERI
ER4
ER1
ER4
BraG-DAG
1.7=2.5 13.824.8
8,34-288
4.0 +33.4
45.2= 10.7
5.1+2.88
102.6+21.2
15.1+6.00
206.6++31.5
GES
32.34-4.3
43.3=124
94.6=9.8
106.6=24.7
254.44:39.3
292.9:33.6
542.6:55.22
PC
8.423
34:2.6
20.136/41:.5
73.1=5.8
44.0::16
183.6::20
95.249.1
35..8:220.6
SF1
SF4
SFI
SF4
SFI
SF4
SFI
SF4
GraN-DAG
1.2-11
9,944.00
7.642.5
36.8=5.1
25.5:6.2
111.3412.3
59.247.7
192.7:19.66
GES
8.1=2.4
17,444.55
26.2:7.55
50.7=6.2
73.9:7.44
178.8+16..5
190.3:22
408.7::2499
PC
4.8=2.4
16.4422.
13.6+2.1
44.4=4.6
43.1+5.7
13.
97.6++6.6
314.2++1755
Table 11: Lower and upper bound on the SID for GES and PC (against GraN-DAG for reference) with Gauss-ANe ANM data. See Appendix A.7 for details on how to compute SID for CPDAGs.
10 modes
20 nodes
50 oodes
100 nodes
ERI
ER4
ERI
ER4
ERI
ERI
ER4
GraN-DAG
1.743.1
21.8-89
17.9-195
165.1421.0
22.411788
1060. 109.4
83.9-46.0
4207.3+-419.7
GES
124.1117.3
68.55100.
162.1144
9301
115091727
6996.6113.1
582.519911
5054
27.2417.51
75471
65.724455)
3193313333
158.1
2032.9488.7
598.94408.61 6]
8124.24440..2
122.61155
[28.1=74
(809=511
1316,72257
1222
[21679488
[630/742409
41786994488.5
PC
23.3++1.11
79.2-5.7)
94.9-46.11
328.7+-5.66
256.7 127.31
80.8)
702.5+
8265.4-1470.2)
GraN-DAG
4.1266.1
1644660
28.8210.4
62.541888
90.02:1899
271.226544
265.4261122
872.0=13044
[171617761
[511.51257.6
1142).7=247.4
GES
16.4+11.7)
43.9-144
57.9+2446]
105. 14-414)
182.7477]
$29.7- 184.5)
5244-552.2
1485,44123.66
18.3240
[363=62
142:2114
(95.6:37
(124.24333
[4832=1133 115.
[414,544244
[1369 2:9002229900
16.8:12.3]
41.7=6.91
59.7=1499
118.5230]
167.124114]
538:143.71
186.5:120.91
1513.7:296.21
10-1
10-2
Entries of Ae
10-3
10
10-5
0
10000
20000
30000
40000
Number gradient steps
Figure Entries of the weighted adjacency matrix Ag as training proceeds in GraN-DAG for synthetic data set ER4 with 10 nodes. Green curves represent edges which appear in the ground truth graph while red ones represent edges which do not. The horizontal dashed line at 10- the threshold introduced in Section 3.4. We can see that GraN-DAG successfully recovers most edges correctly while keeping few spurious edges.
18

---

Published as a conference paper at ICLR 2020
Table 12: LIN
Graph Type
FRI
ER4
ERI
ER4
Metrics
SHID
SID
SHDD
SID
SHID
SID
SHID
SID
GraN-DAG
7.2+-2.0
27.3 8.1
30.7 + 3.3
70.8 6.9
33.9 + 8.6
265.8 + 158.4
181.9 24.0
21135.8 + 137.2
DAG-GNN
39.6
NOTEARS
9.0+33.
35.3 13.4
27.9
4.3
72.1 7.9
45.5++72 7.8
310.7
125.9
126 + 13.0
1971.1
134.3
CAM
1..2.1.6.3
31.2 10.9
33.6
1:
77.5 2.3
36.2
234.8
105.1
182.5
17.8
1948.7
113.5
GSF
9.2
[19.3 1 14.6
38.6
3.7
[738 1: 7.6
46.7 4.1
[176.4
98.8
> hours
31.6 = 17.3
85.2+-8.3 8.3
215.0 + 108.9)
RANDOM
34.4
Table 13: ADD-FUNC
Nodes
10
Graph Type
ERI
ER4
ERI
ER4
Metrics
SHD
SID
SHD
SID
SHD
SID
SHD
SID
Method
GraN-DAG
2.8
32.6
16.6
103.6
DAG-GNN
10.1=3.
23.3 11.5
18.3
661
45.5 7.9
88.8
221.3 31.6
1741.0
136.3
NOTEARS
11.1-5.0
16.9 11.3
20.3
53.5 10.5
53.7 + 9.5
276.1 96.8
201.8 22.1
1813.6
148.4
CAM
2.5
6.4
6.6
29.3
19.3
9.6 1.1
39.0 34.1
42.9 6.6
857.0
184.5
GSF
9.3++3.
113.9
2..5
60.3
11.6
49.5
[151.5 73.6
12 hours
24.1
12.55
75.0 7..0
213.9
62.5
RANDOM
2..0=22
26.9 18.1
76.0 6.2
689.7 1 6.1
3100
2: 113.6
711.5 1:
1916.2 63.8
Table 14: Synthetic post nonlinear data sets
PNI--GP
PNI.-MULT
SHD
STD
SHD
SID
10 nodes ERI
GraN-DAG
1.6+3.0
3,9+8.0
13.1=3.8
35.7+12.3
DAG-GNN
11 11.5+668
32.4+19.3
17.900+6.2
40.700+
NOTEARS
10.7:15.5
34.4:119.1
1..0
38.61:11.9
CAM
1.5+2.6
6.8+12.1
12.0=6.4
36.3+17.7
GSF
6.2+-3.3
[7.7+8.7, 18.9+12.4]
10.7+3.0
[9.8+119, 25.3+11.55
RANDOM
23.8+-2.9
36.8-t19.1
23.7++2.99
37.7++20.7
TO nodes ER4
GraN-DAG
10.916.8
39,8111111
12.1
71.115.9
DAG-GNN
32.31:4.3
8.1..3
37.013.1
82.716.4
NOTEARS
34.11:3.2
8U......55
37.723.0
81.700
CAM
8.4144.8
1
34.41:399
79.6
25.0.16.0
[44.3 66.1
31.315.4
[58.6 8.1, 76.4.9.9]
RANDOM
35.01:3.3
80.015.11
33.62:3.5
76.247.3
50 nodes ERI
GraN-DAG
16.5=7.0
64.11:35.4
38.2=11.4
213.8=114.4
DAG-GNN
56.5=11.1
334.3++-0.3
83,9=23.8
205.7
NOTEARS
50.1+-9.9
319.1+
78.5+:21.5
425.7+197.0
CAM
5,11:2.6
10.7+1224
44.9-9.9
284.3
31.2:+6.0
[59.5+34.1,
122.4=32.0]
46.3-11.11
165.8-62.2,
141.6+72.6]
RANDOM
688.4+44.9
691.3+77.3
488.0+247.8
50 nodes FR4
GraN-DAG
68.7+17.0
1127.0+188% 188.5
211.7+11.66
2047.7+77.7
DAG-GNN
203.8:188.9
2173.1187.77
166.
42.33
NOTEARS
189,5116.0
2134,2.112.66
220.0
58.3
CAM
48.2+-10.3
899.5+195.6
208. 108.1+114.8
2029.7 55.4
GSF
105.2+15..5
[1573.7:121.2
16208
hours
RANDOM
722.319.0
4183.7
710.2.19.5
1935.81.569
A.7 METRICS
SHD takes two partially directed acyclic graphs (PDAG) and counts the number of edge for which the edge type differs both PDAGs. There are four edge types: i f j. i + j. - j and j. Since this distance is defined over the space of PDAGs, DAGs with CPDAGs and CPDAGs with CPDAGs. When comparing a DAG with a CPDAG, having i +
use it compare DAGs with DAGs,
instead of i j counts mistake.
SHD-C is very similar to SHD. The only difference is that both DAGs are first mapped to their respective CPDAGs before measuring the SHD.
19

---

Published as a conference paper at ICLR 2020
Introduced by Peters & Bühlmann (2015), SID counts the number of interventional distribution of the form p(z;| do(z, = 8z)) that would be miscalculated using the parent adjustment for- mada (Pearl, 2009) if we were to the predicted DAG instead of the ground truth DAG to form the parent adjustment set. Some care needs to be taken to evaluate the SID for methods outputring a CPDAG such as GES and PC. Peters & Bühlmann (2015) proposes to report the SID of the DAGs which have approximately the minimal and the maximal SID in the Markov equivalence class given by the CPDAG. See Peters & Bühlmann (2015) for more details.
A.8 HYPERPARAMETERS
All GraN-DAG runs up to this point were performed using the following set of hyperparameters. We used RMSprop optimizer with learning rate of 10 for the first subproblem and 10-4 for all subsequent suproblems, Each NN has two hidden layers with 10 units (except for the real and pseudo-real data experiments of Table 3 which uses only hidden layer). Leaky-Rell is used as activation functions. The NN are initialized using the initialization scheme proposed in Glorot & Bengio (2010) also known Xavier initialization. We used minibatches of 64 samples. This hyper- parameter combination have been selected via a small scale experiment in which many hyperparam- eter combinations have been tried manually on a single data set of type ERI with 10 nodes until one yielding satisfactory SHD was obtained. Of course in practice one cannot select hyperparameters in this way since we do not have access to the ground truth DAG. In Appendix A.9, we explain how one could use held-out data set to select the hyperparameters of score-based approaches and report the results of such a procedure on almost settings presented in this paper.
For NOTEARS, DAG-GNN, and GSF, we used the default hyperparameters found in the authors code. (rarely) happens that NOTEARS and DAG-GNT returns a cyclic graph. those cases, we removed edges starting from the weaker ones to the strongest (according to their respective weighted adjacency matrices), stopping as soon as acyclicity is achieved (similarly to what was explained in Appendix A.2 for GraN-DAG). For GES and PC, we used default hyperparameters of the pealg R package. For CAM, we used the the default hyperparameters found in the CAM R package, with default PNS and DAG pruning.
A.9 HYPERPARAMETER SELECTION VIA HELD-OUT SCORE
Most structure/causal learning algorithms have hyperparameters which must be selected prior to learning. For instance, NOTEARS and GES have regularizing term in their score controlling the sparsity level of the resulting graph while CAM has a thresholding level for its pruning phase (also controlling the sparsity of the DAG). GraN-DAG and DAG-GNN have many hyperparameters such as the learning rate and the architecture choice for the neural networks (i.c. number of hidden layers and hidden units per layer). One approach to selecting hyperparameters in practice consists in trying multiple hyperparameter combinations and keeping the one yielding the best score evaluated on a held-out set (Koller& & Friedman, 2009, p. 960). By doing so, one can hopefully avoid finding DAG which is too dense or too sparse since if the estimated graph contains many spurious edges, the score on the held-out data should be penalized. In the section, we experiment with this approach on almost all settings and all methods covered in the present paper.
Experiments: We explored multiple hyperparameter combinations using random search (Bergstra & Bengio, 2012). Table 15 to Table 23 report results for each dataset types. Each table reports the SHD and SID averaged over 10 data sets and for each data set, we tried 50 hyperparameter combi- nations sampled randomly (see Table 24 for sampling schemes). The hyperparameter combination yielding the best held-out score among all 50 runs is selected per data set (i.e. the average of SHD and SID scores correspond to potentially different hyperparameter combinations on different data sets). 80% of the data was used for training and 20% was held out (GraN-DAG uses the same data for carly stopping and hyperparameter selection). Note that the held-out score is always evaluated without the regularizing term (c.g. the held-out score of NOTEARS was evaluated without its LI regularizer).
The symbols ++ and + indicate the hyperparameter search improved performance against default hyperparameter runs above one standard deviation and within standard deviation, respectively. Analogously for and which indicate a performance reduction. The flag ... indicate that, on average, less than 10 hyperparameter combinations among the 50 tried allowed the method to
20

---

Published as a conference paper at ICLR 2020
converge in less than 12 hours. Analogously, =* indicates between 10 and 25 runs converged and indicates between 25 and 45 runs converged.
Discussion: GraN-DAG and DAG-GNN are the methods benefiting the most from the hyperparam- eter selection procedure (although rarely significantly). This might be explained by the fact that neural networks in general very sensitive to the choice of hyperparameters. However, not all methods improved their performance and no method improves and GSF for instance, often have significantly worse results. This might be due some degree of model misspecification which renders the held-out score a poor proxy for graph quality. Moreover, for some methods the gain from the hyperparameter tuning might be outweighed by the loss due to the 20% reduction in training samples.
performance in all settings. GES
Additional implementation details for held-out score evaluation: GraN-DAG makes use of a fi- nal pruning step to remove spurious edges. One could simply mask the inputs of the NN correspond- ing removed edges and evaluate the held-out score. However, doing so yields an unrepresentative score since some masked inputs have important role in the learned function and once these inputs are masked, the quality of the might greatly suffer. To avoid this, we retrained the whole model from scratch on the training set with the masking fixed to the one recovered after pruning. Then, evaluate the held-out score with this retrained architecture. During this retraining phase. the esti- mated graph is fixed, only the conditional densities relearned. Since NOTEARS and DAG-GNN are not al always guaranteed to return a DAG (although they almost always do). some extra threshold- ing might be needed as mentioned in Appendix A.8. Similarly to GraN-DAG's pruning phase, this step can seriously reduce the quality of the fit. To avoid this, we also perform a retraining phase for NOTEARS and DAG-GNN. The model of CAM also retrained after its pruning phase prior to evaluating its held-out score.
Table 15: Gauss-ANM- 10 nodes with hyperparameter search
Graph Type
SF4
Method
GraN-DAG
1011
1.3
9.6 4.5*
DAG-GNN
2.6
95.5 1: 13.6
38.3 2.9
8.4
30.3
21 2.1
41.0 15.5
NOTEARS
6.9
35.2 10.6*
20.9 6.6*
62.0 62.0+6.77
38.8 + 16.7
26.9 7.4"
61.1 13.8
GSF
5.3+-33
[8.3 13.2+
23.1 7.9
[66.1 20.4
3.3
17.0
14.2 5.6
26..2 11.1
15.4 13.5
66.1 19.3
12.2
36.9 00:121..0
GES
15.4*
7.0
5.2
32.7
28.3 + 18.4]
76.6 4.3]
25.5 18.2]
52.0
Table 16: Gauss-ANM 50 nodes with hyperparameter search
Graph Type
GraN-DAD
1131.7
26.7
NG.5 36.6
1127
288.3
DMG-GNN
0700 7.8"
298.1 11807
100.2
1972 11066
106.3 48.9
148.2 11.6
4T3-4 105.07
NOTEARS
300.5
213.0
722.7 177.0
CAM
4.0
21.1 22.1
20.9
1225.9
15.3"
156.0
1200
121.3
1120.8
35.3
[M9.7
75.
3310.8
45.31
13611 156.7
131.9 30.71
301.9
GES
115011 9.8
1112.
1006.1
11364.5 8 91.8**
110177
8302 211.1
1115.
[1000.7 $10.9
132.0 90.0
11618
330.0 215.4]
1091.0 345.1]
Table 17: Gauss-ANM 20 nodes with hyperparameter search
Graph Type
SF4
Metrios
SIID
SID
SUDD
SID
SHDD
SUD
SID
Method
GrIN-DGG
2.3"
9.6 10.38
35.9 11.8*
120.4 +
2.4
6.55"
86.6 411
64.8 14.3
DAG-GNN
21.0
96.8 42.2*
77.2 17.22:66. 6.5
315.6
18.6
19.1 0.7
55.0 10-1-20.. 20.
30.2 5.4"
118.7 33.2
NOTEARS
101.5 30.6
100.4 1 47.0*
121.0 16.1
26700 06.5*
55.0 28.2
$7.6 $.2.1
60.7 8.3
154. 41.0
CAM
10.2"
64.5 29.3
121.7
73.1*
5.61
19.3 7
36.0 1.1
66.3 28.6
GSF
11.5
26.4
41.8
128
[221
14.11
42.4
[68.9
33.8)
33.67
13.31120.0)
57.44
21

---

Published as a conference paper at ICLR 2020
Table 18: Gauss-ANM- 100 nodes with hyperparameter search
Goph type
BR1
IR4
SII
GuN-DDG
15.4 7.5"
06.1 53.2"
1918 2 17.8
6090. 4 418.0
5165 10.2
210.0 $1.9
258.7 21.1
790.6
150.7
DAG-GNN
103.0
= 216.0
387.1 25.5
774.. 922.5
100.5::22 8.2
301.7 66.0
514.5 2 16.5
1257.8 183.2'
NOTHARS
E2:
107
B33.T 333-7
631-1
6572.4 444
2:4.3 61.
815.6 346.5
187.3
11411
1936-7 383
CAM
12.3
125. 66.3
198.8 22.2
1002.2 520.7
611 994
233.0 62.3
265.7
22.2
=51.4 200.0
GSF
100.2 9.5c
[10.8 242
387.6 23.9
(7585.] 686.2...
67.1 1402
254.5 53
$15.1 Taee
[1214,0 156
[]2.0]
70.4
19644
U18827=
$64.1
1988.2
1073.7
Table 19: PNL-GP with hyperparameter search
NNodes Type
50
EX1
SHD
122+
56.3
1107
DAG GNN
10.6 1.9"
35.8
19.8
1. 2.0
82.2 5.7.7
1.44"
330.4
69.9
152.5
10.2
2779.5
NOTEARS
206 11.4
30.5
30.5-1855 18.5*
6.5
68.4 6.09
102.1 27.3
200.8 65.8*
660.0
266.2
1746.0
g++
CAM
2.7 4.0
6..8118""
8.7
30.9 20.4
4.0 2.4"
10.7 12.4
52.3 8.5
913.9 1: 209
GSF
129.1.99
[1.55
8.7.
19.21...827
4.99
[281.6 70.7."
199.9
18.2
1878.0 122.4..
68.6=255
1.39.)
Table 20: PNL-MULT with hyperparameter search
#Nodes Type
GraN DAG
10.0 4.3*
29.1 9.7 9.7*
32.9 3.3
70.7 4.17
39.8 =: 28.2
213.6 97.3*
272.1 09.4
2021.6 A: 185.
DAG-GNN GNN
14.6 1161..11 5.1
36.9 10.6
3..91.1
85.8 1.2
643 27.8*
508.8 $17.2
212.5 12.3
2216.9 L. 05.6*
NOTEARS
28.8 9.1
30.3 11.8
36.4 3.8
784
160.2 67.6
463.5 205.1
229.2 =: 25.4
2108.8 =: 70.3
CAM
17.2
14.4
6.5
2.268.2
251.0 20.9
2026.2 18.2*
15.6
6..0
76.0+99...
06.4
12.ouu
5.3)
257.0
Table 21: LIN with hyperparameter search
Graph Type
Manics
SHD
SID
Method
GraN DAG
10.1 3.9
28.7 2 14.7
34.7 =
79.5
40.8 10.3
236.3
101.7
256.9 9: 55.7
2151.4 144.1
DAG GNN
9.0 2.7
35.6 11.4
1966 4.6*
63.9 1.5
6.8*
381.7 1 145.6
149 17.2
2070.7 1. 51.1
NOTEARS
14.0 =
7.9.
20.7 0.1
8.0
87.7 44.3.
294.3 99.3
20013 + 67.1
1772.7 #: 143.711
CAM
8.8 6.00"
15.5
33.0 2.8
77.1+44"
7
221.2 96.3
202.2 4: 14.3
1900.8 07.
GSF
10.7 3.5
1.44
33.4
[1.7 11.5"
SL4.66.
[158.1 115.9
195.5 9.9.
7133+
660.9 220.7]
Table 22: ADD-FUNC with hyperparameter search
#Nodes
Graph Type
Metrics
SHD
SHD
SHD
Method
GraN DAG
26-24" 2.4"
4.3 43"
31
37.1 12.4
1877
72.1 65.2*
20.1
1241.7 = 289.8*
DAG GNN
8.7 1 2.8
22.1
25
8.6
46.7 9.7*+
306.9 114.7*
194.0 20.4*
1940.3 1. 107.1*
NOTEARS
21.2 11.5.
155+999
13
41.3 11.6
186.8 $3.0
276.9 92.1
718.4 170.4
1105.1 250.1
3.=222"
8.1.6
622
28.5
10.0 4.6"
44.2 32.1
46.6 9.5
5.5+4..1
17.++123 12.3*
1..1
(4455
42.82"
140.4
[1674.4 133.9...
16.9
16.5)
48.7
145.2)
22

---

Published as a conference paper at ICLR 2020
Table 23:
Results for real and pseudo real data with hyperparameter search
Data Type
Protein signaling data set
SynTReN
20 nodes
Metrics
SIID
SIID-C
SID
SHID
SHID-C
SID
Method
GraN-DAG
12.0
9.0*
48.0
43.7 = 8.3
144.3 61.3
GraN-DAG+
14.0
11.0
57.0
46.9 14.9
49.5 = 14.7
1
DAG-GNN
16.0
14.0
59.0
3 5.0++
3 5.6
191.2 50.2
NOTEARS
15.0+
14.0*
58.0
27.5
45.8 27.7+
183.1
CAM
11.0*
9.0
31.0*
101.7
105.6+361
1
GSF
20.0
14.0
[37.0*
27.8 + 5.44
27.8 + 5.47
227..6+
60.0)
209.6 59.1
GES
47.07
50.0
[37.0
167.5 5.6
172.2 7.0"
[75.3 2
47.0]
97.6 = 30.8)
Table 24: Hyperparameter search spaces for each algorithm
Hyperparameter space
Log(lcarning rate) ~ - (first subproblemm)
Log(learning rate) ~ U -3, -4] (other subproblems)
e ~ (110- , 10-5)
Log(pruning cutoff) ~ Uf [f-5, - -2
GraN-DAG
# hidden units ~ Uf4, 8, 12]
# hidden layers ~ (11,, 2,3)
Constraint convergence tolerance U110 10- 10
8,100,10-10 10-10}
PNS threshold ~ 0.7. 1,2
0.75,
Log(learning rate) U1-4,
# hidden units in encoder ~ Uf16, 32, 64, 128, 256)
DAG-GNN
# hidden units decoder ~ U(16, 32, 64, 128, 256)
Bottleneck dimension (dimension of Z) U[1,5, 10, 100)
Constraint convergence tolerance ~ U { 6,000-,,0000 1001,00-100
LI regularizer coefficient (00.001, 0.005, 0.01, 0.5,
NOTEARS
Final threshold ~ U (0.0011, 0.003, 0.01, 0.00, 0.03, 0.1,0,11
Constraint convergence tolerance ~- U(10
10-8, 10
CAM
Log(Pruning cutoff) -
GSF
Log(RKHS regression regularizer) 0
1.4
GES
Log(Regularizer coefficient) ~ 44
23


Reviews:
[{"review_id": "HyefuvE-5B", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "Summary: \nThe authors propose a prediction model for directed acyclic graphs (DAGs) over a fixed set of vertices based on a neural network. The present work follows the previous work on undirected acyclic graphs, where the key constraint is (3), ensuring the acyclic property. The proposed method performed favorably on artificial/real data compared to previous baselines. \n\nComments: \nI do not understand yet if the proposed formulation really ensures the acyclic condition. More precisely, the condition (3) ensures undirected acyclic property, which also implies the directed ones. However, I am afraid that the condition (3) might also eliminate DAGs containing “undirected” cycles, i.e., cycles when we neglect directions of edges as well. So, I think a formal proof is necessary to show that the proposed formulation can output all possible DAGs, not a subset. \n\nThe present work heavily relies on the previous work of Zheng et al. (2018). The proposed formulation is designed for DAGs, but due to my concern about the correctness raised above, it is difficult for me to evaluate the originality of the present work. \n\nMinor comments:\n-“X_\\pi_j^G denote the random vector containing the variables corresponding to the parents of j in G.” – this definition is not clear. Please clarify the meaning. \n\n\nComments after Rebuttal:\nAfter reading the rebutall comments, I modified my score.", "extra_scores": {"experience_assessment": "I have published in this field for several years.", "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": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HkxS9xbTKH", "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": "################## Response to rebuttal ##################\n\nI would like to thank the authors for their diligent rebuttal. Overall, I consider that the issues raised in my review have been mostly addressed. \n\nIn summary, after reading the other reviews and going through the revised manuscript, I believe the paper would be a meaningful contribution to the field, both in terms of the potential applicability of GraN-DAG to certain problems and its potential to lead to future work. As a result, I have updated my score to support acceptance without reservations.\n\nIn particular, the authors now discuss clearly some of the potential limitations of the approach, such as the scaling of its asymptotic complexity with respect to the number of features or its statistical performance in the small sample size regime. The additional experimental results are, in my opinion, also a good addition to the manuscript, further clarifying the strengths and weaknesses of NOTEARS, CAM and GraN-DAG under different regimes and modelling assumptions. I believe these changes substantially improve the paper, as it will both help practitioners decide when to use GraN-DAG as well as clearly delineate avenues for future work.\n\n# Minor points / additional questions\n\nThe new results regarding convergence in Table 4 suggest that early stopping might indeed be instrumental. I wonder the extent to which NOTEARS could benefit from early stopping, be it in terms of runtime or, less likely due to its use of explicit regularization, the quality of the structures it retrieves. \n\n# Typos\n\nSection 3.6: … GraN-DAG performs fewer iteration(s) than NOTEARS…\nSection 5: GraN-DAGand DAG-GNN -was- were not designed with…\nAppendix A.1: … graphs can be roughly -halfed- halved when executed…\nAppendix A.5: … graphs were generated using the -BarabsiAlbert- Barabási-Albert model …\nAppendix A.9: … to selecting hyperparameters in practice consist(s) in ...\n\n################## Original review ##################\n\n# High-level assessment\n\nThe manuscript presents an alternative to the method in [2] to extend the ideas of [1] to nonlinear DGMs using neural networks. \n\nFrom a methodological perspective, in my opinion the proposed approach is sound and sufficiently novel to be a significant contribution. While some limitations might be present, such as being prone to overfitting without relying on additional (pre/post)-processing heuristics or perhaps being computationally intensive for high-dimensional data ($O(p^{3})$ runtime per SGD step), I believe that addressing those issues (if indeed present) could be left for future work.\n\nThe paper is written clearly, being easy-to-follow while providing enough detail to understand the proposed method in depth and attempt to reproduce the results.\n\nIn my opinion, perhaps the weakest aspect of the paper concerns the experimental analysis, which is somewhat limited in its scope. Most notably, I believe the current experimental setup focuses on cases likely to be favorable to the proposed approach relative to the baselines ($n > p$, nonlinear and non-additive ground-truth models) while presenting no results for other settings commonly encountered in structure learning applications, such as high-dimensional $n << p$ problems or ground-truth models given by linear or generalized linear SCMs.\n\nTaking all those factors into consideration, I currently lean towards a “weak accept” rating, but would gladly modify my score after the author discussion phase as appropriate.\n\n# Summary\n\nIn this paper, the authors propose a novel approach for learning the structure of Directed Graphical Models (DGMs) from observational data that allows to flexibly model nonlinear relationships between variables using neural networks.\n\nThe proposed method builds on the recent breakthrough in [1], which pioneered the use of a smooth characterization of acyclicity for directed graphs. In particular, the authors in [1] use this characterization to learn the structure of **linear** DGMs by solving a (nonconvex) continuous, constrained optimization problem.\n\nThe main contribution of this manuscript is an extension of the approach in [1] to the nonlinear setting by making use of neural networks, and which appears to outperform the only existing method that shares a similar high-level goal [2], presumably by allowing more flexibility in its parametrization.\n\nIn a nutshell, the authors accomplish this by:\n\n1. Flexibly modelling the conditional distribution of each variable $X_{j}$ as a nonlinear function of (possibly) all other variables $\\bf{X_{-j}}$ using an MLP, i.e. $p_{j}(X_{j} \\mid \\bf{X_{-j}}) = f_{j}(X_{j} ; NN_{j}(\\bf{X_{-j}}))$, with $f_{j}$ denoting any probability density differentiable w.r.t. its parameters and $NN_{j}$ being a neural network with parameters $\\theta_{j}$.\n\n2. Deriving a function $C_{i,j}(\\theta_{j})$ such that $C_{i,j}(\\theta_{j}) = 0$ is a sufficient condition for the conditional distribution of $X_{j}$ to be independent of $X_{i}$.\n\n3. Using $C_{i,j}(\\theta_{j})$, which is differentiable w.r.t. $\\theta_{j}$ for most values of $\\theta_{j}$, to define a “connectivity” matrix to be directly plugged into the acyclicity criterion in [1].\n\n4. Following [1], the model is then fit by (constrained) maximum likelihood using an augmented Lagrangian method. However, unlike [1], the unconstrained subproblem in each iteration of the augmented Lagrangian method is solved using stochastic gradient descent, rather than a batch quasi-Newton method, and regularized using early stopping on an external validation set, rather than by using explicit $L_{1}$ regularization to induce sparsity.\n\n5. The authors also incorporate additional “heuristics”, such as borrowing the Preliminary Neighbour Selection (PNS) and feature selection pruning techniques from [3], which are reported to greatly enhance the performance of the proposed approach (Table 4), and adapt the thresholding scheme in [1] to the nonlinear setting.\n\nThe proposed method is evaluated using synthetic and semi-synthetic data, as well as one real-world dataset.  The results suggest that, in low dimensional problems ($n > p$, $p \\le 100$) for which the ground-truth relationships between variables are nonlinear and non-additive, the proposed approach clearly outperforms [2] and also outperforms [3] in many settings, albeit often not by a statistically significant margin in the latter case.\n\n# Major points / suggestions\n\n1. I believe it could be argued that one of the main drawbacks of the smooth acyclicity constraint in [1] is its high computational complexity w.r.t. the number of features. However, [1] uses a batch quasi-Newton method to solve the unconstrained subproblems in their optimization routine, which should in principle lead to relatively few $O(p^3)$ matrix exponential evaluations being necessary. \n\nAs it relies on the same characterization of acyclicity as [1], this drawback might also limit the scalability of the proposed approach for high-dimensional data. Moreover, since the method in this manuscript uses mini-batch gradient descent instead, this problem could in fact be exacerbated by perhaps requiring substantially more $O(p^3)$ matrix exponential evaluations.\n\nTo this end, I would be glad if the authors could clarify whether this could indeed be a limitation or not. In any case, I believe it could also be helpful to include a small discussion regarding computational complexity for the proposed approach and related work in the manuscript or appendix.\n\n2. All experimental results reported in the manuscript use only data generation models with nonlinear and non-additive relationships between variables.\n\nThis setup is well aligned with the contribution, as it emphasises its ability to flexibly model nonlinear relationships between variables learnt using global, gradient-based optimization. In contrast, most baselines are limited to either linear/generalized linear models or use some variant of greedy search.\n\nWhile investigating this setting is certainly the main priority given the problem statement, I believe it would also be important to characterize the extent to which the performance of the proposed approach deteriorates when this extra flexibility is not needed, if at all. To this end, I would encourage the authors to repeat their experiments for synthetic data, but generated under linear and generalized linear SCMs.\n\nIn particular, it is remarkable that CAM [3], which suffers from model misspecification due assuming the data is generated by a generalized linear model, remains so competitive throughout all experiments shown in the current version of the manuscript. It would be interesting to see if the proposed approach can match or still outperform CAM under CAM’s “ideal scenario”.\n\n3. Parallel to the previous point, all experimental results in the manuscript so far concern somewhat low-dimensional data $p \\le 100$ and relatively abundant data $n \\approx 1,000 > p$.\n\nThis setup might also be well aligned with the proposed approach, which could be more prone to overfitting than some of the baselines due to its high capacity and its lack of explicit sparsity-inducing regularization, which was found to boost performance for NOTEARS [1, Fig. 3(b)]. \n\nAlso, in the small sample size regime, the reliance on an external validation set for early stopping might become a liability by (i) reducing the amount of data available for model fitting and (ii) making the fitted DGMs less stable w.r.t. changes in the random seed which, compounded with the nonconvexity of the objective, could be problematic if the DGM is to be interpreted for downstream tasks (e.g. in computational biology).\n\nTo this end, I would likewise encourage the authors to extend their experimental setup for synthetic data by (i) repeating their experiments for one “high-dimensional” setting ($n << p$) and (ii) studying the stability of the DAG w.r.t. different training runs.\n\n4. All experimental results currently focus heavily on accuracy, while the computational aspect of the problem is largely secondary.\n\nTo this end, it would be interesting to study, using synthetic data, the scalability of the proposed approach relative to the best performing baseline approaches with respect to (i) number of features, (ii) sparsity of the ground-truth DGM and (iii) sample size, covering a sufficiently broad range in each case (low-dimensional vs high-dimensional, sparse vs dense, scarce vs abundant sample size).\n\n# Minor points / suggestions\n\n1. Perhaps the decision to permanently set the thresholding masks to 0 when an entry of the connectivity matrix is “small enough” could be justified further. For example, by showing that this heuristic rarely “zeros out” incorrectly any entries, or that such entries would seldom become nonzero again if the model was trained for longer without being masked.\n\n2. I might have missed it, but I could not find the info about how the train/validation split is formed for the experiments following Section A.1.\n\n# References\n\n[1] Zheng, Xun, et al. \"DAGs with no tears: Continuous optimization for structure learning.\" *Advances in Neural Information Processing Systems*. 2018.\n[2] Yu, Yue, et al. \"DAG-GNN: DAG Structure Learning with Graph Neural Networks.\" *International Conference on Machine Learning*. 2019.\n[3] Bühlmann, Peter, Jonas Peters, and Jan Ernest. \"CAM: Causal additive models, high-dimensional order search and penalized regression.\" *The Annals of Statistics* 42.6 (2014): 2526-2556.", "extra_scores": {"experience_assessment": "I have read many 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": {}}, {"review_id": "HJxLg1uTYS", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "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 problem of learning the structure of directed acyclic graphs in the presence of nonlinearities. The proposed approach is an extension of the NOTEARS algorithm which uses a neural network for each node in the graph during structure learning. This adaptation allows for non-linear relationships to be easily modeled. In addition to the proposed adaptation, the authors employ a number of heuristics from the causal discovery literature to improve the efficacy of the search. Empirical results are provided which compare the proposed algorithm to prior art. \n\nOverall, I found the paper to be well written and sensible. However, there are a few items that prevent me from recommending acceptance at this time: \n\n1) The proposed approach seems to necessitate a fair number of hyperparameters. How were these chosen? On what basis should practitioners choose their hyperparameters for real-world application?\n\n2) The authors employ causal language in both the introduction and the experiment section. However, there is no notion of (a) whether the proposed algorithm is sound and complete,  and (b) under what assumptions we can expect a fully directed graph to be reasonable (I assume the additive noise model?). \n\n3) The experiments seem to indicate that prior work (CAM) outperforms or performs similarly to the proposed method. While this shouldn't prevent a paper from being published, it would be nice to see an extended discussion as to why the authors think this is occuring. ", "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 thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: In this paper, the authors propose a novel approach for learning the structure of a directed acyclic graph from observational data that allows to flexibly model nonlinear relationships between variables using neural networks. While the reviewers initially had concerns with respect to the positioning of the paper and various questions regarding theoretical results and experiments, these concerns have been addressed satisfactorily during the discussion period.  The paper is now acceptable for publication in ICLR-2020. 
Decision: Accept (Poster)