DAGMA: Learning DAGs via M-matrices and a Log-Determinant Acyclicity Characterization
Authors: Kevin Bello
Bryon Aragam
Pradeep Kumar Ravikumar
Venue: NeurIPS 2022 Accept (2022)
Keywords: structure learning
causal discovery
graphical models
continuous optimization
non-convex optimization
log-determinant
directed acyclic graphs

Abstract: The combinatorial problem of learning directed acyclic graphs (DAGs) from data was recently framed as a purely continuous optimization problem by leveraging a differentiable acyclicity characterization of DAGs based on the trace of a matrix exponential function. Existing acyclicity characterizations are based on the idea that powers of an adjacency matrix contain information about walks and cycles. In this work, we propose a new acyclicity characterization based on the log-determinant (log-det) function, which leverages the nilpotency property of DAGs. To deal with the inherent asymmetries of a DAG, we relate the domain of our log-det characterization to the set of $\textit{M-matrices}$, which is a key difference to the classical log-det function defined over the cone of positive definite matrices.
Similar to acyclicity functions previously proposed, our characterization is also exact and differentiable. However, when compared to existing characterizations, our log-det function: (1) Is better at detecting large cycles; (2) Has better-behaved gradients; and (3) Its runtime is in practice about an order of magnitude faster. From the optimization side, we drop the typically used augmented Lagrangian scheme and propose DAGMA ($\textit{Directed Acyclic Graphs via M-matrices for Acyclicity}$), a method that resembles the central path for barrier methods. Each point in the central path of DAGMA is a solution to an unconstrained problem regularized by our log-det function, then we show that at the limit of the central path the solution is guaranteed to be a DAG. Finally, we provide extensive experiments for $\textit{linear}$ and $\textit{nonlinear}$ SEMs and show that our approach can reach large speed-ups and smaller structural Hamming distances against state-of-the-art methods. Code implementing the proposed method is open-source and publicly available at https://github.com/kevinsbello/dagma.

DAGMA: Learning DAGs via M-matrices and a Log-Determinant Acyclicity Characterization
Pradeep Ravikumar* 'oooh School of Business, University of Chicago, Chicago. IL 60637 Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA 15213
Kevin Bello't
Bryon Aragam! |
Abstract
The combinatorial problem of learning directed acyclic graphs (DAGs) from data was recently framed as a purely continuous optimization problem by leveraging a differentiable acyclicity characterization of DAGs based on the trace of matrix exponential function. Existing acyclicity characterizations are based on the idea that powers of an adjacency matrix contain information about walks and cycles. In this work, we propose new acyclicity characterization based on the log- log-determinant (log-det) function, which leverages the nilpotency property of DAGs. To deal with the inherent asymmetries of a DAG, we relate the domain of our log-det characterization the set of M-matrices, which is key difference the classical log-det function defined over the cone of positive definite matrices. Similar to acyclicity functions previously proposed, our characterization is also exact and differentiable. However, when compared to existing characterizations, our log-det function: (1) better at detecting large cycles; (2) Has petter-behaveds gradients; and (3) Its runtime is in practice about an order of magnitude faster. From the optimization side, we drop the typically used augmented Lagrangian scheme and propose DAGMA (Directed Acyclic Graphs via M-matrices for Acyclicity), a method that resembles the central path for barrier methods. Each point in the central path of DAGMA is solution to an unconstrained problem regularized by our log-det function, then we show that at the limit the central path the solution is guaranteed to be DAG. Finally, we provide extensive experiments for linear and nonlinear SEMs and show that our approach can reach large speed- ups and smaller structural Hamming distances against state-of ttte-of-thhe-r art methods. Code implementing the proposed method is open-source and publicly available at httpss//gihuu..oo/keeiiistt
Introduction
Structural equation models (SEMs) [44] are a standard modeling tool in several fields such as economics, social sciences, genetics, and causal inference, name few. Under this framework, in its general form, the value of each variable in the model is assigned by a general nonlinear, nonparametrie function that takes as input the values of other variables in the model, thus, every SEM can be associated to a graphical model. In particular, we will consider graphical models that are directed acyclic graphs (DAGs).
long-standing and active research area deals with the problem of learning the graphical structure (DAG) given passively observed data ka.k.a. causal discovery). Computationally, this problem is well-known to he NP-hard in general [10, 12]. mainly due to the combinatorial nature of the space of DAGs. In this work, we will follow score-based approach, popular learning framework where the goal is to find DAG that minimizes a given score [21, 11]. Recently, Zheng et al. [58] proposed an exact smooth nonconvex characterization of acyclicity which opened the door solving the originally
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

