GRADIENT-BASED NEURAL DAG LEARNING
Year: 2020

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.

heading: 1 INTRODUCTION
text: Structure learning and causal inference have many important applications in different areas of science such as genetics (Koller & Friedman, 2009; Peters et al., 2017), biology (Sachs et al., 2005) and economics (Pearl, 2009). Bayesian networks (BN), which encode conditional independencies using directed acyclic graphs (DAG), are powerful models which are both interpretable and computationally tractable. Causal graphical models (CGM) (Peters et al., 2017) are BNs which support interventional 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 learning 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 exponentially 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 a causal graphical model is to predict the effect of various interventions. A CGM can be best estimated when given interventional data, but interventions are often costly or impossible to obtained. As an alternative, one can use exclusively observational data and rely on different assumptions which make the graph identifiable from the distribution (see Section 2.2). This is the approach employed in this paper.
We propose a score-based method (Koller & Friedman, 2009) for structure learning named GraNDAG 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 a linear structural equation model (SEM) (Pearl, 2009) (see Section 2.3) and enforces acyclicity using a constraint which is both efficiently computable and easily differentiable, 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 use an argument similar to what is used in
†Canada CIFAR AI Chair Correspondence to: sebastien.lachapelle@umontreal.ca
Zheng et al. (2018) and apply it first at 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 families of conditional probability distributions, our experiments focus on the special case of nonlinear Gaussian additive noise models since, under specific assumptions, it provides appealing theoretical 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 et al., 2019), a recent nonlinear extension of Zheng et al. (2018) which uses an evidence lower bound as score.
Our second contribution is to 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 (Bühlmann et 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 an implementation of GraN-DAG here.

heading: 2 BACKGROUND
text: Before presenting GraN-DAG, we review concepts relevant to structure and causal learning.

heading: 2.1 CAUSAL GRAPHICAL MODELS
text: We suppose the natural phenomenon of interest can be described by a random vector X ∈ Rd entailed by an underlying CGM (PX ,G) where PX is a probability distribution over X and G = (V,E) is a DAG (Peters et al., 2017). Each node j ∈ V corresponds to exactly one variable in the system. Let πGj denote the set of parents of node j in G and let XπGj 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 to be Markov to G, i.e. we can write the probability density function (pdf) of PX as p(x) = ∏d j=1 pj(xj |xπGj ) where pj(xj |xπGj ) is the conditional pdf of variable Xj given XπGj . A CGM can be thought of as a BN in which directed edges are given a causal meaning, allowing it to answer queries regarding interventional distributions (Koller & Friedman, 2009).

heading: 2.2 STRUCTURE IDENTIFIABILITY
text: In general, it is impossible to recover G given only samples from PX , i.e. without interventional data. It is, however, customary to rely on a set of assumptions to render the structure fully or partially identifiable.
Definition 1 Given a set of assumptions A on a CGM M = (PX ,G), its graph G is said to be identifiable from PX if there exists no other CGM M̃ = (P̃X , G̃) satisfying all assumptions in A such that G̃ 6= G and P̃X = PX .
There are many examples of graph identifiability results for continuous variables (Peters et al., 2014; Peters & Bühlman, 2014; Shimizu et al., 2006; Zhang & Hyvärinen, 2009) as well as for discrete variables (Peters et al., 2011). These results are obtained by assuming that the conditional densities belong to a specific parametric family. For example, if one assumes that the distribution PX is entailed by a structural equation model of the form
Xj := fj(XπGj ) +Nj with Nj ∼ N (0, σ2j ) ∀j ∈ V (1)
where fj is a nonlinear function satisfying some mild regularity conditions and the noises Nj are mutually independent, then G is identifiable from PX (see Peters et al. (2014) for the complete theorem and its proof). This is a particular instance of additive noise models (ANM). We will make use 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 G itself, but the Markov equivalence class to which it belongs (Spirtes
et al., 2000). A Markov equivalence class is a 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 et al., 2017). Some algorithms we use as baselines in Section 4 output only a CPDAG.

