Title: GRADIENT-BASED NEURAL DAG LEARNING

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.

1. INTRODUCTION:

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 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 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 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 (Bü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.

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 be described by a random vector X ∈ R d entailed by an underlying CGM (P X , G) where P X 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 π G j denote the set of parents of node j in G and let X π G j 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 P X is said to be Markov to G, i.e. we can write the probability density function (pdf) of P X as p(x) = d j=1 p j (x j |x π G j ) where p j (x j |x π G j ) is the conditional pdf of variable X j given X π G j . 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) .

2.2. STRUCTURE IDENTIFIABILITY:

In general, it is impossible to recover G given only samples from P X , 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 = (P X , G), its graph G is said to be identifiable from P X if there exists no other CGM M = ( PX , G) satisfying all assumptions in A such that G = G and PX = P X .

There are many examples of graph identifiability results for continuous variables (Peters et al., 2014; Peters & Bühlman, 2014; Shimizu et al., 2006; Zhang & Hyvä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 P X is entailed by a structural equation model of the form

X j := f j (X π G j ) + N j with N j ∼ N (0, σ 2 j ) ∀j ∈ V (1)

where f j is a nonlinear function satisfying some mild regularity conditions and the noises N j are mutually independent, then G is identifiable from P X (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.

2.3. NOTEARS: CONTINUOUS OPTIMIZATION FOR STRUCTURE LEARNING:

Structure learning is the problem of learning G using a data set of n samples {x (1) , ..., x (n) } from P X . Score-based approaches cast this problem as an optimization problem, i.e. Ĝ = arg max G∈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 = [u 1 | . . . |u d ] ∈ R d×d which represents (possibly negative) coefficients in a linear SEM of the form X j := u j X + N i ∀j where N j is a noise variable. Let G U be the directed graph associated with the SEM and let A U be the (binary) adjacency matrix associated with G U . One can see that the following equivalence holds:

EQUATION

To make sure G U is acyclic, the authors propose the following constraint on U :

Tr e U U -d = 0 (3) where e M ∞ k=0 M k

k! is the matrix exponential and is the Hadamard product. To see why this constraint characterizes acyclicity, first note that (A U k ) jj is the number of cycles of length k passing through node j in graph G U . Clearly, for G U to be acyclic, we must have Tr A U 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 2 F -λ U 1 s.t. Tr e U U -d = 0 (4)

where X ∈ R n×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).

3. GRAN-DAG: GRADIENT-BASED NEURAL DAG LEARNING:

We propose a new nonlinear extension to the framework presented in Section 2.3. For each variable X j , 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 ∈ R d , i.e. the vector X with the jth component masked to zero, and outputs θ (j) ∈ R m , the m-dimensional parameter vector of the desired distribution family for variable X j . 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 p j (x j |x -j ; φ (j) ). Note that the product d j=1 p j (x j |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 matrix A φ similar to the one encountered in Section 2.3, which can be directly used inside the constraint of Equation 3 to enforce acyclicity.

3.1. NEURAL NETWORK CONNECTIVITY:

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) | ∈ R m×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 C ki is the sum of all NN path products from input i to output k. This means it is sufficient to have C ki = 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 X i , it is sufficient to have m k=1 C ki = 0.

3.2. A WEIGHTED ADJACENCY MATRIX:

We now define a weighted adjacency matrix A φ that can be used in constraint of Equation 3.

EQUATION

where C (j) denotes the connectivity matrix of the NN associated with variable X j .

As the notation suggests, A φ ∈ R d×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 X i (8) We note G φ to be the directed graph entailed by parameter φ. We can now write our adapted acyclicity constraint:

EQUATION

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).

3.3. A DIFFERENTIABLE SCORE AND ITS OPTIMIZATION:

We propose solving the maximum likelihood optimization problem

EQUATION

where π φ j denotes the set of parents of node j in graph G φ . Note that

d j=1 log p j (X j |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

EQUATION

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 (P X , 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 (P X , G) ∈ C and C = M(A) for a specific set of assumptions A such that G is identifiable from P X , then G φ * = G.

Proof: Let P φ be the joint distribution entailed by parameter φ. Note that the population loglikelihood E X∼P X log p φ (X) is maximal iff P φ = P X . We know this maximum can be achieved by a specific parameter φ * since by hypothesis (P X , G) ∈ C. Since G is identifiable from P X , we know there exists no other CGM ( PX , G) ∈ C such that G = G and PX = P X . 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).

3.4. THRESHOLDING:

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 = 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)

EQUATION