combinatorial problem via set of tools that work in the fully continuous regime (e.g., gradient-based methods).
Let W E Rddd be a weighted adjacency matrix of a graph G of d nodes, and let W o W denote the Hadamard product. The acyclicity function introduced by Zheng et al. [58] is defined as Axppm (W) Tr(cWoW) d, where Tr denotes the trace of matrix, and it was shown that hexpm (W) 0 if and only if W corresponds to a DAG. A follow-up work [56] proposed another acyclicity function, which can be computed slightly faster, defined Npoly is the identity matrix. It was similarly shown by Yu et al. [56] that Apoly =0 if and only if W corresponds to a DAG. While seemingly different, both hexpm and Apoly
TT((( + J W))) d. where
functions of the
form Tr(En=o Gp(W o WY*) d for some Ck > 0, as noted by Wei et al. [55]. To the best of our knowledge, all subsequent work [e.g.. 59, 26, 60, 37, 38, 36, 25, 41, to name a few] that has built upon the idea of using a continuous acyclicity characterization have used either hexpm or /plly. or some acyclicity characterization in the form of trace of a sum of matrix powers. The latter should come as no surprise, after all. nonzero diagonal entry of the matrix power (W o W) reveals the existence of closed walk of length k in W.
Contributions, In this work, we propose a new acyclicity function that, as hexpm and hooly. is both an exact and smooth acyclicity characterization but that also possesses several advantages when compared to hexpm and Aplyy Specifically, we make the following set of contributions:
1. We propose a novel acyclicity characterization based on the log-determinant (log-det) function (see Theorem and Section 3). In contrast the classical log-det function defined over the cone of positive definite matrices, we define the domain of our log-det function be the set of M-matrices due to the inherent asymmetries of DAGs. To our knowledge, we are the first to connect the notion of M-matrices to acyclicity and structure learning for DAGS.
2.
We provide detailed study of the properties of our log-det characterization in Section 3.1. First, we establish the similarities of our log-det function to other existing functions such hexpm and hpplly- Second, we formally argue why these functions can be regarded as acyclicity regularizers, similar in spirit to the classical l and l regularizers. Third, we show that our log-det function is an inver function, i.c., all stationary points are global minimum, moreover, these stationary points correspond to DAGS.
3. In Section 3.2, we present three arguments as to why our log-det characterization could be preferred over other existing acyclicity functions. Briefly, our log-det function is better at detecting large cycles, has better behaved gradients, and can be computed in about an order of magnitude faster than hexpm and hpoly-
4. Motivated by the properties of our log-det function, in Section 4, we present DAGMA (Directed Acyclic Gruphs via M-matrices for Acyclicity), a method that resembles the widely known central path approach for barrier methods [39]. We show that, at the limit of the central path, the solution is guaranteed to be a DAG. In contrast to the commonly adopted augmented Lagrangian scheme (originally proposed in [58]) each point in the central path of DAGMA is a solution to an unconstrained problem regularized by our log-det function.
5. Finally, in Section 5 and Appendix C. we provide extensive experiments for linear and nonlinear SEMs under different score functions (both least squares and log-likeline likelihoodd, where we show that DAGMA is capable of obtaining DAGs with better accuracy, i.e., lower structural Hamming distance (SHD), in much faster way than the state-of-the-arr.
1.1
Related work
The vast majority of methods for learning DAGS can be categorized into two groups: constraint-basee algorithms, which rely on conditional independence tests; and score-based algorithms, which focus on finding a DAG that minimizes a given score/loss function. We briefly mention classical constraint- based methods as we follow a score-based approach. [50] developed the PC algorithm, a popular general method that leams the Markov equivalence class. Other algorithms such as [52, 29] are based on local Markov boundary search. Finally, hybrid approaches that combine constraint-based learning with based learning, such as [53, 15].
2

---

In the line of score-based methods, popular score functions include BDeu [211, BIC [31], and MDL [7]. Works that study linear Gaussian SEMs include [2. 3, 16, 17, 32,45,, and for linear non- Gaussian SEMs [27, 491. For nonlinear SEMs, we note works on additive models [9, 13, 54], additive noise models [23, 46, 35], generalized linear models [43, 42, 19], and general nonlinear SEMs [34, 18].
More closely related to our work is the line of work built on the nonconvex continuous framework of Zheng et al. [58]. such as, [59, 26, 60, 37, 36, 25, 41]. In contrast to our work, all of the aforementioned methods rely on the nonconvex acyclicity functions Aexpm or hpooy. and with the exception of [37], all of these works also use the augmented Lagrangian scheme. Finally, the NoCurl method [57] also departs from using hexpms although no other acyclicity constraint is proposed. Two immediate distinctions can be made to our work. First, we propose a novel acyclicity function based on the log-det function which we show to be prefereable to Rexpm and Apoly Second, we drop the commonly adopted augmented I Lagrangian scheme to solve the constrained problem and instead follow central path approach to leverage the barrier property of our log-det function.
Remark 1. To avoid confusion, we also note that in the GOLEM method Ng et al. [37] the score includes a log-determinant function of the form log det(T - W)l which stems from the Gaussian log-likelihood. While this expression is zero if W corresponds to a DAG, it is not an exact acyclicity characterization (i.e. log det((T W) 0 does not imply W is a DAG). By contrast, the use of M-matrices in our work is crucial to translating the log-det function into a valid acyclicity regularizer. Moreover, it is not obvious how extend GOLEM to arbitrary score functions, as their analysis is specific to the Gaussian likelihood function.
2 Notation and background
Notation. We use [d] to denote the set of integers {1 (1...d}. dd. For a square matrix A, we use (x to denote its i-th minimum eigenvalue, and use p(AA) to denote its spectral radius. Also, we use Tr( (A). and (A) to denote the trace and determinant of A. For matrices A, B, we let A o B represent the element- woinent-wise or Hadamard product, moreover, the expression A 2 B is entrywise, i.e., Aij 2 B Then, we say that matrix A is nonnegative whenever 2 O. For a complex number a + bi, we let R(a + bi) = a denote its real part. We use ||-1Ip to denote the vector f--norm, and ||- is the T--norrm on functions. Lastly, -> j and i ~ j represent an edge from to j and a directed walk from i to j. respectvvely.
Let X (X, Xa) be a d-dimensional random vector. In its general form, a (nonparametriie) structural equation model (SEM) consists of a set of equations of the form:
X₃= fX, 2j), Vj E [d].
where each f R -> R is nonlinear nonparametric function, and Z, is an exogenous variable representing errors due to omitted factors. We consider the Markovian model, which assumes that each Zj is an independent random variable. Note that each Sj depends only on a subset of X (i.e., the parents of X,) and Z,; nonetheless, to simplify notation we ensure that each ʃ is defined on the same space. Then f = (f1..... fa) induces a graphical structure, where we focus on directed acyclic graphs.
For any joint distribution over Z (Z1..... Za). the functions f, define a joint distribution PXX) over the observed data, and a graph G(J) via the dependencies in each fz. Then, our goal is to learn G(f) given n i.i.d. samples from P.(X.. In basre-based learning, given a data matrix X xt] E Raxdd we define a score function Q(;; X) to measure the "quality" of candidate SEM as follows: X) = Σÿ= ( f, (X)). where we adopt the convention that (X) e v) 1 E (u )2)² or the log-likelihood function, often augmented with a penalty such as BIC or l, Given the score function Q and a family of functions F, we seek to find the f E F that minimizes the score, i.c.
[21 (.1,...,
R". Here loss can be any loss function such as least squares loss( u,
min X) subject to G(f) E DAGs. FEF
Similar to [59], we consider that each f, lives in a Sobolev space of square-integrable functions whose derivatives are also square integrable. Then, let 0xf, denote the partial derivative of f, w.r.t. Xx. it is easy to that fj is independent of X₃ if and only if [|0xflllz O. With this observation, we construct the matrix W (f) Idd with entries [W))s
[@.f, I1, which precisely encodes def
3

---

the graphical structure amongst the variables X,. That is G(f) E DAGs = W E DAGs, where W is interpreted as the usual weighted adjacency matrix.
In practice, f is replaced with a flexible family of parametrized functions such as deep neural networks, so that problem (2) is finite dimensional. Finally, note that model (1) includes several models as special cases, c.g., additive noise models, generalized linear models, additive models, polynomial regression, and index models. Previous work has studied the identifiability of several of these models, e.g., [23, 49, 46, 43, 42, 27]. In the sequel, we assume that the model is chosen such that the graph G(f) is uniquely defined from (2).
3 A new characterization of acyclicity via log-determinant and M-matrices
In this section, we present our acyclicity characterization and study its properties, To declutter notation, this section simply write W instead of W (f) to denote the weighted adjacency matrix of a graph; however, it should be clear that W depends on functions ʃJ as explained in the previous section.
We develop our characterization by first noting that for any nonnegative weighted adjacency matrix W. we have that W e DAGs if and only if W is nilpotent matrix, or equivalently, all the cigenvalues of W are zero, i.e., >x(W)= = 0, Vi E [d]. Then, for any W E Rdxdd we have the following obvious implications:
W DAGs + (W o W) E DAGs + 5 N.W o W) s, s,
(3)
d
II 5 (. o W) = det(sl W o W) = 5d
(4)
Implication (4) can be thought of as a relaxation of acyclicity in the sense that all DAGs satisfy (4), but not all W that satisfy (4) are DAGs. For example, let and W W [82]. then it is clear that det(sI- - W W) 1 and, thus, (4) is satisfied; however, clearly W is not DAG.
Thus, two immediate questions arise: (i) Does there exist domain for W such that (4) = (3)? (ii) If so, what is the dexcription of such domain? We answer (i) in the affirmative, and answer (ii) by relating the domain of W to the set of M-matrices, which is defined below.
Definition 1 (M-matrix)), ). An M-matrix is a matrix A E Iddxd of the form A A= sl - B, where B 2 0 and 8 > B)).
M-matrices were introduced by Ostrowski [40] and arise in variety of areas including input-output analysis in economics, linear complementarity problems in operations research, finite difference methods for partial differential equations, and Markov chains in stochastic processes. To the best of our knowledge, we are the first to connect the notion of M-matrices to graphical model structure learning through an acyclicity characterization.
The following proposition is immediate consequence of Definition 1:
Proposition 1 (Berman and Plemmons [5]). Let A E Rddd be an M-matrix, then:
(i) (A)( 0. all E [d]. (ii) A- exists and is nonnegative, i.e., A-1 2 0.
In the above, item (i) states that the eigenvalues of an M-matrix lie in the open right-half plane. Matrices which satisfy the latter property are also known as positive stable matrices. We thus have that M-matrices are special cases of positive stable matrices. It follows that the determinant of any M-matrix is positive. This fact will be used for defining hi our acyclicity characterization given in Theorem 1. Finally, the nonnegativity of the inverse from item (ii) will used to understand some properties of the gradient of h hjce (W).
We now define the domain over which (4) = (3) (see Theorem 1). For any s > 0, define
W = (W E R ddd s > p(W o W)}.
(5)
Moore precisely, we consider the definition of non-singular M-matrix, which is sufficient for the purposes of this work.
NNote that due to asymmetries is possible for an M-matrix to have complex eigenvalues. However, since we work with matrices with real entries, the complex eigenvalues come in conjugate pairs.

---

i.e. W is the set of real matrices whose entry-wise square given by W o W have spectral radius less than s. The following lemma lists some relevant properties of W7.
Lemma 1. Let V defined as in (5). Then, for all 0:
(i) DAGs C CW.
(it) W is path-connctedd
(ii) W C Wt for any L>
In Lemma 1, item (i) implies that the W we look for is in the interior of W" item (ii) indicates that we can find a path from any point in W to any DAG without leaving the set Wo; item (iii) shows that one can vary s to enlarge shrink the set W
Having defined the domain set W", we now define our acyclicity characterization hider : W - R. Recall that item (i) in Proposition implies that applying the logarithm function to the determinant of an M-matrix is always well defined, which motivates our following result.
Theorem (Loog-determinant characterization). Let 8 > O and let hider" W - R be defined as log det(sl I
ideer (W) def
W W) + d dolo s. Then, the following holds:
(i) det (W) 2 O, with hige (W)= = 0 if and only if W is a DAG.
(ii) Vhiiee (W) 2(s1 - W W W, with Vhia (W) 0 if and only if W is DAG.
3.1 Properties of hidec (W)
In this section, we list several properties of our acyclicity characterization. The first property we discuss is related to the entries of Vhjaet (W).
Lemma 2. For all j.j E [d]. Vhidet (W) Wijj = 0 if and only if W., 0 or there is no directed walk from j to i. Equivalently, Whhie (W) )kk F 0 if and only if the edge - is part of some cycle in W. Finally, whenever (V hidet (W) Olk, # O, we have that sign
(W)) sign(W..))
The lemma above characterizes the nonzero entries and their signs of the gradient of hfse This property formally offers a regularizer perspective for ha which we highlight next.
Remark 2 (A regularizer viewpoint). The function hia promotes small parameters values in much the same way that the classical e, and lg regularizers do. In contrast to the latter regularizers, hitet will only shrink the value of a parameter Wi. if and only if the edge (i,j) is part of some cycle in W, as prescribed by Lemma 2.
Recall that hexpm (W) = ToWr - d and hooy (W) = Tr((I+ + ! o W) - d. It was noted by [55] that acyclicity characterizations of the form Tr(_n= (W W)d) for cp > also have the property in Lemma 2. This implies that hexpm and Apoly hold the property above and can also be interpreted as acyclicity regularizers. We note that the interesting part here is that hides holds this property besides being different nature to hexpm and hpplyy-
Next, we state an important consequence of Lemma 2, which is related to the direction of Vhided (W).
Corollary 1. At any W E W", the negative gradient Whider (W) points towards the interior of Ws.
In optimization, the Hessian matrix plays important role it contains relevant information about saddle points and local extrema of a function, and is key to Newton-type methods. Another appealing property of Maet! (W) is that it has Hessian described by a simple closed-form expression.
Lemma 3. The Hessian of Meae (W), which resides in Rddd xd? is given by:
Vhifact (W)- 4 Diag(vec(W))) & TT) (W) + Diag(vec([V) ))),
where N (sI- W W/ denotes the Kronecker product, and Kddd the di de commutation matrix such that dddd vec(A) vec(A] for any d matrix A.
Here we note that among /eeppmy Appuly and heaeet only Mfaet has a tractable expression for the Hessian.
Furthermore, note that V2 (W) is indexed by vertex pairs so that the entry [V (W)) (W))(..)
8 corresponds to the second partial derivative aWe Using Lemma 3, we can characterize the nonzero entries, and their signs, of the Hessian of hidet
5