heading: 2.3 NOTEARS: CONTINUOUS OPTIMIZATION FOR STRUCTURE LEARNING
text: Structure learning is the problem of learning G using a data set of n samples {x(1), ..., x(n)} from PX . Score-based approaches cast this problem as an optimization problem, i.e. Ĝ = arg maxG∈DAG S(G) where S(G) is a 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 5 for a review). We now present the work of Zheng et al. (2018) which proposes to cast this combinatorial optimization problem into a continuous constrained one.
To do so, the authors propose to encode the graph G on d nodes as a weighted adjacency matrix U = [u1| . . . |ud] ∈ Rd×d which represents (possibly negative) coefficients in a linear SEM of the form Xj := u>j X + Ni ∀j where Nj is a noise variable. Let GU be the directed graph associated with the SEM and let AU be the (binary) adjacency matrix associated with GU . One can see that the following equivalence holds:
(AU )ij = 0 ⇐⇒ Uij = 0 (2)
To make sure GU is acyclic, the authors propose the following constraint on U :
Tr eU U − d = 0 (3)
where eM , ∑∞ k=0 Mk k! is the matrix exponential and is the Hadamard product.
To see why this constraint characterizes acyclicity, first note that (AUk)jj is the number of cycles of length k passing through node j in graph GU . Clearly, for GU to be acyclic, we must have TrAU
k = 0 for k = 1, 2, ...,∞. By equivalence (2), this is true when Tr(U U)k = 0 for k = 1, 2, ...,∞ . From there, one can simply apply the definition of the matrix exponential to see why constraint (3) characterizes acyclicity (see Zheng et al. (2018) for the full development).
The authors propose to use a regularized negative least square score (maximum likelihood for a Gaussian noise model). The resulting continuous constrained problem is
max U S(U,X) , − 1 2n ‖X−XU‖2F − λ‖U‖1 s.t. Tr eU U − d = 0 (4)
where X ∈ Rn×d is the design matrix containing all n samples. The nature of the problem has been drastically changed: we went from a combinatorial to a continuous problem. The difficulties of combinatorial optimization have been replaced by those of non-convex optimization, since the feasible set is non-convex. Nevertheless, a 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 no need to 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 iteration. Continuous approaches to combinatorial optimization have sometimes demonstrated improved performance over discrete approaches in the literature (see for example Alayrac et al. (2018, §5.2) where they solve the multiple sequence alignment problem with a continuous optimization method).