(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).

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 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 (Bü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 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 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.

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 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(d 3 ). NOTEARS justifies the use of a batch proximal quasi-Newton algorithm by the low number of O(d 3 ) 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.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(d 2 ) 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 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 (Bü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

X j |X π G j ∼ N (μ (j) , σ2 (j) ) ∀j.

The parameters σ2 (j) are learned as well, but do not depend on the parent variables X π G j (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 Erdős-Ré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 & Bü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(x j |do(X i = 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.

4.1. SYNTHETIC DATA:

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: Erdős-Ré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

X j := f j (X π G j ) + N j with mutually independent noises N j ∼ N (0, σ 2 j

) ∀j where the functions f j are independently sampled from a Gaussian process with a unit bandwidth RBF kernel and with σ 2 j 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 f j 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 f j 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 f j are additive, i.e.

f j (x π G j ) = i∈π G j f ij (x j

) ∀j, it manages to compete with GraN-DAG which does not make this incorrect modelling assumption 2 . 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∈π G j f ij (X i

) 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 j := g j (f j (X π G j ) + N j ) where N j 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 g j , 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.

4.2. REAL AND PSEUDO-REAL 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 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.

In applications, it is not clear whether the conditions of Proposition 1 hold since we do not know whether (P X , 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 π G i 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.

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.

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 & 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-based score. However, RESIT does not scale well to graph of more than 20 nodes. CAM (Bü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. Bü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.

6. CONCLUSION:

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.

A APPENDIX:

A.1 OPTIMIZATION Let us recall the augmented Lagrangian:

EQUATION

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 ) ≈ ∇ φ LB (φ, λ t , µ t ) with LB (φ, λ t , µ t ) 1 |B| x∈B d i=1 log p i (x i |x π φ i ; φ (i) ) -λ t h(φ) - µ t 2 h(φ) 2 (13)

where B is a minibatch sampled from the data set and |B| is the minibatch size. The gradient estimate ∇ φ LB (φ, λ t , µ t ) can be computed using standard deep learning libraries. We consider a subproblem has converged when LH (φ, λ 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:

EQUATION

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 K t be the number of iterations executed to approximately solve it. Let I = T t=1 K t 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. 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 until A φ * 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 X i . 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 X i and not a necessary one.

To avoid this issue, we consider the following alternative. Consider the function L : R d → R d which maps all d variables to their respective conditional likelihoods, i.e. L i (X)

p i (X i | X π φ * T i ) ∀i.

We consider the following expected Jacobian matrix

J E X∼P X ∂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 J ij can be loosely interpreted as the strength of edge (i, j). We propose removing edges starting from the lowest J ij 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.

A.3 PRELIMINARY NEIGHBORHOOD SELECTION AND DAG PRUNING:

PNS: For graphs of 50 nodes or more, GraN-DAG performs a preliminary neighborhood selection (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 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, GraN-DAG performs a pruning step identical to CAM (Bühlmann et 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 i , a generalized additive model is fitted against the current parents of X i 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 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 GraN-DAG, 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. 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: 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 = 2e d 2 -d were e is the expected number of edges in the resulting DAG. Scale-free (SF) graphs were generated using the Barabási-Albert model (Barabá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(k i ) = ki j kj to be chosen, where k i 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. Barabá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) .

A.6 SUPPLEMENTARY EXPERIMENTS:



Gauss-ANM::

The results for 20 and 100 nodes are presented in Table 8 and 9 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 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. 14 contains the performance of GraN-DAG and other baselines on post nonlinear data discussed in Section 4.1. Introduced by Peters & Bühlmann (2015) , SID counts the number of interventional distribution of the form p(x i | do(x j = xj )) 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 & 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.

PNL-GP & PNL-MULT: Table:

A.8 HYPERPARAMETERS 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.

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 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. 

:

Not all parameter vectors need to have the same dimensionality, but to simplify the notation, we suppose mj = m ∀j

Although 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.

This was the default value of thresholding used in NOTEARS and DAG-GNN.

Note that GSF results are missing for two data set types in Tables9 and 14. This is because the search algorithm could not finish within 12 hours, even when the maximal in-degree was limited to

All other methods could run in less than 6 hours.

Figures and Tables:

Figure1: Entries of the weighted adjacency matrix A φ as training proceeds in GraN-DAG for a 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 -4 is the threshold introduced in Section 3.4. We can see that GraN-DAG successfully recovers most edges correctly while keeping few spurious edges.

Results for ER and SF graphs of 10 nodes with Gauss-ANM data

Table content: <table><tr><td/><td>ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td></tr><tr><td/><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td colspan="2">GraN-DAG 1.7±2.5</td><td>1.7±3.1</td><td>8.3±2.8</td><td>21.8±8.9</td><td>1.2±1.1</td><td>4.1±6.1</td><td>9.9±4.0</td><td>16.4±6.0</td></tr><tr><td colspan="3">DAG-GNN 11.4±3.1 37.6±14.4</td><td colspan="2">35.1±1.5 81.9±4.7</td><td>9.9±1.1</td><td colspan="3">29.7±15.8 20.8±1.9 48.4±15.6</td></tr><tr><td colspan="3">NOTEARS 12.2±2.9 36.6±13.1</td><td colspan="2">32.6±3.2 79.0±4.1</td><td>10.7±2.2</td><td colspan="3">32.0±15.3 20.8±2.7 49.8±15.6</td></tr><tr><td>CAM</td><td>1.1±1.1</td><td>1.1±2.4</td><td colspan="2">12.2±2.7 30.9±13.2</td><td>1.4±1.6</td><td>5.4±6.1</td><td>9.8±4.3</td><td>19.3±7.5</td></tr><tr><td>GSF</td><td>6.5±2.6</td><td>[6.2±10.8</td><td colspan="3">21.7±8.4 [37.2±19.2 1.8±1.7</td><td>[2.0±5.1</td><td>8.5±4.2</td><td>[13.2±6.8</td></tr><tr><td/><td/><td>17.7±12.3]</td><td/><td>62.7±14.9]</td><td/><td>6.9±6.2]</td><td/><td>20.6±12.1]</td></tr><tr><td colspan="3">RANDOM 26.3±9.8 25.8±10.4</td><td colspan="2">31.8±5.0 76.6±7.0</td><td colspan="4">25.1±10.2 24.5±10.5 28.5±4.0 47.2±12.2</td></tr></table>

Results for ER and SF graphs of 50 nodes with Gauss-ANM data

Table content: <table><tr><td/><td>ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td/><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td colspan="2">GraN-DAG 5.1±2.8</td><td>22.4±17.8</td><td>102.6±21.2</td><td>1060.1±109.4</td><td>25.5±6.2</td><td>90.0±18.9</td><td>111.3±12.3</td><td>271.2±65.4</td></tr><tr><td colspan="2">DAG-GNN 49.2±7.9</td><td colspan="2">304.4±105.1 191.9±15.2</td><td>2146.2±64</td><td>49.8±1.3</td><td>182.8±42.9</td><td>144.9±13.3</td><td>540.8±151.1</td></tr><tr><td colspan="2">NOTEARS 62.8±9.2</td><td colspan="2">327.3±119.9 202.3±14.3</td><td>2149.1±76.3</td><td>57.7±3.5</td><td>195.7±54.9</td><td>153.7±11.8</td><td>558.4±153.5</td></tr><tr><td>CAM</td><td>4.3±1.9</td><td>22.0±17.9</td><td>98.8±20.7</td><td>1197.2±125.9</td><td>24.1±6.2</td><td>85.7±31.9</td><td>111.2±13.3</td><td>320.7±152.6</td></tr><tr><td>GSF</td><td>25.6±5.1</td><td>[21.1±23.1</td><td>81.8±18.8</td><td>[906.6±214.7</td><td>31.6±6.7</td><td>[85.8±29.9</td><td>120.2±10.9</td><td>[284.7±80.2</td></tr><tr><td/><td/><td>79.2±33.5]</td><td/><td>1030.2±172.6]</td><td/><td>147.3±49.9]</td><td/><td>379.9±98.3]</td></tr><tr><td colspan="5">RANDOM 535.7±401.2 272.3±125.5 708.4±234.4 1921.3±203.5</td><td colspan="4">514.0±360.0 381.3±190.3 660.6±194.9 1198.9±304.6</td></tr></table>

Results on real and pseudo-real data

Table content: <table><tr><td/><td colspan="3">Protein signaling data set</td><td colspan="2">SynTReN (20 nodes)</td><td/></tr><tr><td/><td colspan="3">SHD SHD-C SID</td><td>SHD</td><td>SHD-C</td><td>SID</td></tr><tr><td>GraN-DAG</td><td>13</td><td>11</td><td>47</td><td>34.0±8.5</td><td>36.4±8.3</td><td>161.7±53.4</td></tr><tr><td colspan="2">GraN-DAG++ 13</td><td>10</td><td>48</td><td>33.7±3.7</td><td>39.4±4.9</td><td>127.5±52.8</td></tr><tr><td>DAG-GNN</td><td>16</td><td>21</td><td>44</td><td>93.6±9.2</td><td>97.6±10.3</td><td>157.5±74.6</td></tr><tr><td>NOTEARS</td><td>21</td><td>21</td><td>44</td><td colspan="3">151.8±28.2 156.1±28.7 110.7±66.7</td></tr><tr><td>CAM</td><td>12</td><td>9</td><td>55</td><td>40.5±6.8</td><td>41.4±7.1</td><td>152.3±48</td></tr><tr><td>GSF</td><td>18</td><td>10</td><td colspan="2">[44, 61] 61.8±9.6</td><td>63.3±11.4</td><td>[76.7±51.1, 109.9±39.9]</td></tr><tr><td>GES</td><td>26</td><td>28</td><td colspan="2">[34, 45] 82.6±9.3</td><td>85.6±10</td><td>[157.2±48.3, 168.8±47.8]</td></tr><tr><td>PC</td><td>17</td><td>11</td><td colspan="2">[47, 62] 41.2±5.1</td><td>42.4±4.6</td><td>[154.8±47.6, 179.3±55.6]</td></tr><tr><td>RANDOM</td><td>21</td><td>20</td><td>60</td><td>84.7±53.8</td><td>86.7±55.8</td><td>175.8±64.7</td></tr></table>

Total number of iterations (×10 3 ) before augmented Lagrangian converges on Gauss-ANM

Table content: <table><tr><td>data.</td><td/><td/><td/></tr><tr><td colspan="4">20 nodes ER1 20 nodes ER4 100 nodes ER1 100 nodes ER4</td></tr><tr><td>GraN-DAG 27.3 ± 3.6</td><td>30.4 ± 4.2</td><td>23.1 ± 0.7</td><td>23.1 ± 0.8</td></tr><tr><td>NOTEARS 67.1 ± 35.3</td><td>72.3 ± 24.3</td><td>243.6 ± 12.3</td><td>232.4 ± 12.9</td></tr><tr><td colspan="2">A.2 THRESHOLDING TO ENSURE ACYCLICITY</td><td/><td/></tr></table>

PNS and pruning ablation study for GraN-DAG (averaged over 10 datasets from ER1 with

Table content: <table><tr><td>50 nodes)</td><td/><td/><td/></tr><tr><td colspan="2">PNS Pruning SHD</td><td>SID</td><td colspan="2">NLL (train) NLL (validation)</td></tr><tr><td>False False</td><td colspan="3">1086.8±48.8 31.6±23.6 0.36±0.07</td><td>1.44±0.21</td></tr><tr><td>True False</td><td>540.4±70.3</td><td colspan="2">17.4±16.7 0.52±0.08</td><td>1.16±0.17</td></tr><tr><td>False True</td><td>11.8±5.0</td><td colspan="2">39.7±25.5 0.78±0.12</td><td>0.84±0.12</td></tr><tr><td>True True</td><td>6.1±3.3</td><td colspan="2">29.3±19.5 0.78±0.13</td><td>0.83±0.12</td></tr></table>



Table content: <table><tr><td colspan="5">: PNS and pruning ablation study for DAG-GNN and NOTEARS (averaged over 10 datasets</td></tr><tr><td colspan="2">from ER1 with 50 nodes)</td><td/><td/></tr><tr><td>Algorithm</td><td colspan="2">PNS Pruning SHD</td><td>SID</td><td>Score (train) Score (validation)</td></tr><tr><td colspan="2">DAG-GNN False False</td><td colspan="3">56.8±11.1 322.9±103.8 -2.8±1.5</td><td>-2.2±1.6</td></tr><tr><td/><td>True False</td><td colspan="3">55.5±10.2 314.5±107.6 -2.1±1.6</td><td>-2.1±1.7</td></tr><tr><td/><td>False True</td><td>49.4±7.8</td><td colspan="2">325.1±103.7 -1.8±1.1</td><td>-1.8±1.2</td></tr><tr><td/><td>True True</td><td>47.7±7.3</td><td colspan="2">316.5±105.6 -1.9±1.6</td><td>-1.9±1.6</td></tr><tr><td colspan="2">NOTEARS False False</td><td>64.2±9.5</td><td colspan="2">327.1±110.9 -23.1±1.8</td><td>-23.2±2.1</td></tr><tr><td/><td>True False</td><td colspan="3">54.1±10.9 321.5±104.5 -25.2±2.7</td><td>-25.4±2.8</td></tr><tr><td/><td>False True</td><td>49.5±8.8</td><td colspan="2">327.7±111.3 -26.7±2.0</td><td>-26.8±2.1</td></tr><tr><td/><td>True True</td><td>49.0±7.6</td><td colspan="2">326.4±106.9 -26.23±2.2</td><td>-26.4±2.4</td></tr></table>

Results for ER and SF graphs of 20 nodes with Gauss-ANM data

Table content: <table><tr><td>ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr></table>

SHD for GES and PC (against GraN-DAG for reference) with Gauss-ANM data

Table content: <table><tr><td/><td>10 nodes</td><td/><td>20 nodes</td><td/><td>50 nodes</td><td/><td>100 nodes</td></tr><tr><td/><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td></tr><tr><td colspan="2">GraN-DAG 1.7±2.5</td><td>8.3±2.8</td><td>4.0 ±3.4</td><td colspan="2">45.2±10.7 5.1±2.8</td><td colspan="2">102.6±21.2 15.1±6.0</td><td>206.6±31.5</td></tr><tr><td>GES</td><td colspan="3">13.8±4.8 32.3±4.3 43.3±12.4</td><td>94.6±9.8</td><td colspan="4">106.6±24.7 254.4±39.3 292.9±33.6 542.6±51.2</td></tr><tr><td>PC</td><td>8.4±3</td><td>34±2.6</td><td colspan="2">20.136.4±6.5 73.1±5.8</td><td>44.0±11.6</td><td>183.6±20</td><td>95.2±9.1</td><td>358.8±20.6</td></tr><tr><td/><td>SF1</td><td>SF4</td><td>SF1</td><td>SF4</td><td>SF1</td><td>SF4</td><td>SF1</td><td>SF4</td></tr><tr><td colspan="2">GraN-DAG 1.2±1.1</td><td>9.9±4.0</td><td>7.6±2.5</td><td>36.8±5.1</td><td>25.5±6.2</td><td colspan="2">111.3±12.3 59.2±7.7</td><td>262.7±19.6</td></tr><tr><td>GES</td><td>8.1±2.4</td><td colspan="2">17.4±4.5 26.2±7.5</td><td>50.7±6.2</td><td>73.9±7.4</td><td colspan="2">178.8±16.5 190.3±22</td><td>408.7±24.9</td></tr><tr><td>PC</td><td>4.8±2.4</td><td colspan="2">16.4±2.8 13.6±2.1</td><td>44.4±4.6</td><td>43.1±5.7</td><td colspan="2">135.4±10.7 97.6±6.6</td><td>314.2±17.5</td></tr></table>

Lower and upper bound on the SID for GES and PC (against GraN-DAG for reference) with Gauss-ANM data. See Appendix A.7 for details on how to compute SID for CPDAGs.

Table content: <table><tr><td>10 nodes</td><td/><td>20 nodes</td><td/><td>50 nodes</td><td/><td>100 nodes</td><td/></tr><tr><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td><td>ER1</td><td>ER4</td></tr></table>



Table content: <table><tr><td>: LIN</td></tr></table>



Table content: <table><tr><td>: ADD-FUNC</td></tr></table>

Gauss-ANM -10 nodes with hyperparameter search

Table content: <table><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">GraN-DAG 1.0 ± 1.6 +</td><td>0.4 ± 1.3 ++</td><td>5.5 ± 2.8 +</td><td>9.7 ± 8.0 ++</td><td>1.3 ± 1.8 -</td><td>3.0 ± 3.4 +</td><td>9.6 ± 4.5 +</td><td>15.1 ± 6.1 +</td></tr><tr><td colspan="2">DAG-GNN 10.9 ± 2.6 +</td><td>35.5 ± 13.6 +</td><td colspan="2">38.3 ± 2.9 --84.4 ± 3.5 -</td><td>9.9 ± 1.7 +</td><td colspan="2">30.3 ± 18.8 -21.4 ± 2.1 -</td><td>44.0 ± 15.5 +</td></tr><tr><td colspan="3">NOTEARS 26.7 ± 6.9 --35.2 ± 10.6 +</td><td colspan="5">20.9 ± 6.6 ++ 62.0 ± 6.7 ++ 20.4 ± 9.6 --38.8 ± 16.7 -26.9 ± 7.4 -</td><td>61.1 ± 13.8 -</td></tr><tr><td>CAM</td><td>3.0 ± 4.2 -</td><td>2.2 ± 5.7 -</td><td colspan="2">7.7 ± 3.1 ++ 23.2 ± 14.7 +</td><td>2.4 ± 2.5 -</td><td>5.2 ± 5.5 +</td><td>9.6 ± 3.1 +</td><td>20.1 ± 6.8 -</td></tr><tr><td>GSF</td><td>5.3 ± 3.3 +</td><td>[8.3 ± 13.2 +</td><td>23.1 ± 7.9 -</td><td>[56.1 ± 20.4 -</td><td>3.3 ± 2.5 -</td><td>[7.0 ± 11.6 -</td><td colspan="2">14.2 ± 5.6 --[26.2 ± 11.1 -</td></tr><tr><td/><td/><td>15.4 ± 13.5]</td><td/><td>65.1 ± 19.3]</td><td/><td>12.2 ± 11.0]</td><td/><td>36.9 ± 21.6]</td></tr><tr><td>GES</td><td colspan="2">38.6 ± 2.1 --[20.3 ± 15.4 +</td><td>33.0 ± 3.4 -</td><td>[66.2 ± 7.0 +</td><td colspan="2">38.3 ± 2.4 --[8.8 ± 5.2 -</td><td colspan="2">33.6 ± 4.8 --[32.7 ± 12.7 -</td></tr><tr><td/><td/><td>28.3 ± 18.4]</td><td/><td>76.6 ± 4.3]</td><td/><td>25.5 ± 18.2]</td><td/><td>52.0 ± 14.0]</td></tr></table>

Gauss-ANM -50 nodes with hyperparameter search

Table content: <table><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">GraN-DAG 3.8 ± 3.3 +</td><td>15.0 ± 14.0 +</td><td>105.6 ± 16.5 -</td><td>1131.7 ± 91.0 -</td><td>24.7 ± 6.4 +</td><td>86.5 ± 34.6 +</td><td>112.7 ± 15.5 -</td><td>268.3 ± 85.8 +</td></tr><tr><td colspan="2">DAG-GNN 47.0 ± 7.8 +</td><td colspan="2">268.1 ± 118.0 + 196.2 ± 14.4 -</td><td>1972.8 ± 110.6 ++</td><td>51.8 ± 5.6 -</td><td>166.5 ± 48.9 +</td><td>144.2 ± 11.6 +</td><td>473.4 ± 105.4 +</td></tr><tr><td colspan="3">NOTEARS 193.5 ± 77.3 --326.0 ± 99.1 +</td><td>369.5 ± 81.9 --</td><td>2062.0 ± 107.7 +</td><td colspan="2">104.8 ± 22.4 --290.3 ± 136.8 -</td><td>213.0 ± 35.1 --</td><td>722.7 ± 177.3 -</td></tr><tr><td>CAM</td><td>4.0 ± 2.7 +</td><td>21.1 ± 22.1 +</td><td>105.6 ± 20.9 -</td><td>1225.9 ± 205.7 -</td><td>23.8 ± 6.0 +</td><td>81.5 ± 15.3 +</td><td>112.2 ± 14.0 -</td><td>333.8 ± 156.0 -</td></tr><tr><td>GSF</td><td>24.9 ± 7.4 +  *</td><td>[40.0 ± 26.3 - *</td><td>129.3 ± 20.4 -- *</td><td>[1280.8 ± 202.3 -- *</td><td>35.3 ± 6.9 - *</td><td>[99.7 ± 41.7 - *</td><td>121.6 ± 11.7 - *  *  *</td><td>[310.8 ± 108.1 - *  *  *</td></tr><tr><td/><td/><td>77.5 ± 45.3]</td><td/><td>1364.1 ± 186.7]</td><td/><td>151.9 ± 59.7]</td><td/><td>391.9 ± 93.3]</td></tr><tr><td>GES</td><td colspan="2">1150.1 ± 9.8 --[112.7 ± 71.1 +</td><td colspan="2">1066.1 ± 11.7 --[1394.3 ± 81.8 ++</td><td colspan="2">1161.7 ± 7.0 --[322.8 ± 211.1 -</td><td colspan="2">1116.1 ± 14.2 --[1002.7 ± 310.9 --</td></tr><tr><td/><td/><td>132.0 ± 89.0]</td><td/><td>1464.8 ± 63.8]</td><td/><td>336.0 ± 215.4]</td><td/><td>1094.0 ± 345.1]</td></tr></table>

Gauss-ANM -20 nodes with hyperparameter search NOTEARS 101.5 ± 39.6 --100.4 ± 47.0 + 124.0 ± 16.3 --267.0 ± 46.5 ++ 55.0 ± 28.2 --87.6 ± 26.9 - 66.7 ± 8.3 -- 154.6 ± 43.0 - CAM 2.8 ± 2.2 - 11.5 ± 10.2 -64.3 ± 29.3 - 121.7 ± 73.1 + 5.5 ± 1.6 + 19.3 ± 7.8 + 36.0 ± 5.1 - 66.3 ± 28.6 - GSF 11.6 ± 3.0 + [26.4 ± 13.3 - 49.8 ± 26.5]

Table content: <table><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">GraN-DAG 2.7 ± 2.3 +</td><td>9.6 ± 10.3 +</td><td>35.9 ± 11.8 +</td><td colspan="2">120.4 ± 37.0 ++ 6.5 ± 2.4 +</td><td>17.5 ± 6.3 ++</td><td>35.6 ± 4.1 +</td><td>54.8 ± 14.3 +</td></tr><tr><td colspan="2">DAG-GNN 21.0 ± 6.1 +</td><td>98.8 ± 42.2 +</td><td>77.2 ± 6.5 -</td><td>345.6 ± 18.6 -</td><td>19.1 ± 0.7 +</td><td>55.0 ± 20.1 +</td><td>50.2 ± 5.4 -</td><td>118.7 ± 33.2 -</td></tr><tr><td/><td/><td/><td>46.2 ± 12.6 -</td><td>[172.7 ± 40.8 -</td><td>12.8 ± 2.1 --</td><td>[32.1 ± 14.0 --</td><td>42.3 ± 5.1 -</td><td>[68.9 ± 27.7 -</td></tr><tr><td/><td/><td/><td/><td>213.5 ± 38.6]</td><td/><td>56.2 ± 13.8]</td><td/><td>95.1 ± 33.8]</td></tr><tr><td>GES</td><td>169.9 ± 5.0 --</td><td>[45.4 ± 29.2 +</td><td>142.8 ± 7.7 --</td><td/><td/><td/><td/><td/></tr><tr><td/><td/><td>57.2 ± 36.6]</td><td/><td/><td/><td/><td/><td/></tr></table>

Gauss-ANM -100 nodes with hyperparameter search DAG 15.1 ± 7.5 + 65.1 ± 33.2 + 191.6 ± 17.8 + 4090.7 ± 418.0 + 51.6 ± 10.2 + 210.6 ± 51.9 ++ 255.7 ± 21.1 + 790.5 ± 159.7 + DAG-GNN 103.9 ± 9.1 + 757.6 ± 215.0 + 387.1 ± 25.3 - 7741.9 ± 522.5 + 103.5 ± 8.2 - 391.7 ± 60.0 + 314.8 ± 16.3 + 1257.3 ± 185.2 + NOTEARS 421.3 ± 207.0 --945.7 ± 339.7 - 631.1 ± 136.6 --8272.4 ± 444.2 - 244.3 ± 63.8 -- 815.6 ± 346.5 - 482.3 ± 114.1 --1929.7 ± 363.1 --

Table content: <table><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>SF1</td><td/><td>SF4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>GraN-CAM</td><td>12.3 ± 4.9 ++</td><td>128.0 ± 66.3 -</td><td>198.8 ± 22.2 -</td><td>4602.2 ± 523.7 -</td><td>51.1 ± 9.4 +</td><td>233.6 ± 62.3 -</td><td>255.7 ± 22.2 -</td><td>851.4 ± 206.0 -</td></tr><tr><td>GSF</td><td>100.2 ± 9.9 -- *  *</td><td>[719.8 ± 242.1 -- *  *</td><td>387.6 ± 23.9 *  *  *</td><td>[7535.1 ± 595.2 *  *  *</td><td>67.3 ± 14.0 +  *  *  *</td><td>[254.5 ± 35.4 - *  *  *</td><td>315.1 ± 16.7 -- *  *  *</td><td>[1214.0 ± 156.4 -- *  *  *</td></tr><tr><td/><td/><td>721.1 ± 242.9]</td><td/><td>7535.1 ± 595.2]</td><td/><td>340.4 ± 70.4]</td><td/><td>1214.0 ± 156.4]</td></tr><tr><td>GES</td><td colspan="2">4782.5 ± 22.9 --[362.3 ± 267.7 +</td><td colspan="2">4570.1 ± 27.9 --[5400.7 ± 299.2 ++</td><td colspan="2">4769.1 ± 26.7 --[1311.1 ± 616.6 --</td><td colspan="2">4691.3 ± 47.3 --[3882.7 ± 1010.6 --</td></tr><tr><td/><td/><td>384.1 ± 293.6]</td><td/><td>5511.5 ± 308.5]</td><td/><td>1386.2 ± 713.9]</td><td/><td>3996.7 ± 1075.7]</td></tr></table>

PNL-GP with hyperparameter search DAG 1.2 ± 2.2 + 1.9 ± 4.2 + 9.8 ± 4.9 + 29.0 ± 17.6 + 12.8 ± 4.9 + 55.3 ± 24.2 + 73.9 ± 16.8 - 1107.2 ± 144.7 + DAG-GNN 10.6 ± 4.9 + 35.8 ± 19.6 - 38.6 ± 2.0 --82.2 ± 5.7 -- 48.1 ± 8.4 + 330.4 ± 69.9 + 192.5 ± 19.2 + 2079.5 ± 120.9 + NOTEARS 20.6 ± 11.4 -30.5 ± 18.8 + 24.2 ± 6.5 ++ 66.4 ± 6.9 ++ 102.1 ± 27.3 --299.8 ± 85.8 + 660.0 ± 258.2 --1744.0 ± 232.9 ++ CAM 2.7 ± 4.0 - 6.4 ± 11.8 + 8.7 ± 4.5 - 30.9 ± 20.4 -4.0 ± 2.4 + 10.7 ± 12.4 + 52.3 ± 8.5 - 913.9 ± 209.3 - GSF 12.9 ± 3.9 -- * * *

Table content: <table><tr><td>#Nodes</td><td>10</td><td/><td/><td/><td>50</td><td/><td/><td/></tr><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>ER1</td><td/><td>ER4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="3">GraN-[10.5 ± 8.7 -- *  *  *</td><td>40.7 ± 1.3 -- *  *</td><td>[79.2 ± 3.8 -- *  *</td><td>48.8 ± 3.9 -- *  *</td><td>[281.6 ± 70.7 -- *  *</td><td>199.9 ± 15.2 -- *  *  *</td><td>[1878.0 ± 122.4 -- *  *  *</td></tr><tr><td/><td/><td>53.6 ± 23.8]</td><td/><td>79.2 ± 3.8]</td><td/><td>281.6 ± 70.7]</td><td/><td>1948.4 ± 139.6]</td></tr></table>

PNL-MULT with hyperparameter search DAG-GNN 14.6 ± 3.1 ++ 36.9 ± 10.6 + 38.9 ± 2.0 - 85.8 ± 1.2 -- 64.3 ± 27.8 + 508.8 ± 317.2 - 212.5 ± 12.3 ++ 2216.9 ± 95.6 + NOTEARS 28.8 ± 9.1 --30.3 ± 11.8 + 35.4 ± 3.8 + 78.4 ± 7.5 + 160.2 ± 67.5 --443.5 ± 205.1 -

Table content: <table><tr><td>#Nodes</td><td>10</td><td/><td/><td/><td>50</td><td/><td/><td/></tr><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>ER1</td><td/><td>ER4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">GraN-DAG 10.0 ± 4.5 +</td><td>29.1 ± 9.7 +</td><td>32.9 ± 3.3 -</td><td>76.7 ± 4.1 +</td><td>59.8 ± 28.2 -</td><td>213.6 ± 97.3 +</td><td>272.1 ± 69.4 -</td><td>2021.6 ± 185.8 +</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td>229.2 ± 25.4 -</td><td>2158.8 ± 70.3 +</td></tr><tr><td>CAM</td><td>17.2 ± 8.0 -</td><td>33.7 ± 14.4 +</td><td>32.3 ± 6.5 +</td><td>76.6 ± 8.2 +</td><td>97.5 ± 71.1 -</td><td colspan="2">282.3 ± 123.8 + 251.0 ± 25.9 --</td><td>2026.2 ± 58.2 +</td></tr><tr><td>GSF</td><td>15.6 ± 4.4 -- *  *</td><td>[10.0 ± 6.3 -- *  *</td><td>39.3 ± 2.2 -- *  *</td><td>[76.0 ± 9.6 -- *  *</td><td>66.4 ± 14.4 -- *  *  *</td><td>[145.1 ± 96.1 -- *  *  *</td><td>&gt; 12 hours</td><td/></tr><tr><td/><td/><td>60.1 ± 17.2]</td><td/><td>79.9 ± 5.3]</td><td/><td>618.8 ± 257.0]</td><td/><td/></tr></table>

LIN with hyperparameter search DAG 10.1 ± 3.9 - 28.7 ± 14.7 -34.7 ± 2.9 -- 79.5 ± 4.4 - 40.8 ± 10.3 -236.3 ± 101.7 + 256.9 ± 55.7 -- 2151.4 ± 144.3 - DAG-GNN 9.0 ± 2.7 ++ 35.6 ± 11.4 -19.6 ± 4.6 + 63.9 ± 7.5 - 48.3 ± 6.8 + 381.7 ± 145.4 - 149.7 ± 17.2 ++ 2070.7 ± 51.9 -- NOTEARS 14.0 ± 4.1 -- * 32.2 ± 7.9 + * 20.7 ± 5.1 ++ 63.1 ± 8.0 ++ 87.7 ± 44.3 - * 294.3 ± 99.3 + * 200.3 ± 67.1 -- 1772.7 ± 143.7 ++ CAM 8.8 ± 6.0 + 25.8 ± 13.5 + 33.9 ± 2.8 - 77.1 ± 4.5 + 34.8 ± 7.0 + 221.2 ± 98.3 + 202.2 ± 14.3 -- 1990.8 ± 97.5 -

Table content: <table><tr><td>#Nodes</td><td>10</td><td/><td/><td/><td>50</td><td/><td/><td/></tr><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>ER1</td><td/><td>ER4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>GraN-GSF</td><td>10.7 ± 3.5 -</td><td>[15.8 ± 8.4 -</td><td>33.4 ± 3.3 ++</td><td>[71.7 ± 11.5 +</td><td>54.4 ± 6.5 -- *</td><td>[158.1 ± 115.9 - *</td><td>195.6 ± 9.9 *  *</td><td>[2004.9 ± 85.2 *  *</td></tr><tr><td/><td/><td>45.2 ± 20.2]</td><td/><td>77.3 ± 6.1]</td><td/><td>560.9 ± 220.7]</td><td/><td>2004.9 ± 85.2]</td></tr></table>

ADD-FUNC with hyperparameter search DAG 2.6 ± 2.4 + 4.3 ± 4.3 + 7.0 ± 3.1 ++ 37.1 ± 12.4 ++ 13.2 ± 6.7 + 72.1 ± 55.2 + 90.1 ± 25.6 - 1241.7 ± 289.8 + DAG-GNN 8.7 ± 2.8 ++ 22.3 ± 9.4 + 25.3 ± 3.8 ++ 63.6 ± 8.6 ++ 44.7 ± 9.7 ++ 306.9 ± 114.7 + 194.0 ± 20.4 + 1949.3 ± 107.1 + NOTEARS 21.2 ± 11.5 - * 15.5 ± 9.9 + * 13.3 ± 4.3 ++ 41.3 ± 11.5 ++ 186.8 ± 83.0 --276.9 ± 92.1 - 718.4 ± 170.4 --1105.9 ± 250.1 ++ CAM 3.0 ± 2.2 - 8.1 ± 6.3 - 6.2 ± 5.5 - 28.5 ± 21.5 + 10.0 ± 4.6 - 44.2 ± 32.1 - 46.6 ± 9.5 - 882.5 ± 186.5

Table content: <table><tr><td>#Nodes</td><td>10</td><td/><td/><td/><td>50</td><td/><td/><td/></tr><tr><td colspan="2">Graph Type ER1</td><td/><td>ER4</td><td/><td>ER1</td><td/><td>ER4</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td><td>SHD</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>GraN-</td><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Results for real and pseudo real data sets with hyperparameter search

Table content: <table><tr><td>Data Type</td><td colspan="3">Protein signaling data set</td><td colspan="2">SynTReN -20 nodes</td><td/></tr><tr><td>Metrics</td><td>SHD</td><td colspan="2">SHD-C SID</td><td>SHD</td><td>SHD-C</td><td>SID</td></tr><tr><td>Method</td><td/><td/><td/><td/><td/><td/></tr><tr><td>GraN-DAG</td><td colspan="2">12.0 + 9.0 +</td><td>48.0 -</td><td>41.2 ± 9.6 -</td><td>43.7 ± 8.3 -</td><td>144.3 ± 61.3 +</td></tr><tr><td colspan="3">GraN-DAG++ 14.0 -11.0 -</td><td>57.0 -</td><td>46.9 ± 14.9 -</td><td>49.5 ± 14.7 -</td><td>158.4 ± 61.5 -</td></tr><tr><td>DAG-GNN</td><td>16.0</td><td>14.0 +</td><td>59.0 -</td><td>32.2 ± 5.0 ++</td><td>32.3 ± 5.6 ++</td><td>194.2 ± 50.2 -</td></tr><tr><td>NOTEARS</td><td colspan="2">15.0 + 14.0 +</td><td>58.0 -</td><td>44.2 ± 27.5 ++</td><td>45.8 ± 27.7 ++</td><td>183.1 ± 48.4 --</td></tr><tr><td>CAM</td><td colspan="2">11.0 + 9.0</td><td>51.0 +</td><td colspan="3">101.7 ± 37.2 --105.6 ± 36.6 --111.5 ± 25.3 ++</td></tr><tr><td>GSF</td><td colspan="2">20.0 -14.0 -</td><td>[37.0 +</td><td>27.8 ± 5.4 ++  *</td><td>27.8 ± 5.4 ++  *</td><td>[207.6 ± 55.4 -- *</td></tr><tr><td/><td/><td/><td>60.0]</td><td/><td/><td>209.6 ± 59.1]</td></tr><tr><td>GES</td><td colspan="2">47.0 -50.0 -</td><td>[37.0 +</td><td/><td/><td/></tr></table>

Peer review (iclr2020):
################## Response to rebuttal ################## I would like to thank the authors for their diligent rebuttal. Overall, I consider that the issues raised in my review have been mostly addressed. In 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. In 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. # Minor points / additional questions The 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. # Typos Section 3.6: … GraN-DAG performs fewer iteration(s) than NOTEARS… Section 5: GraN-DAGand DAG-GNN -was- were not designed with… Appendix A.1: … graphs can be roughly -halfed- halved when executed… Appendix A.5: … graphs were generated using the -BarabsiAlbert- Barabási-Albert model … Appendix A.9: … to selecting hyperparameters in practice consist(s) in ... ################## Original review ################## # High-level assessment The manuscript presents an alternative to the method in [2] to extend the ideas of [1] to nonlinear DGMs using neural networks. From 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 ( runtime per SGD step), I believe that addressing those issues (if indeed present) could be left for future work. The 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. In 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 ( , nonlinear and non-additive ground-truth models) while presenting no results for other settings commonly encountered in structure learning applications, such as high-dimensional problems or ground-truth models given by linear or generalized linear SCMs. Taking 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. # Summary In 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. The 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. The 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. In a nutshell, the authors accomplish this by: 1. Flexibly modelling the conditional distribution of each variable as a nonlinear function of (possibly) all other variables using an MLP, i.e. , with denoting any probability density differentiable w.r.t. its parameters and being a neural network with parameters . 2. Deriving a function such that is a sufficient condition for the conditional distribution of to be independent of . 3. Using , which is differentiable w.r.t. for most values of , to define a “connectivity” matrix to be directly plugged into the acyclicity criterion in [1]. 4. 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 regularization to induce sparsity. 5. 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. The 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 ( , ) 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. # Major points / suggestions 1. 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 matrix exponential evaluations being necessary. As 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 matrix exponential evaluations. To 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. 2. All experimental results reported in the manuscript use only data generation models with nonlinear and non-additive relationships between variables. This 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. While 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. In 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”. 3. Parallel to the previous point, all experimental results in the manuscript so far concern somewhat low-dimensional data and relatively abundant data . This 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)]. Also, 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). To 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 ( ) and (ii) studying the stability of the DAG w.r.t. different training runs. 4. All experimental results currently focus heavily on accuracy, while the computational aspect of the problem is largely secondary. To 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). # Minor points / suggestions 1. 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. 2. 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. # References [1] Zheng, Xun, et al. *DAGs with no tears: Continuous optimization for structure learning.* *Advances in Neural Information Processing Systems*. 2018. [2] Yu, Yue, et al. *DAG-GNN: DAG Structure Learning with Graph Neural Networks.* *International Conference on Machine Learning*. 2019. [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.