---

Corollary 2. The entries of the Hessian ha (W) are described as follows:
M (W) (k.0).(p-))
AW.LN N Wp, 4.P 4(WL)) 2N
if (k, + (p. 9), if (k, = (p.g),
where (s1 W o w) Moreover, an off-diagonal entry [V2 (W) (k..).......) is nonzero if and only if there exists cycle in W of the form q p 1 + k ~n q, and has sign equal to Lastly, a diagonal entry ( (1..).))) is nonzero if and only if there exists Wi.WWa.)) sign(Wi,
directed from k to I, and its sign is always positive.
Recall from Theorem that all DAGs attain the minimum value and are critical points of hjaes thus, DAGs are local (and global) minimum of Miet and the Hessian matrix evaluated a DAG must be positive (semi)definite Let us corroborate the latter, when W is a DAG, from Corollary 2 we have that the off- off-diagonal elements of Maee (W) are zero, while the diagonal entries nonnegative. That is, the Hessian is positive semidefinite whenever W is a DAG. This implies that all its stationary points are global minima: Such functions called invex [20, 30].
Corollary 3. Let s > O. Then, hide (W) is an invex function,, i.e., all its stationary points are global minima, and these correspond DAGs.
We note that even though hexpm and poly are also invex functions, they were not explicitly considered as invex functions before. In fact, it was noted in [58] that DAGs were global minima of hexpm but no characterization of its stationary points were given. Wei et al. [55] noted that DAGs were stationary points of hoxpm and hpoly but the notion of invexity was not explicitly stated.
Remark 3 (A dynamical system perspective). The importance of invexity here is that in the eyes of hey all DAGs are the same. That is, DAGs correspond to the set of attractors in hider and depending the initial condition, the system will converge different attractor. This offers the following viewpoint for the role of the score function Qf;;X) namely, "use the score Q to find a basin of attraction such that the force field of Mase will dictate the trajectory towards DAG that is equal or close to the ground-truth".
In Figure 1, we illustrate in a toy example the properties discussed in this subsection.
(a) Made (W)
(b) Contours of fiaes (W)
(c) Vector field of V (W)
Here clearly W is a DAG whenever one of wi or We (or both) are zero. In particular, for (c) note the perspective given in Remark 3, i.e. starting at any point in region R2 will converge to attractors (DAGs) of the form X, + X2 (red line); while starting any point in region R1 will converge attractors (DAGs) of the form X₂ X, (blue line).
Figure 1: Behavior of hider for W
wa w2 0
0
3.2 Why the log-determinant regularizer is preferable to existing acyclicity regularizers
In this section, we present three arguments as to why one should use hijee instead of existing functions such hexpm and hpply- We invite the reader to look at Appendix B for additional details.
Argument (i). hiet does not diminish cycles of any length. Let us expand the functions hexpm and hpoly in their sum of matrix powers form, that is, hexpm (W) Dx 1/ke Tr((W o W)) - d and hpoly (W)- Ed=o (2)/z- Tr(W wY²)
d. Recall also that the entry [(W w)* represents

---

the sum of weighted walks from node i to node i of length k, where each edge has weight w Thus, one can notice that if W has cycles of length their contribution to Aexpm and Apoly are diminished by 1/kl and (2)/a", respectively, Numerically, the latter can be problematic for the following reason: Cycles of length & can go undetected even for small values of d and k. In practice, value of hexpm, Apoly E [10-8 8, 10 is typically regarded as zero [55, 58J. Consider a cycle graph nodes where each edge weight +1 or 1.3 The plot Figure shows how the values of hexpm and pplyy decay much faster than that of bid in fact, at d 13 we already observe hexpm (W) R 10-9 and Ppoly (W) 10-14, i.e., cycles of length at least 13 would be numerically undetected by hexpm and ppoly- In contrast, we observe that the value of the log-det function remains bounded away from zero and is able detect larger cycles.
AP(W)
[W7(W))
180
(W) -..
V (W)
5
Vhhaoy (W7)
k(W)
120
80
60
20
nodes
Figure 2: The values of hexpm and poly get very close to zero for a number of nodes as small as ten. In contrast, letting 8 1.001, we observe that heee can stay away from zero even for a cycle graph of 100 nodes. Finally, we note a similar pattern for the entrywise Eox norm of the gradients.
At first glance, it might seem difficult directly compare the value of hidet to hexpm and hpoly. We next show that when 1, Midet is upper bound to hexpm and hpooy.
Lemma 4. For all W E W==1, we have (W)) hexpm (W) S Lide (W).
The lemma above shows that in spite of ===1 Pidet hesppm and Pplyy being exact acyclicity characteriza-
tions, hi will attain the largest value.
Argument (ii). higer has better behaved gradients. Similar to argument (i), we show in Appendix B that heep and Apply are susceptible to vanishing gradients even when the graph contains cycles (see Figure 2). The following lemma states that the magnitude of each entry of Vhidet least as large as the magnitude of the corresponding entry of V Teeppm and Vh Vhpoly, and hence, hidet has larger gradients to guide optimization.
Lemma 5. For any walk of length k. its contribution to the gradients hexppm (W) and V hh diminished by 1/st and (*z (/a"///-1>", respectively. In contrast, Vhiae (W) does not diminish any walk of any length. This implies that /V (W)l < Whexpp (W)) S Wha (W)l.
Argument (iil). Computing hiatt and Vhidet is empirically faster. Even though hidet hexpm and hpoly all three share the same computational complexity of O(d²), in practice hidet can be computed in about an order of magnitude faster than hexpm and pplyy- In Figure 3, we compare the runtimes of hideer hexpm and liplly for randomly generated matrices, where observe that computing Midet can be 10x faster than hexpm and hppll-- See Appendix B for further details.
Optimization
In the previous section argued why the log-det function should be preferred in practice. Let fo denote a model with parameters 0 for the functions fj in (1), e.g., neural networks as in [59]. In this section we turn to the problem of minimizing given score function Q(fo; X) constrained to Ma (W) = 0. That is, we aim to solve:
min Qffo: X) + |111 subject to (W)))) O,
where we include the & regularizer to promote sparse solutions.
Note that here the sign of edge is not important since W W will have all edge weights equal to +1.
7