heading: 3 GRAN-DAG: GRADIENT-BASED NEURAL DAG LEARNING
text: We propose a new nonlinear extension to the framework presented in Section 2.3. For each variable Xj , we learn a fully connected neural network with L hidden layers parametrized by φ(j) , {W (1) (j) , . . . ,W (L+1) (j) } where W (`) (j) is the `th weight matrix of the jth NN (biases are omitted for clarity). Each NN takes as input X−j ∈ Rd, i.e. the vector X with the jth component masked to zero, and outputs θ(j) ∈ Rm, the m-dimensional parameter vector of the desired distribution family
for variable Xj .1 The fully connected NNs have the following form
θ(j) ,W (L+1) (j) g(. . . g(W (2) (j) g(W (1) (j)X−j)) . . . ) ∀j (5)
where g is a nonlinearity applied element-wise. Note that the evaluation of all NNs can be parallelized on GPU. Distribution families need not be the same for each variable. Let φ , {φ(1), . . . , φ(d)} represents all parameters of all d NNs. Without any constraint on its parameter φ(j), neural network j models the conditional pdf pj(xj |x−j ;φ(j)). Note that the product∏d j=1 pj(xj |x−j ;φ(j)) does not integrate to one (i.e. it is not a joint pdf), since it does not decompose according to a DAG. We now show how one can constrain φ to make sure the product of all conditionals outputted by the NNs is a joint pdf. The idea is to define a new weighted adjacency matrixAφ similar to the one encountered in Section 2.3, which can be directly used inside the constraint of Equation 3 to enforce acyclicity.

heading: 3.1 NEURAL NETWORK CONNECTIVITY
text: Before defining the weighted adjacency matrix Aφ, we need to focus on how one can make some NN outputs unaffected by some inputs. Since we will discuss properties of a single NN, we drop the NN subscript (j) to improve readability.
We will use the term neural network path to refer to a computation path in a NN. For example, in a NN with two hidden layers, the sequence of weights (W (1)h1i,W (2) h2h1 ,W (3) kh2
) is a 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 (1)h1i||W (2) h2h1 ||W (3)kh2 | ≥ 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 to output k are inactive (i.e. the sum of their path products is zero), then output k does not depend on input i anymore since the information in input i will never reach output k. The sum of all path products from input i to output k for all input i and output k can be easily computed by taking the following matrix product.
C , |W (L+1)| . . . |W (2)||W (1)| ∈ Rm×d≥0 (6)
where |W | is the element-wise absolute value of W . Let us name C the neural network connectivity matrix. It can be verified that Cki is the sum of all NN path products from input i to output k. This means it is sufficient to have Cki = 0 to render output k independent of input i.
Remember that each NN in our model outputs a parameter vector θ for a conditional distribution and that we want the product of all conditionals to be a valid joint pdf, i.e. 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 θ not dependent on certain inputs of the NN. To have θ independent of variable Xi, it is sufficient to have ∑m k=1 Cki = 0.

heading: 3.2 A WEIGHTED ADJACENCY MATRIX
text: We now define a weighted adjacency matrix Aφ that can be used in constraint of Equation 3.
(Aφ)ij ,
{ ∑m k=1 ( C(j) ) ki , if j 6= i
0, otherwise (7)
where C(j) denotes the connectivity matrix of the NN associated with variable Xj .
As the notation suggests,Aφ ∈ Rd×d≥0 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 (Aφ)ij = 0 =⇒ θ(j) does not depend on variable Xi (8) We note Gφ to be the directed graph entailed by parameter φ. We can now write our adapted acyclicity constraint:
h(φ) , Tr eAφ − d = 0 (9) 1Not all parameter vectors need to have the same dimensionality, but to simplify the notation, we suppose
mj = m ∀j
Note that we can compute the gradient of h(φ) w.r.t. φ (except at points of non-differentiability arising from the absolute value function, similar to standard neural networks with ReLU activations (Glorot et al., 2011); these points did not appear problematic in our experiments using SGD).

heading: 3.3 A DIFFERENTIABLE SCORE AND ITS OPTIMIZATION
text: We propose solving the maximum likelihood optimization problem
max φ EX∼PX d∑ j=1 log pj(Xj |Xπφj ;φ(j)) s.t. Tr e Aφ − d = 0 (10)
where πφj denotes the set of parents of node j in graph Gφ. Note that ∑d j=1 log pj(Xj |Xπφj ;φ(j)) is a valid log-likelihood function when constraint (9) is satisfied.
As suggested in Zheng et al. (2018), we apply an augmented Lagrangian approach to get an approximate solution to program (10). Augmented Lagrangian methods consist of optimizing a sequence of subproblems for which the exact solutions are known to converge to a stationary point of the constrained problem under some regularity conditions (Bertsekas, 1999). In our case, each subproblem is
max φ L(φ, λt, µt) , EX∼PX d∑ j=1 log pj(Xj |Xπφj ;φ(j))− λth(φ)− µt 2 h(φ)2 (11)
where λt and µt are 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 a well-known stochastic gradient algorithm popular for NN in part for its implicit regularizing effect (Poggio et al., 2018). See Appendix A.1 for details regarding the optimization procedure.
In the current section, we presented GraN-DAG in a general manner without specifying explicitly which distribution family is parameterized by θ(j). In principle, any distribution family could be employed as long as its log-likelihood can be computed and differentiated with respect to its parameter θ. However, it 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 ,G).
Proposition 1 Let φ∗ and Gφ∗ be the optimal solution to (10) and its corresponding graph, respectively. Let M(A) be the set of CGM (P ′,G′) 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 ,G) ∈ C and C = M(A) for a specific set of assumptions A such that G is identifiable from PX , then Gφ∗ = G.
Proof: Let Pφ be the joint distribution entailed by parameter φ. Note that the population loglikelihood EX∼PX log pφ(X) is maximal iff Pφ = PX . We know this maximum can be achieved by a specific parameter φ∗ since by hypothesis (PX ,G) ∈ C. Since G is identifiable from PX , we know there exists no other CGM (P̃X , G̃) ∈ C such that G̃ 6= G and P̃X = PX . Hence Gφ∗ has to be equal to G. In Section 4.1, we empirically explore the identifiable setting of nonlinear Gaussian ANMs introduced in Section 2.2. In practice, one should keep in mind that solving (10) exactly is hard since the problem is non-convex (the augmented Lagrangian converges only to a stationary point) and moreover we only have access to the empirical log-likelihood (Proposition 1 holds for the population case).

heading: 3.4 THRESHOLDING
text: The solution outputted by the augmented Lagrangian will satisfy the constraint only up to numerical precision, thus several entries of Aφ might not be exactly zero and require thresholding. To do so, we mask the inputs of each NN j using a binary matrix M(j) ∈ {0, 1}d×d initialized to have (M(j))ii = 1 ∀i 6= j and zeros everywhere else. Having (M(j))ii = 0 means the input i of NN
j has been thresholded. This mask is integrated in the product of Equation 6 by doing C(j) , |W (L+1)(j) | . . . |W (1) (j) |M(j) without changing the interpretation of C(j) (M(j) can be seen simply as an extra layer in the NN). During optimization, if the entry (Aφ)ij is smaller than the threshold = 10−4, the corresponding mask entry (M(j))ii is set to zero, permanently. The masks M(j) are never updated via gradient descent. We also add an iterative thresholding step at the end to ensure the estimated graph Gφ is acyclic (described in Appendix A.2).

heading: 3.5 OVERFITTING
text: 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 a held-out data set and declare convergence once it has stopped improving. This approach is known as early stopping (Prechelt, 1997). Second, to optimize (11), we use a stochastic gradient algorithm variant which is now known to have an implicit regularizing effect (Poggio et 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 (Bühlmann et al., 2014) to 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) before running the optimization procedure as was also recommended in Bü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 L1 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-GNN with PNS and pruning, without noting a significant improvement.

heading: 3.6 COMPUTATIONAL COMPLEXITY
text: 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 a standard choice when it comes to NN training because it scales well with both the sample size and the number of parameters and it implicitly regularizes learning. Similarly to NOTEARS, GraN-DAG requires the evaluation of the matrix exponential of Aφ at each iteration costing O(d3). NOTEARS justifies the use of a batch proximal quasi-Newton algorithm by the low number of O(d3) iterations required to converge. Since GraNDAG 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.1). We hypothesize this is due to early stopping which avoids having to wait until the full convergence of the subproblems hence limiting the total number of iterations. Moreover, for the graph sizes considered in this paper (d ≤ 100), the evaluation of h(φ) 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 approximate gradient of the log-likelihood (costing O(d2) assuming a fixed minibatch size, NN depth and width) appears to be of greater importance for d ≤ 100.

heading: 4 EXPERIMENTS
text: In this section, we compare GraN-DAG to various baselines 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 (Bühlmann et al., 2014), GSF (Huang et al., 2018), GES (Chickering, 2003) and PC (Spirtes et 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 µ̂(j), i.e. θ(j) := µ̂(j) and Xj |XπGj ∼ N (µ̂(j), σ̂ 2 (j)) ∀j. The parameters σ̂ 2 (j) are learned as well, but do not depend on
the parent variables XπGj (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 Erdős-Rényi (ER) scheme described in Appendix A.5 (denoted by RANDOM). For each approach, we evaluate the estimated graph on two metrics: the structural hamming distance (SHD) and the structural interventional distance (SID) (Peters & Bühlmann, 2013). The former simply counts the number of missing, falsely detected or 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(xj |do(Xi = x̄)) would be miscalculated if we were to use the estimated graph to form the parent adjustement set. Note that GSF, GES and PC output only a CPDAG, hence the need to report a 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 a procedure on almost every settings discussed in the present section.

heading: 4.1 SYNTHETIC DATA
text: We have generated different data set types which vary along four dimensions: data generating process, number of nodes, level of edge sparsity and graph type. We consider two graph sampling schemes: Erdős-Rényi (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 randomly sampled following either the ER or the SF scheme. Then, the data is generated according to a specific sampling scheme.
The first data generating process we consider is the nonlinear Gaussian ANM (Gauss-ANM) introduced in Section 2.2 in which data is sampled following Xj := fj(XπGj ) + Nj with mutually independent noises Nj ∼ N (0, σ2j ) ∀j where the functions fj are independently sampled from a Gaussian process with a unit bandwidth RBF kernel and with σ2j sampled uniformly. As mentioned in Section 2.2, we know G to be identifiable from the distribution. Proposition 1 indicates that the modelling choice of GraN-DAG together with this synthetic data ensure that solving (10) to optimality 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 fj 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 ER1 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.
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 performs 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 a strong form of parameter sharing between the fj 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. fj(xπGj ) = ∑ i∈πGj
fij(xj) ∀j, it manages to compete with GraN-DAG which does not make this incorrect modelling assumption2. 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 to NOTEARS) and nonlinear Gaussian data with additive functions (better adjusted to CAM) named LIN and ADD-FUNC, respectively. See Appendix A.5 for the details of their generation. We report the results of GraN-DAG and other baselines 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 ER1 graph). However, GraN-DAG outperforms all other methods which can be explained by the fact that the conditions of Proposition 1 are satisfied (supposing the functions∑ i∈πGj fij(Xi) 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 & Hyvärinen, 2009), PNL-GP and PNL-MULT (see Appendix A.5 for details about their generation). A post nonlinear model has the form Xj := gj(fj(XπGj ) + Nj) where Nj is a noise variable. Note that GraN-DAG (with the current modelling choice) and CAM do not have the representational power to 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 gj , 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 PNL-GP but not on PNL-MULT where it performs better in SID.

heading: 4.2 REAL AND PSEUDO-REAL DATA
text: 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 dataset 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 Bulcke, 2006). This generator was 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.
2Although it is true that GraN-DAG does not wrongly assume that the functions fj are additive, it is not clear whether its neural networks can exactly represent functions sampled from the Gaussian process.
In applications, it is not clear whether the conditions of Proposition 1 hold since we do not know whether (PX ,G) ∈ 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 σ̂2(i) to depend on the parent variables XπGi through the NN, i.e. θ(i) := (µ̂(i), log σ̂ 2 (i)). Note that this is violating the additive noise assumption (in ANMs, the noise is independent of the parent variables).
In addition to metrics used in Section 4.1, 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 in Table 3.
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 on par in terms of SHD and SHD-C. On this data set, DAG-GNN has a reasonable performance, beating GraN-DAG in SID, but not in SHD. On SynTReN, GraN-DAG obtains the best SHD but not the best SID. Overall, GraN-DAG is always competitive with the best methods of each task.

heading: 5 RELATED WORK
text: 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.
Score-based methods (Koller & Friedman, 2009; Peters et al., 2017) cast the problem of structure learning as an 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 information 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 G fully identifiable. For example, Peters & Bühlman (2014) relies on a linear Gaussian model with equal variances to render the DAG identifiable. RESIT (Peters et al., 2014), assumes nonlinear relationships with additive Gaussian noise and greedily maximizes an independence-based score. However, RESIT does not scale well to graph of more than 20 nodes. CAM (Bühlmann et 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 in 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 is a drastic departure from previous combinatorial approaches which enables the application of well studied numerical solvers for continuous optimizations. Recently, Yu et al. (2019) proposed DAG-GNN, a graph neural network architecture (GNN) which can be used to learn DAGs via the maximization of an evidence lower bound. By design, a 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
learning to support nonlinear relationships. Although Yu et al. (2019) provides empirical comparisons 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-based approaches which uses integer linear programming (ILP) (Jaakkola et al., 2010; Cussens, 2011) which internally solve continuous linear relaxations. Connections between such methods and the continuous constrained approaches are yet to 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-dimensional setting in which the number of samples is significantly smaller than the number of nodes. Bühlmann 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 a 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 et al., 2018), when used for multivariate data, requires an initial skeleton to learn the different functional relationships.
GraN-DAG has strong connections with MADE (Germain et al., 2015), a method used to learn distributions using a masked NN which enforces the so-called autoregressive property. The autoregressive 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, it would be similar to MADE, except the variable ordering would be learned from data instead of fixed a priori.

heading: 6 CONCLUSION
text: The continuous constrained approach to structure learning has the advantage of being more global than other approximate greedy methods (since it 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 score-based approach for structure learning supporting nonlinear relationships while leveraging 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 is 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.

heading: ACKNOWLEDGMENTS
text: 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 Rémi Le Priol, Tatjana Chavdarova, Charles Guille-Escuret, Nicolas Gagné and Yoshua Bengio for insightful discussions as well as Alexandre Drouin and Florian Bordes for technical support. The experiments were in part enabled by computational resources provided by Calcul Québec, Compute Canada and Element AI.

heading: A APPENDIX

heading: A.1 OPTIMIZATION
text: Let us recall the augmented Lagrangian:
max φ L(φ, λt, µt) , EX∼PX d∑ i=1 log pi(Xi|Xπφi ;φ(i))− λth(φ)− µt 2 h(φ)2 (12)
where λt and µt are the Lagrangian and penalty coefficients of the tth subproblem, respectively. In all our experiments, we initialize those coefficients using λ0 = 0 and µ0 = 10−3. We approximately solve each non-convex subproblem using RMSprop (Tieleman & Hinton, 2012), a stochastic gradient descent variant popular for NNs. We use the following gradient estimate:
∇φL(φ, λt, µt) ≈ ∇φL̂B(φ, λt, µt)
with L̂B(φ, λt, µt) , 1 |B| ∑ x∈B d∑ i=1 log pi(xi|xπφi ;φ(i))− λth(φ)− µt 2 h(φ)2
(13)
where B is a minibatch sampled from the data set and |B| is the minibatch size. The gradient estimate ∇φL̂B(φ, λt, µt) can be computed using standard deep learning libraries. We consider a subproblem has converged when L̂H(φ, λt, µt) evaluated on a held-out data set H stops increasing. Let φ∗t be the approximate solution to subproblem t. Then, λt and µt are updated according to the following rule:
λt+1 ← λt + µth (φ∗t ) µt+1 ← { ηµt, if h (φ∗t ) > γh ( φ∗t−1 ) µt, otherwise
(14)
with η = 10 and γ = 0.9. Each subproblem t is initialized using the previous subproblem solution φ∗t−1. The augmented Lagrangian method stops when h(φ) ≤ 10−8. Total number of iterations before augmented Lagrangian converges: In GraN-DAG and NOTEARS, every subproblem is approximately solved using an iterative algorithm. Let T be the number of subproblems solved before the convergence of the augmented Lagrangian. For a given subproblem t, let Kt be the number of iterations executed to approximately solve it. Let I = ∑T t=1Kt be the total number of iterations before the augmented Lagrangian converges. Table 4 reports the total number of iterations I for GraN-DAG and NOTEARS, averaged over ten data sets. Note that the matrix exponential is evaluated once per iteration. Even though GraN-DAG uses a stochastic gradient algorithm, it requires less iterations than NOTEARS which uses a batch proximal quasi-Newton method. We hypothesize early stopping avoids having to wait until full convergence before moving to 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 to evaluate (total runtime ≈ 10 minutes against ≈ 1 minute for 20 nodes graphs and ≈ 4 hours against ≈ 1 hour for 100 nodes graphs). GraN-DAG runtime on 100 nodes graphs can be roughly halved when executed on GPU.

heading: A.2 THRESHOLDING TO ENSURE ACYCLICITY
text: The augmented Lagrangian outputs φ∗T where T is the number of subproblems solved before declaring convergence. Note that the weighted adjacency matrix Aφ∗T will most likely not represent an acyclic graph, even if we threshold as we learn, as explained in Section 3.4. We need to remove
additional edges to obtain a 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 (Aφ∗T )ij up to the edge with the highest (Aφ∗T )ij . This amounts to gradually increasing the threshold untilAφ∗T is acyclic. However, this approach has the following flaw: It is possible to have (Aφ∗T )ij significantly higher than zero while having θ(j) almost completely independent of variable Xi. This can happen for at least two reasons. First, the NN paths from input i to output k might end up cancelling each others, rendering the input i inactive. Second, some neurons of the NNs might always be saturated for the observed range of inputs, rendering some NN paths effectively inactive without being inactive in the sense described in Section 3.1. Those two observations illustrate the fact that having (Aφ∗T )ij = 0 is only a sufficient condition to have θ(j) independent of variable Xi and not a necessary one.
To avoid this issue, we consider the following alternative. Consider the functionL : Rd 7→ Rd which maps all d variables to their respective conditional likelihoods, i.e. Li(X) , pi(Xi | X
π φ∗ T i ) ∀i. We consider the following expected Jacobian matrix
J , EX∼PX ∣∣∣∣ ∂L∂X ∣∣∣∣> (15) where ∣∣ ∂L ∂X
∣∣ is the Jacobian matrix of L evaluated at X , in absolute value (element-wise). Similarly to (Aφ∗T )ij , the entry Jij can be loosely interpreted as the strength of edge (i, j). We propose removing edges starting from the lowestJij to the highest, stopping as soon as acyclicity is achieved. We believe J is better than Aφ∗T 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 Aφ∗T yields better results, and thus we report the results with J in this paper.

heading: A.3 PRELIMINARY NEIGHBORHOOD SELECTION AND DAG PRUNING
text: PNS: For graphs of 50 nodes or more, GraN-DAG performs a preliminary neighborhood selection (PNS) similar to what has been proposed in Bü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 fitting an extremely randomized trees (Geurts et al., 2006) (using ExtraTreesRegressor from scikit-learn) for each variable against all the other variables. For each node a feature importance score based on the gain of purity is calculated. Only nodes that have a feature importance score higher than 0.75 · mean are kept as potential parent, where mean is the mean of the feature importance 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.
Pruning: Once the thresholding is performed and a DAG is obtained as described in A.2, GraNDAG performs a pruning step identical to CAM (Bühlmann et al., 2014) in order to remove spurious edges. We use the implementation of Bühlmann et al. (2014) based on the R function gamboost from the mboost package. For each variable Xi, a generalized additive model is fitted against the current parents of Xi 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 Bühlmann et al. (2014) observed, this pruning phase generally has the effect of greatly reducing the SHD without considerably changing the SID.
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 it affects the likelihood of the model (hence the need to retrain).
A first observation is that adding PNS and pruning improve the NLL on the validation set while deteriorating the NLL on the training set, showing that those two steps are indeed reducing overfitting. Secondly, the effect on SHD is really important while the effect on SID is almost nonexistent. This can be explained by the fact that SID has more to do with the ordering of the nodes than with false positive edges. For instance, if we have a complete DAG with a node ordering coherent with the ground truth graph, the SID is zero, but the SHD is not due to 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 performance could also be improved. Table 6 reports an ablation study for DAG-GNN and NOTEARS. First, the SHD improvement is not as important as for GraN-DAG and is almost not statistically significant. The improved SHD does not come close to performance of GraN-DAG. Second, PNS and pruning do not have a significant effect of SID, as was the case for GraN-DAG. The lack of improvement for those methods is probably due to the fact that they are not overfitting like GraNDAG, as the training and validation (unregularized) scores shows. NOTEARS captures only linear relationships, thus it will have a hard time overfitting nonlinear data and DAG-GNN uses a strong form of parameter sharing between its conditional densities which possibly cause underfitting in a setup where all the parameters of the conditionals are sampled independently.
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 a 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, j) for which (Aφ)ij < 0.33 for GraN-DAG and obtained an SHD of 29.4±15.9 and an SID of 85.6±45.7, showing a significant improvement over NOTEARS and DAG-GNN, even without PNS and pruning.

heading: A.4 LARGE SAMPLE SIZE EXPERIMENT
text: 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 fj 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 it 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 is substantial. The present paper focuses on sample sizes typically encountered in the structure/causal learning litterature and leave this question for future work.

heading: 500 CAM 123.5 ± 13.9 1181.2 ± 160.8

heading: A.5 DETAILS ON DATA SETS GENERATION
text: Synthetic data sets: For each data set type, 10 data sets are sampled with 1000 examples each. As the synthetic data introduced in 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 U [−1, 1].
• Gauss-ANM is generated following Xj := fj(XπGj ) + Nj ∀j with mutually independent noises Nj ∼ N (0, σ2j ) ∀j where the functions fj are independently sampled from a Gaussian process with a unit bandwidth RBF kernel and σ2j ∼ U [0.4, 0.8]. Source nodes are Gaussian with zero mean and variance sampled from U [1, 2]
• LIN is generated following Xj |XπGj ∼ w T j XπGj + 0.2 · N (0, σ2j ) ∀j where σ2j ∼ U [1, 2] and wj is a vector of |πGj | coefficients each sampled from U [0, 1].
• ADD-FUNC is generated following Xj |XπGj ∼ ∑ i∈πGj
fj,i(Xi) + 0.2 · N (0, σ2j ) ∀j where σ2j ∼ U [1, 2] and the functions fj,i are independently sampled from a Gaussian process with bandwidth one. This model is adapted from Bühlmann et al. (2014).
• PNL-GP is generated following Xj |XπGj ∼ σ(fj(XπGj ) + Laplace(0, lj)) ∀j with the functions fj independently sampled from a Gaussian process with bandwidth one and lj ∼ U [0, 1]. In the two-variable case, this model is identifiable following the Corollary 9 from Zhang & Hyvärinen (2009). To get identifiability according to this corollary, it is important to use non-Gaussian noise, explaining our design choices.
• PNL-MULT is generated following Xj |XπGj ∼ exp(log( ∑ i∈πGj
Xi) + |N (0, σ2j )|) ∀j where σ2j ∼ U [0, 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: //bioinformatics.intec.ugent.be/kmarchal/SynTReN/index.html) using the
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 a 20 nodes graph.
Graph types: Erdős-Rényi (ER) graphs are generated by randomly sampling a topological order and by adding directed edges were it is allowed independently with probability p = 2ed2−d were e is the expected number of edges in the resulting DAG. Scale-free (SF) graphs were generated using the Barabási-Albert model (Barabási & Albert, 1999) which is based on preferential attachment. Nodes are added one by one. Between the new node and the existing nodes, m edges (where m is equal to d or 4d) will be added. An existing node i have the probability p(ki) = ki∑
j kj to be chosen,
where ki represents the degree of the node i. While ER graphs have a degree distribution following a Poisson distribution, SF graphs have a degree distribution following a power law: few nodes, often called hubs, have a high degree. Barabási (2009) have stated that these types of graphs have similar properties to real-world networks which can be found in many different fields, although these claims remain controversial (Clauset et al., 2009).

heading: A.6 SUPPLEMENTARY EXPERIMENTS
text: Gauss-ANM: The results for 20 and 100 nodes are presented in Table 8 and 9 using the same GaussANM 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 Aφ as training proceeds in a typical run for 10 nodes.
LIN & ADD-FUNC: Experiments with LIN and ADD-FUNC data is reported in Table 12 & 13. The details of their generation are given in Appendix A.5.
PNL-GP & PNL-MULT: Table 14 contains the performance of GraN-DAG and other baselines on post nonlinear data discussed in Section 4.1.
4Note that GSF results are missing for two data set types in Tables 9 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.

heading: A.7 METRICS
text: SHD takes two partially directed acyclic graphs (PDAG) and counts the number of edge for which the edge type differs in both PDAGs. There are four edge types: i ← j, i → j, i −− j and i j. Since this distance is defined over the space of PDAGs, we can use it to compare DAGs with DAGs, DAGs with CPDAGs and CPDAGs with CPDAGs. When comparing a DAG with a CPDAG, having i← j instead of i −− j counts as a 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.
Introduced by Peters & Bühlmann (2015), SID counts the number of interventional distribution of the form p(xi| do(xj = x̂j)) that would be miscalculated using the parent adjustment formula (Pearl, 2009) if we were to use 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 outputting a CPDAG such as GES and PC. Peters & Bü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 & Bühlmann (2015) for more details.

heading: A.8 HYPERPARAMETERS
text: All GraN-DAG runs up to this point were performed using the following set of hyperparameters. We used RMSprop as optimizer with learning rate of 10−2 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 1 hidden layer). Leaky-ReLU is used as activation functions. The NN are initialized using the initialization scheme proposed in Glorot & Bengio (2010) also known as Xavier initialization. We used minibatches of 64 samples. This hyperparameter combination have been selected via a small scale experiment in which many hyperparameter combinations have been tried manually on a single data set of type ER1 with 10 nodes until one yielding a 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 a 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. It (rarely) happens that NOTEARS and DAG-GNN returns a cyclic graph. In 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 pcalg R package. For CAM, we used the the default hyperparameters found in the CAM R package, with default PNS and DAG pruning.

heading: A.9 HYPERPARAMETER SELECTION VIA HELD-OUT SCORE
text: Most structure/causal learning algorithms have hyperparameters which must be selected prior to learning. For instance, NOTEARS and GES have a 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.e. 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 a DAG which is too dense or too sparse since if the estimated graph contains many spurious edges, the score on the held-out data set 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 combinations 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 early stopping and hyperparameter selection). Note that the held-out score is always evaluated without the regularizing term (e.g. the held-out score of NOTEARS was evaluated without its L1 regularizer).
The symbols ++ and + indicate the hyperparameter search improved performance against default hyperparameter runs above one standard deviation and within one 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
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 hyperparameter selection procedure (although rarely significantly). This might be explained by the fact that neural networks are in general very sensitive to the choice of hyperparameters. However, not all methods improved their performance and no method improves its performance in all settings. GES and GSF for instance, often have significantly worse results. This might be due to 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.
Additional implementation details for held-out score evaluation: GraN-DAG makes use of a final pruning step to remove spurious edges. One could simply mask the inputs of the NN corresponding to removed edges and evaluate the held-out score. However, doing so yields an unrepresentative score since some masked inputs have an important role in the learned function and once these inputs are masked, the quality of the fit 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, we evaluate the held-out score with this retrained architecture. During this retraining phase, the estimated graph is fixed, only the conditional densities are relearned. Since NOTEARS and DAG-GNN are not always guaranteed to return a DAG (although they almost always do), some extra thresholding 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 is also retrained after its pruning phase prior to evaluating its held-out score.

Review summaries:
The 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. 

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.