Dagma-DCE: Interpretable, Non-Parametric Differentiable Causal Discovery

We introduce Dagma-DCE, an interpretable and model-agnostic scheme for differentiable causal discovery. Current non- or over-parametric methods in differentiable causal discovery use opaque proxies of “independence” to justify the inclusion or exclusion of a causal relationship. We show theoretically and empirically that these proxies may be arbitrarily different than the actual causal strength. Juxtaposed with existing differentiable causal discovery algorithms, Dagma-DCE uses an interpretable measure of causal strength to define weighted adjacency matrices. In a number of simulated datasets, we show our method achieves state-of-the-art level performance. We additionally show that Dagma-DCE allows for principled thresholding and sparsity penalties by domain-experts. The code for our method is available open-source at https://github.com/DanWaxman/DAGMA-DCE, and can easily be adapted to arbitrary differentiable models.


I. INTRODUCTION
The discovery of causal relationships is the subject of many scientific questions, in medicine, economics, the social sciences, and more, and is therefore a question of utmost importance in science and engineering.Given an observational signal, one may therefore be interested in learning the causal structure underlying the data.Under Pearl's formalism of causality [1], the causal structure is captured by a causal graph, which is a directed acyclic graph (DAG) with edges denoting direct causation of one variable to another.
A core principle of causal inference is the notion of intervention, where one can interact with the system that produces signals.This is in contrast to causal inference from observational data, where we can only observe the outputs of a system, as is typically the case in signal processing or machine learning.Unfortunately, causal discovery using observational data is generally ill-posed and requires strong assumptions about the system we observe [2, pp. 135].Different methods may use different assumptions, where the validity of the assumptions may be extremely domain-or problem-specific.Even when assumptions are made to ensure (at least partial) identifiability of the causal graph, causal discovery can be extremely expensive, in terms of both the computation and the amount of required data.
One major class of methods, referred to as constraintbased methods, is based on the observation that a causal graph encodes a set of conditional independencies (CIs).By assuming that all CIs are encoded in the graph, an assumption known as causal faithfulness, constraint-based methods utilize statistical CI tests to determine which causal graphs are consistent with the data.
While some of the earliest successful causal discovery algorithms are constrained-based, such as the PC Algorithm and Fast Causal Inference [3], they can be very costly.Namely, they are both expensive computationally, scaling exponentially with the number of variables in the worst case, and they are expensive in terms of data.Additionally, CI tests tend not to be very powerful, and constrained-based causal discovery methods require conditional data for many different values and combinations of random variables.Moreover, causal faithfulness is not a testable assumption in general, though some weaker versions are [4].
An alternative approach to causal discovery using observational data is to use score-based methods.These methods attempt to find a causal model M and its associated causal graph G, which maximizes a score function S(M, D) for observational data D = {x 1 , x 2 , . . ., x N } [2, pp. 148].For example, one can search for the causal model which maximizes the likelihood or minimizes the mean square error (MSE).In practice, one typically optimizes a version of these objectives which is penalized for model complexity, such as the Bayesian information criterion.For certain classes of causal models, this approach is simple and effective -indeed, early attempts at score-based causal discovery were made contemporaneously to early constraint-based approaches [5].Despite their simplicity, score-based methods suffer from the super-exponential search space of DAGs with d variables, as well as identifiability issues [2, pp. 149].
Fundamentally troubling to score-based methods is the difficulty of optimization in discrete spaces, which tend to be combinatorially more difficult than continuous objectives.To address this, NOTEARS (short for Noncombinatorial Optimization via Trace Exponential and Augmented lagRangian for Structure learning) [6] proposes a novel algebraic constraint on adjacency matrices, such that the constraint is satisfied if and only if the adjacency matrix A represents a DAG.NOTEARS then performs constrained continuous optimization on a penalized MSE objective to learn a causal graph without explicitly exploring the discrete space of all DAGs.
NOTEARS was originally defined for linear models, with the adjacency matrix entries being coefficients of the corresponding linear functions.NOTEARS+ [7] extends NOTEARS to nonlinear functions and achieves state-of-the-art performance on simulated data.In lieu of linear coefficients, NOTEARS+ uses opaque proxies to define A, providing advantages in optimization but at the cost of interpretability and generality of the results.
Several methods have provided alternative acyclicity constraints which are faster to compute or are otherwise preferable in optimization.One striking example is DAGMA (short for DAGs via M-matrices with Acyclicity) [8], which in practice can optimize up to an order of magnitude faster than NOTEARS.
Our contribution is to provide a novel learning objective that results in obtaining interpretable weighted causal graphs.Our method has several important strengths.It is model-agnostic in nature, with the ability to use any twice-differentiable model trainable with gradient descent.Furthermore, the sparsity penalty used in the scoring function S(•, •) is also model-agnostic.Lastly, we can interpret the causal graphs with the differentiable causal effect (DCE) metric [9].We base our optimization on the method of DAGMA, leading us to name our method DAGMA-DCE.
The rest of this paper is structured as follows: in Section II, we review score-based causal discovery and the current stateof-the-art in differentiable causal discovery.In Section III, we show that the adjacency matrix returned by NOTEARS and DAGMA is not interpretable, first by a theoretical argument and then by a toy example.We propose a new graph definition and optimization problem in Section IV, forming DAGMA-DCE, and discuss its interpretability in Section V. Finally, we run experiments on synthetic data in Section VI showing that DAGMA-DCE performs competitively with other methods and returns legitimately different adjacency matrices than DAGMA.Concluding remarks are provided in Section VII.

II. BACKGROUND
We begin by reviewing some concepts of causal inference, namely structural causal models (SCMs), which formalize causality.We then introduce differentiable causal discovery methods, which utilize continuous optimization for scorebased causal discovery.

A. Structural Causal Models and Causal Discovery
SCMs provide a formalism for causality, where causal relationships are described by a DAG and accompanying functional assignments.To make this more clear, let x = [x 1 , . . . ,x d ] be d-dimensional random vector, and G be a DAG whose vertices are the components of x.Then, according to G, an SCM describes each component x j as a function of its parents Pa(x j ) and independent noise z j , f j (Pa(x j ), z j ), i.e., x j := f j (Pa(x j ), z j ) , j = 1, . . ., d. ( We use the notation := to reinforce the somewhat philosophical point that these are assignments, and not merely equivalence [2].We will frequently refer to the SCM of a system as simply the "causal model" M.
The goal of structure learning is then to identify the SCM of a system, given data D. If one cannot interact with the system, then we must make assumptions on the form of each f j and the distribution of the noise z j to uniquely determine the SCM [2, pp. 139].These assumptions can be parametric, for example, assuming each f j is linear and each z j is non-Gaussian [10], or they can be non-parametric, for example, assuming the f j are strictly nonlinear with the z j being additive and Gaussian [11,Corollary 31].Once we make such a restriction to our model, we call it identifiable.
After establishing identifiability, there are several families of techniques for learning the SCM from observational data.Here, we focus on score-based methods, though the interested reader is pointed to [12] for a recent, comprehensive review.In score-based causal discovery, we find a causal model M which maximizes a score function S(M, D).In practice, we parameterize models by some vector θ and then find the model M θ which optimizes S(•, •): The main drawback of naïve score-based methods is that they have to search over the space of DAGs with d nodes, which grows extremely quickly.For d = 5, this number is already 29,281; for d = 14, it is already greater than 10 36 [13, A003024].Therefore, for even small graphs, one must resort to greedy or approximate techniques [2, pp. 150 ( The earliest differentiable causal discovery algorithm with a DAG constraint is NOTEARS [6], which uses a trace of a nonlinear function on the adjacency matrix A θ , i.e., where ⊙ is the Hadamard (elementwise) product.The motivation for this constraint is that powers of the adjacency matrix A k θ reflect how many closed walks of length k are present in the corresponding graph G θ .Since DAGs have no closed walks of any length, h(A θ ) is 0 precisely when A θ is a DAG.
Quickly, several simplifications on the same constraint were derived.For example, one need not calculate the entire matrix exponential, with a polynomial sufficing instead [14].
The alternative constraint that we use is DAGMA [8]; with s > 0, the DAGMA constraint function is where 1 d is the d × d identity matrix.This constraint is motivated by the fact that A θ describes a DAG if and only if it is nilpotent.While the constraint given by Eq. ( 4) is, in general, a relaxation of nilpotency, the authors of DAGMA characterize the set of matrices for which zeroness of Eq. ( 4) is equivalent to acyclicity, called M-matrices.They then show that the set of M-matrices includes all DAGs, and that the gradient of h DAGMA (•) is sufficiently well-behaved, to enable constrained optimization within M-matrices.
It remains to specify the model M θ , the score function S, and how to derive G θ from M θ .The approach introduced in NOTEARS+ and used in DAGMA1 is to define where f j is the function x → x j according to M θ .To guarantee Eq. ( 5) is well-defined, it is required that ∥∂ i f j ∥ L 2 is finite.Accordingly, the NOTEARS+ authors restrict M θ to a space of functions called the Sobolev space H 1 (R d ).
However, the NOTEARS+ authors argue that for causal discovery, the zeroness of A θ ij is what matters, and so a condition equivalent to A θ ij = 0 is derived for several classes of models.For example, in a multilayer perception (MLP) network, the L 2 norm of the first layer in the neural network is chosen.The score function is defined to be a sparsity-penalized MSE 2 , where λ 1 is a hyperparameter and g(M θ ) defines a sparsity penalty.For MLPs, NOTEARS+ and DAGMA use the ℓ 1 norm of the first layer MLP weights.
In practice, it is very difficult to exactly follow the constraint in Eq. ( 3).For this reason, as well as to avoid false positives, an additional thresholding step is performed, where A ij is set to 0 if it is less than some predetermined hyperparameter C.

III. NON-INTERPRETABILITY OF DAGMA
A key motivation for DAGMA-DCE is that existing methods utilize a non-interpretable weighted adjacency matrix A θ .In this section, we argue theoretically that the adjacency matrix of DAGMA is not interpretable and subsequently show empirically that DAGMA fails to return intuitivelyinterpretable results on linear SCMs.
For the purposes of this section, we restrict ourselves to models where f j is modeled with an MLP network, though our arguments can be adjusted for other cases.The theoretical section can safely be skipped, but shows explicitly that adjacency matrices derived from the first layer of a neural network may be quite different than those defined by derivatives.

A. Theoretical Argument
For the purposes of this section, let f j : R d → R be an MLP with weight matrices A (1) , . . ., A (M ) and activation function σ(•).Then the (M −1)-hidden layer MLP f j can be expressed as Despite the exact analytic result that differentiable causal discovery methods rely on thresholding small, nonzero values.The following lemma shows that when using sigmoid activation functions, as is done in NOTEARS+ and DAGMA, the two norms in Eq. ( 8) may be arbitrarily different.
Then for any δ, ϵ > 0, there exists an MLP f j with weight matrices A (1) , . . ., A (M ) such that ∥A (1) Proof: 2 DAGMA uses the log-likelihood instead for benchmarks, but we use the MSE to simplify comparisons, and avoid assumptions on the data-generating process.

VOLUME ,
As ( 7) is a repeated composition of smooth functions, it suffices to show the result for the 1-hidden layer MLP g(x) = A (2) σ A (1) x .
hj x j .

For notational clarity, let
hj x j .
Additionally, note that σ ′ (•) > 0. Then so long as hi is nonzero for some h, we have by chain rule that for some ϵ ′ > 0. Then by directly calculating the norm, we conclude that it is unbounded.Therefore, we can choose A (2) such that the function norm is arbitrarily large, and in particular, larger than δ.While this result may not be particularly "generic," in the sense that the resulting function f j may be unlikely to be learned in practice, it does show that ∂ i f j can be arbitrarily different than any nonzero A (1) .This situation only gets worse for the more popular rectified linear unit (ReLU) activation function, where ∥∂ i f j ∥ L 2 is arbitrary for a fixed f j .This is formalized in the following lemma: Lemma 2. Let σ(•) denote the ReLU activation function.Then for any s > 0 and A (1) , . . ., A (M ) , there exist matrices B (1) , . . ., B (M ) such that (1) x , and ∥B (1) Proof: Note that for any real z, the ReLU function is invariant to positive rescaling: σ(z) = σ(az)/a for any a > 0. This implies, more generally, that S −1 A (2) σ(SA (1) x) = A (2) σ(A (1) x), where S is a diagonal matrix [15].The lemma then immediately follows by setting B (1) = SA (1) and , with the appropriate choice of S.

B. Empirical Result with a Linear SCM
In the lemmas above, we have shown that the adjacency matrix of NOTEARS+ and DAGMA may be arbitrarily different than the L 2 function norm.Here, we illustrate this fact in practice by applying DAGMA to a linear structural equation model (SEM), and showing that the returned adjacency matrix does not correspond to the intuitive notion of causal strength.
For this purpose, we simulated an Erdös-Rényi (ER) graph with 10 random variables and an expectation of 20 edges, creating a DAG as in [6].Given a DAG, a SEM assigns functional relationships between a variable and its parents.In this case, a linear SEM is assigned by randomly sampling linear coefficients from U ((−2.0, −0.5) ∪ (0.5, 2.0)).
The use of a linear SEM allows for an obvious ground-truth weighted graph to compare to, namely the linear coefficients.We then ran DAGMA with each f j being modeled as an MLP with one hidden layer of 10 units and compared the resulting weighted adjacency matrix to the ground truth in Fig. 1a.For the purposes of this demonstration, we tested with λ 1 ∈ {0, 10 −4 , 10 −3 , 10 −2 , 10 −1 } and chose the value which minimized the Frobenius norm ∥A est − A true ∥ F to give a favorable comparison.This resulted in using λ 1 = 10 −2 , with all other hyperparameters set to their default values.
Two things are visually apparent: first, DAGMA learns the graph topology quite well, including one erroneous edge in an otherwise perfect prediction.Second, however, is that the weights recovered by DAGMA do not obviously correspond to the interpretable, linear coefficients, i.e., it cannot easily be interpreted as a "strength of causation."Even worse, it is not easy to predict when DAGMA is over-or under-estimating the strength of an edge.While we defer a definition of DAGMA-DCE and discussion of Fig. 1b to the subsequent sections, a view of Fig. 1b shows a much closer resemblance to linear coefficients.It is important to note that this does not indicate "errors" in DAGMA, but instead that the magnitude of the weighted graph is not defined to be the true "causal strength" -their very definition has lost meaning compared to the linear case.

IV. PROPOSED SOLUTION
In this section, we provide a definition of the learning objective for DAGMA-DCE, and subsequently discuss practical details of implementing such a scheme.In particular, we discuss how DAGMA-DCE is model-agnostic in both formulation and implementation, and we also explore simple methods for accelerating and pre-training models.

A. Definition of the Optimization Problem
An initial solution would be to estimate Eq. ( 5) directly; however, we argue H 1 (R d ) is the wrong space to define models over.Not only are many interesting functions outside of H 1 (R d ), but also simple functions, for example linear functions, do not have a finite Sobolev norm over all of R d .Additionally, from a data-centric perspective, there is little reason to care about the behavior of the model where there is a vanishingly small probability mass.
Instead, we choose to define the adjacency matrix as the L 2 norm of the derivative ∂ i f j with respect to the probability distribution of X, which we denote P X .While the theory of Sobolev spaces plays no further role in our work, to adapt the language of [7], models are defined over the weighted Sobolev space H 1 (R d , P X ) [16].This is simply to say that the models' function and first (weak) derivatives are squareintegrable with respect to P X .We then offer the alternative definition and Monte Carlo approximation where p X (x) is the density of P X with respect to the Lebesgue measure, and x n ∼ p X (x).By noting that the entries of A are the root-mean-square entries of the Jacobian matrix, it is then straightforward to calculate A ij for any differentiable model, with the above being a Monte Carlo approximation over empirical Jacobian matrices.This definition of A also allows for a principled sparsity penalty -indeed, if ∂ i f j is square-integrable with respect to P X , then Jensen's inequality implies that it is also integrable, which allows us to define and approximate Combining Eqs. ( 3), ( 4), ( 6), ( 9) and (10), we arrive at the optimization problem Similar to DAGMA, we use a central path optimization technique to solve the constrained optimization problem.
Our approach has the same convergence properties as NOTEARS and DAGMA, using similar assumptions.Namely, we assume causal sufficiency, i.e., a lack of latent confounders, and identifiability, which ensures that the optima of the score are the true causal graph.There are many options for ensuring identifiability, for example, thrice-differentiable strictlynonlinear additive noise models suffice [11].Additionally, while it is not necessary for the optimization problem to be well-defined, the interpretability of DAGMA-DCE relies on once-differentiability of the underlying f j .Namely, this excludes discrete or mixed-type data, since how a model interpolates functions between discrete labels is arbitrary.

B. DAGMA-DCE is Model Agnostic
The derivatives in Eq. ( 9) may be available analytically for many models.For the simplest case, take M θ to be linear, i.e., Then the linear coefficients β ij are the derivatives used in Eq. ( 9).Note that Eq. ( 11) then becomes the problem of linear DAGMA, meaning DAGMA-DCE is a legitimate generalization.
Simple (semi-)analytic expressions are also available for additive models: in which case Analytic expressions of coefficients are even available for several classes of non-parametric models.For example, by the representer theorem [17, pp. 132], a Gaussian process (GP) posterior can be written in the form of Eq. ( 12), where g(x) depends on the kernel κ(•, •) and β j1 , . . ., β jn are easily calculated: Closed-form expressions for several common kernels, including the automatic relevance detection squared exponential (ARD-SE) and Matérn-3/2 kernels, are derived in [9].Furthermore, even when analytic expressions become cumbersome to derive by hand, the derivative can often be obtained via automatic differentiation (AD) tools.Among many other differentiable models, a notable case is with MLPs.It is also the case for extensions of NOTEARS+-like objectives to time series, for example the one-dimensional convolutional neural networks used in NTS-NOTEARS [18].

C. Practical Considerations
Several general comments are made in Section 4.1 of [8] which also hold for DAGMA-DCE, particularly concerning the optimization details.For example, the hyperparameter s can be changed freely, even between optimization iterations.This can affect performance, as s controls the volume of the feasible region.The authors of [8] recommend starting with s = 1 and slowly reducing its value.Additionally, different optimizers can be used to solve sub-problems along the central path; DAGMA uses Adam [19], which we also use.Finally, the use of central path methods also creates additional hyperparameters, including central path coefficients and weight initializations, which are discussed.
While backward passes over the Jacobian is indeed expensive and a drawback of the DAGMA-DCE, the Monte Carlo approximation can be effectively parallelized with respect to N , particularly on GPUs.Therefore, while scaling with d is inherently difficult, scaling with N is quite inexpensive.
In practice, we additionally see benefits from pre-training for DAGMA-DCE, namely by training a single iteration of DAGMA.Intuitively, while DAGMA and DAGMA-DCE are fundamentally different optimization problems, their underlying task is the same.Therefore, a single iteration of DAGMA serves as a good starting point, and saves costly backward passes over the Jacobian.
Finally, DAGMA-DCE has several hyperparameters, including the strength of L 1 regularization λ 1 and thresholds.In the next section, we discuss how the formulation of DAGMA-DCE makes the threshold hyperparameter interpretable.However, the issue of setting λ 1 remains.One option, explored in [7] and [18], is to optimize λ 1 with respect to some metric using cross-validation if labeled data/causal graph pairs exist.How to choose λ 1 when these labeled pairs are not available is an interesting task for future research.

V. INTERPRETABILITY OF DAGMA-DCE
The interpretability of the DAGMA-DCE approach stems from the observation that the partial derivatives, ∂ i f j , can be interpreted through the lens of causal strength.We begin this section with general notes on causal strengths, followed by the differential causal effect, which forms the basis of interpretability for DAGMA-DCE.We revisit the empirical example of Section III-B and show that DAGMA-DCE recovers the intuitive result.Finally, we discuss how thresholding in DAGMA-DCE is a fundamentally interpretable operation, allowing the incorporation of prior knowledge from experts and decision-makers.

A. Measuring Causal Strength
The quantification of causal strength is a highly nontrivial task, with several competing definitions [20].Despite this, we find it helpful to ground our reasoning in the linear case, where an unambiguous notion of causal strength is available.In particular, the linear model coefficients quantify exactly how a change in one variable will induce changes in its children.Consider a linear SCM with coefficients given by a matrix [B] ij = β ij and additive noise z j : Then the entry β ij represents the direct influence of x i to x j .Note that B can be viewed as a causal graph, provided that it is acyclic, or equivalently, nilpotent.Occasionally, these direct coefficients summarize all causal information that x i holds about x j .Indeed, if B is viewed as a graph, this is the case if the only path x i and x j is the direct one, in which case B ij is the total causal strength.Otherwise, in general, we must sum across all paths in the graph, and multiply the coefficients along each path.This can be expressed compactly by a geometric sum, Causal sufficiency of the system then tells us that the entries T ij of the matrix T are sufficient to describe how any change to x i will affect x j , and thus, T can be useful to ask counterfactual questions or to understand the impact of various causes to a single effect.

B. Differential Causal Effect and DAGMA-DCE
The DCE provides a nonlinear version of the linear model causal strength by considering the local behavior of the model [9].The Jacobian matrix, J f , replaces the role of B, describing how changes to a parent affect the child node when there are no other paths in the system.Similarly, the geometric sum of the Jacobian also defines the total causal effect from source to effect, for small changes.
Relating to DAGMA-DCE, the adjacency matrix defined in (9) can be interpreted as the root-mean-squared DCEas a result, A ij ̸ = 0 if and only if the DCE, ∂ i f j , is not identically zero.The magnitude of A ij quantifies how strong the interaction is, in terms of the energy of the DCE functions.

C. Empirical Results with a Linear SCM
We quickly revisit the experiment of Section III-B, where DAGMA was fit to a linear SCM, and its weighted adjacency matrix compared to the ground truth.The setup is nearly identical, with the same parameterization of each function f j .The main difference is that we do not tune the sparsity penalty hyperparameter, instead we simply set λ 1 = 0. Like DAGMA, the DAGMA-DCE method obtains an accurate estimate of the graph topology, recovering the underlying graph perfectly.Unlike DAGMA, we additionally recover the "intuitive" notion of causal strength.

D. Thresholding in DAGMA-DCE
An ad-hoc component of NOTEARS+-like methodologies is the thresholding step.Indeed, across different datasets, the thresholding hyperparameter can have large consequences to accuracy, as shown in [7].However, many conceivable causal inference tasks will not have any validation set to calibrate the thresholding parameter with, meaning the user will have to make a heuristic estimate.
On the contrary, DAGMA-DCE's weighted adjacency matrix entries are interpretable, so the choice of a threshold carries direct meaning to the modeler or expert.This decision can even be made on an edge-by-edge basis, or across (a possibly rescaled version of) A.
For example, one may decide to normalize the entry A ij of by the marginal variance of x j for the purposes of thresholding, reminiscent of ANOVA [21].An alternative scheme could be normalizing by the sum of the jth column; while subadditivity of norms means this does not actually scale by the total derivative ∂f j , it can be viewed as scaling by the total causal effect.Finally, since the entries of A are now unconstrained, one could attempt thresholding methods developed for Granger causality, for example via clustering [22] or denoising [23].How this thresholding step should be interpreted and executed forms an interesting question for future work.

VI. EXPERIMENTS
We perform several experiments similar to those conducted by the authors of DAGMA [8], providing comparisons to DAGMA and NOTEARS+.We begin with a brief discussion of metrics used to evaluate results, followed by experimental details, results, and discussion.

A. Metrics
Evaluating results in causal discovery is complex in the sense that several different metrics exist with different goals and calculations.The two most basic metrics we report are the structural Hamming distance (SHD), which calculates the number of incorrect edges a graph contains, and the F 1 score, which is the geometric mean of precision and recall.
Additionally, we include the structural intervention distance (SID) [24], which is arguably more well-suited to causal inference.The SID is motivated by the idea that, when calculating causal effects via "adjustment sets," some errors are more severe than others.For example, including spurious relations does not affect causal effects, while reversing the order of an edge does.The SID calculates the number of edges i → j such that the adjustment set corresponding to the calculation of the causal effect of x i on x j is not valid.
Finally, we report the time elapsed in minutes in each of our synthetic benchmarks.

B. Experimental Details
We replicate a common synthetic dataset in the causal discovery literature, generating ER-4 DAGs according to the procedure of [7], with functional relationships assigned in one of two ways.The first way is using additive GPs: where each g ij is sampled from an RBF GP with lengthscale 1.The other is with a random MLP network with hidden size 100 and a sigmoid activation.Crucially, all weights are samples from U ((−2.0, −0.5) ∪ (0.5, 2.0)).
As noted in Section IV-C, the hyperparameters of each method may generally be tuned with cross-validation.Nonetheless, to facilitate a fair comparison between methods, we use a set of "defaults" as in [7], [8], [18].For DAGMA and NOTEARS+, we use the hyperparameters noted in Appendix C.2.2 of [8].However, since the sparsity penalty of DAGMA-DCE is fundamentally different than that of other methods, it is difficult to find a "default" setting for a fair comparison.We chose λ 1 = 3.5 × 10 −2 , which provided a sparsity penalty of similar magnitude in the MLP experiments.This is imperfect, but allows for the initial score functions of DAGMA and DAGMA-DCE to be similar without unfairly optimizing λ 1 .Similarly, we used 0.25 as thresholds for comparison, which provides a similar true positive rate in the MLP experiments.Additionally, due to the added cost of DAGMA-DCE, we set the maximum number of iterations to 10% of the corresponding values in DAGMA.
All experiments were run on a NVIDIA Titan RTX GPU.The authors' implementations of NOTEARS+ 3 and DAGMA 4 were used, with slight modifications to allow GPU acceleration.All methods (including DAGMA-DCE) used the same MLP architecture used in the experiments of [7], [8].

C. Results & Discussion
Our results show that DAGMA-DCE achieves the stateof-the-art on common synthetic datasets: in additive GPs (Fig. 2a), the median SID is universally lower than DAGMA or NOTEARS+.Similarly, the F 1 score is higher for all but D = 50 nodes, and the SHD is significantly lower for d ≤ 30 nodes.While DAGMA-DCE is indeed more computationally expensive than DAGMA, it is still substantially faster than NOTEARS across all graph sizes.In the MLP experiments (Fig. 2b), our results are somewhat less impressive, but still competitive.This is despite what we argue is a generating process that gives a distinct advantage to DAGMA and NOTEARS+: by recalling that these methods perform thresholding on the norm of layer weights, sampling weights from U ((−2.0, −0.5) ∪ (0.5, 2.0)) builds-in identifiability in a perhaps unrealistic manner.
Finally, without asserting what the "correct" magnitude of weighted adjacency matrices is, it is interesting to note that DAGMA and DAGMA-DCE return adjacency matrices in which the relative magnitudes differ considerably.We measure this with Kendall's τ -b and the Spearman correlation coefficient, which measure rank correlation.These both take on values in [−1, 1], where a value of ±1 indicates complete agreement or disagreement, and a value of 0 indicates that orderings are unrelated.The resulting values for our synthetic experiments can be found in Table 1.

VII. CONCLUSION
We reformulated the objectives of NOTEARS+ and DAGMA in order to extract a principled, interpretable weighted adjacency matrix.With this comes principled, interpretable notions of an ℓ 1 sparsity penalty, and thresholding the resulting graph.A drawback of this method is its computational complexity, requiring a backward pass over the Jacobian; however, empirical results indicate it is still much faster than some other competing methods that use an augmented Lagrangian scheme for optimization.Furthermore, we were able to achieve stateof-the-art performance on synthetic datasets with small-tomedium-sized graphs.Interesting directions for future work include applications to real datasets and a formal workflow for determining sparsity and thresholding hyperparameters.

FIGURE 1 :
FIGURE 1: The difference between the magnitude of the true derivatives in a linear causal model to the magnitude of the weighted graph in DAGMA for a random 10 × 10 Erdös-Rényi directed graph with 20 expected edges.Gray boxes surrounding each cell denote the magnitude of the ground-truth linear coefficient.

FIGURE 2 :
FIGURE 2: Resulting SID (top left), SHD (top right), F 1 Score (bottom left), and time elapsed (bottom right) for random data generated from (a) the ER-4 GP-additive model and (b) the ER-4 MLP model, as detailed in Section B. Boxes show the median and quartiles across T = 10 trials for DAGMA and DAGMA-DCE, and T = 5 trials for NOTEARS, with whiskers showing the minimum and maximum values.

TABLE 1 :
Rank correlations between adjacency matrix entries of DAGMA and DAGMA-DCE across all trials for synthetic datasets, given as the mean plus-minus one standard deviation.