---

Runtime ((re)
100.0
10.0
d (number of nodes)
Methods
-
logdot
-
2xpm
poly
Figure For each d, 30 matrices were randomly sampled from standard Gaussian distribution.
Since the inception of the purely continuous framework for learning DAGs in [58], almost all follow- up work uses the augmented Lagrangian (AML) scheme to tackle problem (6), and L-BFGS-B [39] for solving the sequence of unconstrained problems. Motivated by the properties of hidet given in
3.1, we propose a simpler scheme named DAGMA, based on solving a sequence of unconstrained problems in which Maet is simply
regularizer. DAGMA resembles the central path approach of barrier methods [8, 39], or the classical path-following approach for solving lasso problems [e.g. 14]. Our method is given in Algorithm 1.
Algorithm 1 DAGMA
Require: Data matrix X, initial central path coefficient (9) (e.g., 1), decay factor a E (0, (e.g.. 0.1), & parameter 31 > 0 (e.g.. 0.01), log-det parameter s >> 0 (e.g.. 1), number of iterations T. 1: Initialize ((0) so that (() W
2: for 0, 1, 2,
do
3: Starting at ((1), solve g(++1) = arg ming wiinn,,' (Q X) (111011) + hidet (W)))
4: Set (1(11) = (()(
Ensure: W(B(T)
The following lemma states that DAGMA will return DAG at the limit of the central path. This is critical distinction against existing methods, many of which rely on some type of post-processing (e.g. thresholding) to ensure that the solution is DAG.
Lemma 6. Algorithm is guaranteed to return DAG whenever pl() - O.
4.1 Practical Considerations
1. As in barrier methods, where it is required to start at the interior of the feasibility region, in Algorithm 1, require that the initial point W W((()) be inside We This is very easy to achieve since the zero matrix is the interior of W for any s > 0; therefore, in our
experiments we simply set ((0) =
2. Note that in Algorithm 1, we let (0) decrease by a constant factor at each iteration; however, it is possible to specify explicitly the value of each (((),,, c.g., for = 4, we can let / {1, 0.1, 0..
3. Regarding the choice of in principle could take any value greater than zero since DAGs are inside 772 for any 0. Similar to μı, is also possible to let s vary at each iteration, e.g., for T 4, can set decreasing can help to obtain larger gradients as gets closer to DAG. Note, however, that letting be equal or close to is generally easier to optimize than setting closer to zero, the reason being that for smaller values of s the volume of W is smaller and will require much smaller learning rates to stay inside W hence affecting convergence.
{1, 0.9, 0.8, 0.8}. In practive, we observe that slightly
4. Finally, we do not specify how to solve line 3 in Algorithm 1, this is because we leave the door open for different solvers to be used. For our experiments in the next section, we solve line 3 by using first-order method with the ADAM optimizer [24], which works remarkably well as shown in our experiments. remains as future work to exploit the Hessian structure of fudet given in Lemma 3 for second-orde methods.

---

5 Experiments
We compare our method against GES [11], PC [SI], NOTEARS [58], and GOLEM [37] on both linear and nonlitear SEMs. In Appendix C, we specify which existing implementation used for each of the aforementioned methods. Consistent with previous work in this area (e.g., NOTEARS and follow up work), we have not performed any hyperparameter optimization: This is to avoid presenting unintentionally biased results. As a concrete example, for each of the SEM settings, we simply chose reasonable value for the l penalty coefficient and used that same value for all graphs across many different numbers of nodes.
Our experimental setting is similar to [58, 59]. For the main text, we present only small fraction of all our experiments. Moreover, since the accuracy of certain methods were significantly lower than other methods, report results only against the most competitive ones; full results for all settings and methods can be found in Appendix C.
Linear Models. In Appendix C.I, we report results for linear SEMs with Gaussian, Gumbel, and exponential noises, and use the least squares loss. For small to moderate number of nodes, see Appendix C.I.1; for large number of nodes, see Appendix C. for denser graphs, see Appendix C.I.4; and for a comparison against GOLEM for sparser graphs, see Appendix C.I.3.
Nonlinear Models. In Appendix C.2, we report results for nonlinear SEMs with binary and con- tinuous data. For binary data, we use a logistic model for each structural equation, and use the log-likelihood loss as the score, we report results for small to large number of nodes in Appendix C.2.1. For continuous data, we consider the continuous additive noise model with Gaussian noise [9], where each nonlinear relationship is modeled by a multilayer perceptron, and use the oog-iikeidood loss as the score, we report results for small to moderate number of nodes in Appendix C.2.2.
In the following figures, ER4 and SF4 denote Erdös-Rényí and scale-free graphs, respectively, where for each number of nodes d, each graph has in expectation 4d edges. It is worth noting that the empirical settings by Zheng et al. [58, 59] consider graph models such as ERI, ER2. SF1, and SF2. Here we focus on the hardest setting, i.c., ER4 and SF4 graphs. For linear SEMs, Figure 4 shows results for graphs with d E [20, 100), and Figure 5 shows results for graphs with d E [200, 1000]. In both regimes, we note that DAGMA obtains significant speedups and improvements terms of structural accuracy (SHD) against NOTEARS and GOLEM, even though GOLEM is specific to and specialized for linear Gaussian SEMs. For nonlinear SEMs, Figure shows results for logistic models, we similarly observe that DAGMA attains major speedups and improvements on SHD against NOTEARS [58]. Finally, for nonlinear models using neural networks, we observe that DAGMA is comparable in SHD to the NONLINEAR NOTEARS [59] but obtains significant speedups. Again, we invite the reader to look at Appendix C for more details and additional experiments.
Final Remarks
relevant assumption in this work is that of sufficiency. that is, there are no hidden variables that are a common cause of at least two observed variables. While this assumption is widely used for structure learning, we nonetheless highlight that in practice it is very difficult to find scenarios where such assumption holds. As with all work that assumes sufficiency, our work is an important necessary step to understanding settings with hidden variables. Finally, we note that the work by [6] proposes differentiable approach for ADMG for the semi- Mari-Maroovian case using hexpm It is left for future work to explore the performance of such method using hidet-
Another important limitation of this and previous work on the continuous framework for learning DAGs is that of providing guarantees on the learned structure. As in real-life applications one does not have access to the mouund-rruth DAG, there is much uncertainty as to whether an edge in the predicted DAG actually corresponds to a causal relation. Thus, there is still a need for formal guarantees under the continuous framework.
A Acknowledgments and Disclosure of Funding
K. B. was supported by NSF under Grant # 2127309 to the Computing Research Association for the CIFellows 2021 Project. B.A. was supported by NSF 115-1956330, NIH
and the
Robert Topel Faculty Research Fund at the University of Chicago Booth School of Business. P.R. was supported by ONR via N000141812861, and NSF via IIS- 190989, IIS- 1955532, IIS 11922119077

---

100
500-
140
300
170
300
100
100
30-
10-
040
(seconds))
220
300
100
100
30
Structural
10
140
1000
500
100
100
300
100
50-
10
20
60
00
10
105
50
60
80
80
100
20
40
05
80
100
70 05
80
100
[number nodes)
[number of nades)
Methods
Figure 4: Experiments on linear SEMs for E [20, 100). Each point in the plot is estimated over 10 repetitions, where error bars the standard error. Wall time limit to 36 hours.
5000
3020
30.00
2000
10.00
1000
380
500
350
1.00
200
150
8.30
30
Structunl temmiing Dstanne
50
10
0.10
3838
2003
1000
1000
368
Runtrr (hourgg)
3.00
100-00
8.30
703
10
0.10-
3050
2000
1000
10.00
380
500
1.00
300
38
8.0
200
0.05
200
400
630
800
1200
200 400
600
802
1000
230
400
000
300
1500
250
400
650
800
1000
(number of nodes)
(number of nodes)
Figure 5: Experiments on linear SEMs for d E [200, 1000]. Each point in the plot is estimated over 10 repetitions, where error bars the standard error. Wall time limit was set to 36 hours.
Heaa Salarra
frandoro (ods)
d of of
Figure 6: Experiments on nonlinear SEMs. Each point in the plot is estimated over 10 repetitions. where eтrог bars are the standard error. Wall time limit was set to 36 hours.
10

---

References
[1] AI-Mohy, A. H. and Higham, N. J. [2011], 'Computing the Action of the Matrix Exponen- tial, with an Application to Exponential Integrators', SIAM Journal on Scientific Computing 33(2), 488-511.
[2] Aragam, B., Amini, A. A. and Zhou, Q- [2019], 'Globally optimal score-based learning of directed acyelic graphs in high-dimensions Advances in Neural Information Processing Systems.
[3] Aragam, B. and Zhou, Q. [2015], 'Concave Penalized Estimation of Sparse GaussianBayesian Networks', Journal of Machine Learning Research.
[4] Barabási, A.-L. and Albert, R. [1999]. 'Emergence of scaling in random networks', science 286(5439), 509-512.
[5] Berman, A. and Plemmons, R. [1994], Nonnegative Matrices in the Mathematical Sciences, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics.
[6] Bhattacharya, R., Nagarajan, T., Malinsky, D. and Shpitser, I. [2020], 'Differentiable Causal Discovery Under Unmeasured Confounding', arXiv
[7] Bouckaert, R. R. [1993]. Probabilistic network construction using the minimum description length principle, in 'European conference on symbolic and quantitative approaches to reasoning and uncertainty", Springer, pp. 41-48.
[8] Boyd, S. and Vandenberghe, L. [2004], Convex Optimization, Cambridge University Press.
[9] Bühlmann, P., Peters, J. and Ernest, J. [2014], 'Cam: Causal additive models, high-dimensional order search and penalized regression', The Annals of Statistics 42(6), 2526-2556.
[10] Chickering, D. M. [1996], Learning bayesian networks is up-complete, 'Learning from data',
Springer, pp. 121-130.
[11] Chickering, D. M. [2003], 'Optimal structure identification with greedy search', JMLR
3, 507-554.
[12] Chickering, D. M., Heckerman, D. and Meek, C. [2004], 'Large-sample learning of Bayesian networks is NP-hard', Journal of Machine Learning Research 5, 1287-1330.
[13] Ernest, J., Rothenhäusler, D. and Bühlmann, P. [2016]. 'Causal inference in partially linear structural equation models: identifiability and estimation', arXiv preprint 1607. 0507005980..
[14] Friedman, J., Hastie, T., Höffing, H. and Tibshirani, R. [2007), 'Pathwise coordinate optimiza- tion', The annals of applied statistics 1(2), 302-332.
[15] Gámez, A., Mateo, L.. and Puerta, J. M. [2011], 'Learning bayesian networks by hill climbing: efficient methods based on progressive restriction of the neighborhood', Data Mining and Knowledge Discovery 22(1), 106-148.
[16] Ghoshal, A. and Honorio, J. [2017], Learning identifiable gaussian bayesian networks in polynomial time and sample complexity, in 'Proceedings of the 31 International Conference on Neural Information Processing Systems', pp. 6460-6469
[17] Ghoshal, A. and Honorio. J. [2018], Learning linear structural equation models in polynomial time and sample complexity, 'Proceedings of the Twenty-Firs International Conference on Artificial Intelligence and Statistics', Vol. 84 of Proceedings of Machine Learning Research, PMLR, pp. 1466-143 1475.
[18) Goudet, O., Kalainathan, D., Caillou, P., Guyon, Lopez-Paz, D. and Sebag, M. [2018), Learn ing functional causal models with generative neural networks, in 'Explainable and interpretable models in computer vision and machine learning'. Springer, pp. 39-80.
[19] Gu, J., Fu, F. and Zhou, Q [2019], 'Penalized estimation of directed acyclic graphs from discrete data', Statistics and Computing 29(1), 161-176.
[20] Hanson, M. A. [1981], 'On sufficiency of the Kuhn-Tucker conditions' canditions",Juunnne ,Journal of Mathematical
Analysis and Applications 80(2), 545-550.
[21] Heckerman, D... Geiger, D. and Chickering, D. M. [1995]. 'Learning bayesian networks: The combination of knowledge and statistical data', Machine learning 20(3), 197-243.
[22] Higham, N. [2005], 'The scaling and squaring method for the matrix exponential revisited',
SIAM Journal on Matrix Analysis and Applications 26(4), 1179-1193.
11

---

[23] Hoyer. P., Janzing, D., Mooij, J. M., Peters, J. and Schölkopf, B. [2008], 'Nonlinear causal discovery with additive noise models', Advances in neural information processing systems 21.
[24] Kingma, D. P. and Ba, J. [2014], 'Adam: A method for stochastic optimization', arXiv preprint arXiv:1412.6980
[25] Kyono, T., Zhang, Y. and van der Schaar, M. [2020], 'Castle: Regularization via auxiliary causal
graph discovery", Advances in Neural Information Processing Systems 33, 1501-1512.
[26] Lachapelle, S., Brouillard, P., Deleu, T. and Lacoste- Jullenn S. [2020], Gradient-basee neural dag learning, in 'nternationaa Conference on Learning Representations'
[27] Loh, P.-L. and Buhlmann, P. [2014), Diigh-Dimensional Learning of Linear Causal Networks via Inverse Covariance Estimation' Journal of Machine Learning Research.
(28) Magnus, J. R. and Neudecker, H. [1985], 'Matrix differential calculus with applications to simple, hadamard, and kronecker products', Journal of Mathematical Psychology 29(4). 474- 492.
[29] Margaritis, D. and Thrun, S. [1999], Bayesian network induction via local neighborhoods, in 'Proceedings of the 12th International Conference on Neural Information Processing Systems', pp. 505-511.
[30] Martin, D. H. [1985], 'The essence of invexity', Journal of Optimization Theory and Applica- tions 47(1), 65-76.
[31] Maxwell Chickering, D. and Heckerman, D. [1997], 'Effficient approximations for the marginal likelihood of bayesian networks with hidden variables', Machine learning 29(2), 181-212.
[32] Meinshausen, N. and Bühlmann, P. [2006], niigh-dimensional graphs and variable selection with the Lasso', The Annals of Statistics 34(3).
[33] Moler, C. and Van Loan, C. [2003], 'Nineteen dubious ways to compute the exponential of matrix, twenty-five years later', SIAM review 45(1), 3-49.
[34] Monti, P., Zhang, K. and Hyvärinen, A. [2020], Causal discovery with general non-linear relationships using non-linear ica, in 'Uncertainty in Artificial Intelligence', PMLR, pp. 186- 195.
[35] Mooij. J. M., Peters, J., Janzing, D., Zscheischler, and Schölkopf, B. [2016], 'Distinguishing cause from effect using observational data: methods and benchmarks', The Journal of Machine Learning Research 17(1), 1104-1204.
[36] Moraffah, R., Moraffah, B., Karami, M., Raglin, A. and Liu, H. [2020], 'Causal adversarial network for learning conditional and interventional distributions', arXiv:2008. 10:::2008.1136....
[37] Ng. I., Ghassami, A. and Zhang, K. [2020], 'On the Role of Sparsity and DAG Constraints for Learning Linear DAGs', arXiv.
[38] Ng. I., Lachapelle, S., Ke, N. R. and Lacoste- Julien, S. [2022], "On the Convergence of Continuous Constrained Optimization for Structure Learning', AISTATS.
[39] Nocedal, J. and Wright, S. J. [2006), Numerical Optimization, Springer Series in Operations Research and Financial Engineering.
[40] Ostrowski, A. [1937], *Über die determinanten mit überwiegender Hauptdiagonale', Commen- tarii Mahematici Helvetici 10(1), 69-96.
[41] Pamfil, R., Sriwattanaworachaii, N., Desai, S., Pilgerstorfer, P., Georgatzis, K., Beaumont, P. and Aragam, B. [2020], Dynotears: Structure learning from time-series data, in 'International Conference on Artificial Intelligence and Statistics', PMLR, pp. 1595-1605.
[42] Park, G. and Park, S. [2019], 'High-dimensional poisson structural equation model learning via ell I-regularized regression'', Mach. Learn. Res. 20, 95-1.
[43] Park, G. and Raskutti, G. [2017], 'Learning quadratic variance function (qvf) dag models via overdispersion scoring (ods)., J. Mach. Learn. Res. 18, 224-1.
[44] Pearl, J. [2009], CAUSALITY: Models, Reasoning, and Inference, 2nd edn, Cambridge University Press.
[45] Peters, J. and Bühlmann, P. [2014], 'Identifiability of gaussian structural equation models with equal error variances', Biometrika 101(1), 219-228.
12

---

[46] Peters, J.. Mooij. J. M., Janzing, D. and Bernhard [2014], 'Causal Discovery with Continuous Additive Noise Models', Journal of Machine Learning Research
[47]
Petersen, K. B. and Pedersen, M. S. [2012], The Matrix Cookbook.
[48] Ramsey, J., Glymour, M., Samchez-Romero, R. and Glymour, C. [2017], 'A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application functional magnetic resonance images', International journal of data science and analytics 3(2), 121-129.
[49] Shimizu, S., Hoyer, P. O.. Hyvärinen, A., Kerminen, A. and Jordan, M. [2006], 'A linear non non-gaussian acyclic model for causal discovery'' Journal of Machine Learning Research 7(10).
[50] Spirtes, P. and Glymour, C. [1991], 'An algorithm for fast recovery of sparse causal graphs', Social Science Computer Review 9(), 62-72.
[51] Spirtes, P., Glymour, C. and Scheines, R. [2000], Causation, Prediction, and Search, The MIT Press.
[52] Tsamardinos, I., Aliferis, C. F., Statnikov, A. R. and Statnikov, E. [2003], Algorithms for large scale markov blanket discovery, 'FLAIRS conference", Vol. pp. 376-380.
[53] Tsamardinos, I., Brown, L.. E. and Aliferis, C. F. [2006], 'The max-min hill-climbing Bayesian network structure learning algorithm', Machine Learning 65(1), 31-78.
[54] Voorman, A., Shojaie, A. and Witten, D. [2014], 'Graph estimation with joint additive models", Biometrika 101(1), 85-101.
[55] Wei, D., Gao, T. and Yu, Y. [2020], "DAGs with No Fears: A Closer Look at Continuous Optimization for Learning Bayesian Networks', arXiv.
[56] Yu, Y., Chen, Gao, T. and Yu, M. [2019], 'DAG-GNN: DAG Structure Learning with Graph Neural Networks', arXiv.
[57] Yu, Y., Gao, T., Yin, N. and Ji, Q. [2021], 'DAGs with No Curl: An Efficient DAG Structure Learning Approach', arXiv.
[58] Zheng, X., Aragam, B.. Ravikumar, P. and Xing, E. P. [2018]. 'DAGs with NO TEARS: Continuous Optimization for Structure Learning', arXiv.
[59] Zheng, X., Dan, C., Aragam, B., Ravikumar, P. and Xing, E. [2020], 'Learning Sparse
Nonparametric DAGs', International Conference on Artificial Intelligence and Statistics.
[60] Zhu, S., Ng, I. and Chen, Z. [2020], Causal discovery with reinforcement learning, in 'Interna- tional Conference Learning Representations'
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] (See Section 6)
(c) Did you discuss any potential negative societal impacts of your work? [Yes] (See Section D)
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] (See Section and 3)
(b) Did you include complete proofs of all theoretical results? [Yes] (See Section A)
3. If you experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL.? [Yes] (See Section
C)
13

---

(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] (See Section C)
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] (See Section C)
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] (See Section C)
4. If you are using existing assets (c.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes] (See Section C)
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "lggm4rIfvr", "reviewer": "Reviewer_UUUC", "summary": "The paper introduces a new penalty for enforcing acyclicity of a weighted adjacency matrix for learning DAGs with differentiable optimization. The penalty is zero if and only if the adjacency matrix is acyclic. In contrast to previous proposals based on matrix polynomials, the introduced penalty is based on the log-determinant and is connected to the theory of M-matrices. The paper studies the properties of the new penalty and compares the penalty with previous proposals. From theoretical analysis and empirical observations, it is argued that the new penalty is more suited for learning DAGs / SEMs. Through numerical experiments, it is shown that the new penalty leads to speedups and higher precisions compared to previous proposals. ", "questions": "1. On discounting long cycles.\n\nArgument (i) claims that the new penalty does not discount a long cycle in the same way some previous proposals do. From Figure 2, it looks persuasive that this could be a big advantage. I would like to see some asymptotic analysis that makes this comparison precise. \n\n2. On vanishing gradients. \n\nLemma 5 seems insufficient to show that the new penalty does not suffer from vanishing gradients (when the other penalties do). Is there a more compelling argument?\n\n3. On computational speedups. \n\nDespite the same computational complexity, argument (iii) curiously claims that the new penalty \"can be computed in about an order of magnitude faster\" than the other two in practice. Is there an explanation?\n\n4. Minor issues\n\n(1) line 42: missing power in the expression?\n\n(2) Lemma 1: I would rather state (i) as \"DAGs are in the interior of $\\mathbb{W}^s$\".", "limitations": "I do not foresee any potential negative societal impact of their work", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n1. The paper makes progress on recent developments for differentiable learning of DAGs. \n2. The proposed penalty seems novel and has interesting connections to M-matrices. \n3. The paper compares the new penalty with previous ones in the literature from various aspects. \n4. The paper is very well-written. \n\n### Weaknesses\n1. Several aspects in the comparison are based on heuristic arguments or empirical observations, which weakens the claim that the new penalty should be preferred over the earlier ones. A theoretical analysis is desirable to show that the advantage is real or \"significant\" in typical scenarios or when the problem size grows large. \n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "T57WmLBTF6B", "reviewer": "Reviewer_kNuE", "summary": "The authors proposed a new DAG constraints based on the log-determinant of the adjacency matrices. Compare to the original polynomial based DAG constraints, the new DAG constraints do have some good properties and it achieves better performance on ER4 and SF4 graphs with node from 200 to 1000.", "questions": "Is the performance from different implementation? In NOTEARS, the authors does not decrease mu but increase the coefficients of h(W dot W). In GOLEM, the authors use fixed coefficient for L1. If put the DAG constraint in the same optimisation framework as NOTEARS, or if we put the NOTEARS in the same framework as the work, what is the performance?", "limitations": "It is discussed.", "rating": 5, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Pros:\n\n1. The proof and derivation of the new DAG constraint is very clear.\n2. The illustration of the behaviour of h_det is very good.\n\n\nCons:\n\n1. The properties (ii) in Theorem 1 is actually as bad as the Exponential or Polynomial based DAG constraint. See [1,2] for more details. With  this property you will not be able to obtain an DAG unless you let mu in Algorithm 1 converge to zero. \n\n2. Corollary 2 may not be correct and Corollary 1 maybe meaningless. For Corollary 2, it only considers h_det, without consider the score function. By [1,2] without attain the limit point it would not be possible to obtain a DAG, and this will require infinite many steps of optimisation.\n\n3. In the experiment part, the hyper parameter for GOLEM may not be quite right. I have tried GOLEM on SF{1,2,3,4}, and ER{1,2,3}, for numbers of nodes in {10, 20, 50, 100}. The performance of GOLEM performance is consistently better than NOTEARS in a large margin, which is quite different from the result reported in the paper. Also for all cases it would be good if the authors can provide results on graphs with nodes in  {10, 20, 50, 100}. This is because with the same number of expected edges, the more the nodes, the sparser the graph is. Also it would be good if the authors can provide results on ER2, ER3 and SF2, SF3, SF4 with all number of nodes from 10 to 1000. This can provide a fair comparison for different methods in different situations. I have also attached the GOLEM code I have used, you only need to replace the HTorch part with the exponential DAG constraints from NOTEARS.\n\n\n\n[1] Wei, Dennis, Tian Gao, and Yue Yu. \"DAGs with No Fears: A closer look at continuous optimization for learning Bayesian networks.\" Advances in Neural Information Processing Systems 33 (2020): 3895-3906.\n[2] Ng, Ignavier, et al. \"On the convergence of continuous constrained optimization for structure learning.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2022.\n\n\nThe GOLEM code I have used.\n\n\n\"\"\"\n\n    import numpy as np\n    import torch\n    import torch.nn as nn\n    from .dag import HTorch\n    h_torch = HTorch.apply\n    class GolemEVModel(nn.Module):\n        \"\"\"\n        Set up the objective function of GOLEM.\n        Hyperparameters\n        (1) GOLEM-NV: equal_variances=False, lambda_1=2e-3, lambda_2=5.0.\n        (2) GOLEM-EV: equal_variances=True, lambda_1=2e-2, lambda_2=5.0.\n        \"\"\"\n        def __init__(self, d, lambda_1, lambda_2, eps=1e-6, h_type='exponential'):\n            \"\"\"\n            Initialize self.\n            Parameters\n            d: int\n                Number of nodes.\n            lambda_1: float\n                Coefficient of L1 penalty.\n            lambda_2: float\n                Coefficient of DAG penalty.\n            equal_variances: bool\n                Whether to assume equal noise variances\n                for likelibood objective. Default: True.\n            \"\"\"\n            super().__init__()\n\n            self.d = d\n            self.lambda_1 = lambda_1\n            self.lambda_2 = lambda_2\n            self.eps = eps\n            self.h_type = h_type\n            self._B = nn.Parameter(torch.zeros(self.d, self.d))\n\n        def forward(self, cov_emp):\n            # Placeholders and variables\n            self.cov_emp = cov_emp\n            self.B = self._preprocess(self._B)\n\n            # Likelihood, penalty terms and score\n            self.likelihood = self._compute_likelihood()\n            self.L1_penalty = self._compute_L1_penalty()\n            self.h = self._compute_h()\n            self.score = (self.likelihood + self.lambda_1 * self.L1_penalty +\n                      self.lambda_2 * self.h)\n\n        def _preprocess(self, B):\n            \"\"\"\n            Set the diagonals of B to zero.\n            Parameters\n            B: tf.Tensor\n                [d, d] weighted matrix.\n\n            Return\n            torch.Tensor: [d, d] weighted matrix.\n            \"\"\"\n            return (1. - torch.eye(self.d)) * B\n    \n        def _compute_likelihood(self):\n            \"\"\"\n            Compute (negative log) likelihood in the linear Gaussian case.\n            Return\n            torch.Tensor: Likelihood term (scalar-valued).\n            \"\"\"\n            I = torch.eye(self.d)\n            return 0.5 * self.d * torch.log(\n                torch.trace((I - self.B).T @ self.cov_emp\n                        @ (I - self.B))) - torch.linalg.slogdet(I - self.B)[1]\n\n        def _compute_L1_penalty(self):\n            \"\"\"\n            Compute L1 penalty.\n            Return\n            tf.Tensor: L1 penalty term (scalar-valued).\n            \"\"\"\n            return torch.norm(self.B, p=1)\n\n        def _compute_h(self):\n            \"\"\"\n            Compute DAG penalty.\n\n            Return\n            torch.Tensor: DAG penalty term (scalar-valued).\n            \"\"\"\n            return h_torch(self.B * self.B, self.h_type, self.eps)\n    def golem_ev(X,\n             lambda_1=2e-2,\n             lambda_2=5.0,\n             learning_rate=1e-3,\n             num_iter=2e+4,\n             graph_thres=0.3,\n             eps=1e-6,\n             h_type='exponential'):\n        n, d = X.shape\n        cov_emp = np.cov(X.T, bias=True)\n        cov_emp = torch.Tensor(cov_emp)\n        model = GolemEVModel(d, lambda_1, lambda_2, eps, h_type)\n        train_op = torch.optim.Adam(model.parameters(), lr=learning_rate)\n\n        for i in range(int(num_iter)):\n            model(cov_emp)\n            score, likelihood, h, B_est = model.score, model.likelihood, model.h, model.B\n            loss = score\n            train_op.zero_grad()\n            loss.backward()\n            train_op.step()\n   \n        return B_est.detach().numpy()\n\"\"\"", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "PbIfP2QdwUS", "reviewer": "Reviewer_vTCA", "summary": "* This paper considers the problem of learning a structure (DAG) over a set of observed variables with gradient-based methods. The approach designs a new continuous optimization problem capturing the acyclic property of DAGs, which defines the domain of the log-det function to be M-matrices. \n* The authors provide a detailed technical description of the approach as well as extensive experiments. \n* The authors describe technically that the proposed approach has more well behaved gradients (perhaps indicative of the improved performance empirically) as well as several other theoretical properties of the approach. \n* Empirically, the authors provide extensive experiments which demonstrate the efficiency and effectiveness of the proposed approach. ", "questions": "* Rather than sampling a DAG structure by ordering the nodes in an undirected graph, why not sample DAGs from something like Price's model? I realize that your experiment setup is following earlier work however. I am curious what if any properties we should be aware of for your generative model of dags. How do in-degree / out-degree distributions look like? How do they differ from degree distribution from undirected graph?\n* When / why would quality $Q(\\cdot)$ make sense / not-make sense for an evaluation measure of the DAGs? Could it be evaluated on held-out data of the same underlying distribution? \n* When / would it make sense to evaluate not just wall clock time, but number of samples $n$? \n* It seems like hyperparameter tuning could play a large role in the performance of these approaches. Do I understand correctly that previous work has tuned hyperparameters for these same exact datasets in the cases of re-use? Is there clarity on how all hyperparemters were selected? I wonder what you think about evaluating methods using bayesian hyperparameter optimization to find, for each dataset (e.g., setting of generative model hyperparameters) best setting of hyperparameters. \n* There is a note about the proposed approach being better at detecting large cycles. For give me if I have missed something, but was there analysis in the experiments section, which noted how much of the improved performance (e.g. SHD) was due to better recovery of such large cycles? \n* If these methods were used as part of some downstream approach, do you think the empirical analysis here is predictive of performance? Why/why not / what other measurements would be interesting to consider?", "limitations": "Yes, the authors have done a very nice job with this.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall I believe that this is a strong paper with an interesting approach for what is becoming a more well studied area (continuous optimization approaches for structure learning). \n\nIn particular strengths include: \n\n* **Technical Contribution** - I believe the authors have clearly explained the problem they are solving, why their methodical approach solves the problem, and what are the challenges and insights that are required in the contribution. I find the contribution to be interesting and meaningful. \n* **Empirical benefits** - The authors have demonstrated the performance characteristics of their approach with extensive experiments, which would allow a practitioner to consider such results and then determine which methods are most suitable for the application (based on metrics such as time, SHD, TPR, etc). And most often, DAGMA would be a top selection, especially in terms of efficiency compared to methods like NOTEARS. \n\nI have a few questions which indicate places for improvement regarding setup of experiments and measurements, please refer to questions section for detailed description of these weaknesses. \n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "2eEjhG7YFQ", "reviewer": "Reviewer_wwF9", "summary": "This paper proposes a new formulation promoting acyclicity. Going beyond the widely used approaches based on the trace of polynomials of the hadamard square of the adjacency matrix, this paper rather takes a different angle by considering a reduced set of matrices, M-matrices. This M-matrix based acyclicity regularizer is thoroughly investigated so that its empirically superior behavior is also understood theoretically.", "questions": "* In Corollary 1, what does it mean by 'point towards the interior'. Is this term 'interior' is the topological term or a general term for somewhere inside? It seems that $\\mathbb{W}^s$ is an open subset of $R^{d \\times d}$. Does this mean that its direction is always directly toward a stationary point?\n\n* It seems that in eq(5) the condition that $W \\ge 0$, that is, all entries are nonnegative is missing. It seems that some Theorem may hold without this nonnegative entry condition but others do not. Theorems/Lemma critically relying on Proposition 1 should require the nonnegative entry condition. However, for example, Lemma 2 seems to consider the cases without nonnegative entry conditions, otherwise, $sign(W_{i,j})$ is always nonnegative. It would be better if the authors clarify this.\n\n", "limitations": "* Except for general limitations of differentiable approaches for DAG discovery, there is no other specific limitation. To discuss such such general limitations some other references can be added. In Differentiable Causal Discovery Under Unmeasured Confounding; Rohit Bhattacharya, Tushar Nagarajan, Daniel Malinsky, Ilya Shpitser; AISTATS 2021. there are many interesting discussions and some solutions on such general limitations.\n  * Differentiable approach for ADMG for semi-Markovian case\n  * How to find MEC not a single DAG", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n* A novel approach based on M-matrix is proposed for an acyclicity regularizer free from optimization hassle in exiting approaches.\n* The behavior of the log-det acyclicity regularizer -- the gradient and other regularization properties -- is analyzed to a quite detailed extent with supporting simulation results.\n* DAG discovery with the proposed method outperforms others even with the advantage in the run-time.\n\nWeaknesses\n* For some theoretical results, it would be better if some technical detail is provided at a level to help to understand. \n  * In line 177, $(4) \\Rightarrow (3)$ required a bit of thought to understand so some brief comments would help instantaneous understanding.\n  * In Lemma 3, the statement that $\\otimes$ is the Kronecker product may be added so that reader can look up its definition.\n* It seems that sometimes theorems assumes slightly different definition for $\\mathbb{W}^2$ with/without nonnegative entry condition (details below)", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept