Constraint-Free Structure Learning with Smooth Acyclic Orientations
Authors: Riccardo Massidda
Francesco Landolfi
Martina Cinquini
Davide Bacciu
Venue: ICLR 2024 poster (2024)
Keywords: Structure Learning
Continuous Acyclic Optimization
Directed Acyclic Graphs
DAG
Primary area: general machine learning (i.e., none of the above)

Abstract: The structure learning problem consists of fitting data generated by a Directed Acyclic Graph (DAG) to correctly reconstruct its arcs. In this context, differentiable approaches constrain or regularize an optimization problem with a continuous relaxation of the acyclicity property. The computational cost of evaluating graph acyclicity is cubic on the number of nodes and significantly affects scalability. In this paper, we introduce COSMO, a constraint-free continuous optimization scheme for acyclic structure learning. At the core of our method lies a novel differentiable approximation of an orientation matrix parameterized by a single priority vector. Differently from previous works, our parameterization fits a smooth orientation matrix and the resulting acyclic adjacency matrix without evaluating acyclicity at any step. Despite this absence, we prove that COSMO always converges to an acyclic solution. In addition to being asymptotically faster, our empirical analysis highlights how COSMO performance on graph reconstruction compares favorably with competing structure learning methods.

Published as a conference paper at ICLR 2024
CONSTRAINT-FREB STRUCTURE LEARNING WITH SMOOTH ACYCLIC ORIENTATIONS
Riccardo Massidda, Francesco Landolfi, Martina Cinquini, Davide Bacciu Department of Computer Science
Università di Pisa, Italy
{riccardo mascardoomasslooooooooooo mass dda, fraaaaaannn francesco. landol fi, martmma..ainnaang...... in [ iit davide www.ciuuuuiuiii......
ABSTRACT
The structure learning problem consists of fitting data generated by Directed Acyclic Graph (DAG) to correctly reconstruct ares. In this context, differentiable approaches constrain relaxation of the acyclicity property. The computational cost of evaluating graph acyclicity is cubic the number of nodes and significantly affects scalability. In this paper, we introduce COSMO, a constraint-fres continuous optimization scheme for acyclic structure learning. At the core of our method lies novel differentiable approximation of an orientation matrix parameterized by a single priority vector. Differently from previous works, our parameterization fits smooth orientation matrix and the resulting acyelic adjacency matrix without evaluating acyclicity at any step. Despite this absence. we prove that COSMO always converges an acyclic solution. In addition to being asymptotically faster, our empirical analysis highlights how COSMO performance on graph reconstruction compares favorably with competing structure learning methods.
regularize optimization problem with continuous
INTRODUCTION
Directed Acyclic Graphs (DAGs) are a fundamental tool in several fields to represent probabilistic or causal information about the world (Koller& & Friedman, 22099] Pearl 2009). A fundamental problem in this context concerns the retrieval of the underlying structure between a set of variables, i.c., the problem of identifying which ares exist between nodes associated to the variables of interest (Spirtes) et al. (0000). In recent years, applications of structure learning to causal discovery led to growing interest tackling the problem using gradient- baseed methods that optimize a smooth representation of DAG (Vowels et al. 2022). For instance, while not suitable for causal discovery per (Reisach et al.,
202T), acyclic structure learners are fundamental components of most state-of-the-art continuous causal discovery algorithms (Lachapelle et al., 2020 Brouillard et al. 2020 Lorch et al., (002).. A
well-established technique, popularized by NOTEARS (Zheng et al. 2018), consists of computing
the trace of the matrix-exponential of the adjacency matrix, which differentiable and provably zero if and only if the corresponding graph is acyclic. However, despite their widespread adoption, NOTEARS-liK acyclicity constraints impose cubic number of operations in the number of nodes per optimization step and substantially prevent scalable and applicable continuous discovery algorithms.
In this context, we propose novel formulation and optimization scheme for learning acyclic graphs that avoids evaluating the acyclicity of the solution in any optimization step. Notably, our proposal does not sacrifice theoretical guarantees of asymptotic convergence acyclic solutions which apply to existing structure learning methods (Ng et al. 2022a). At the core of our scheme lies novel definition of smooth orientation matrix, i.e., differentiable approximation of an orientation matrix parameterized by priority vector on the graph nodes. The priority vector represents discrete orientation where each node has an outgoing are to all nodes with higher priority. We define our smooth orientation matrix by applying tempered sigmoid to the pair-wise priority differences, which equals the discrete orientation in the limit of the sigmoid temperature to zero. By annealing effectively decreasing upper bound the
temperature during training, we prove that acyelicity of the solution. Further, we show that the parameterization represents the space of DAGs as differentiable function of directed graph and our smooth orientation matrix. Since approach

---

Published as a conference paper at ICLR 2024
[4.2 4.1 2.0 -0.4 1.9)
Ste
X1
X2
X
X2
X1
X2
X
X3
X
X4
X
X
Xs
X4
X
Directed Graph H
Smooth Acyclic Orientation
Directed Acyclic Graph
In COSMO propose to optimize acyclic adjacency matrix W E Raxd (right) by learning a directed graph H Rdd (leff) and a priority vector p E R on the nodes (top). To this end, introduce smooth acyclic orientation matrix (P) (center), differentiable function of the priority vector where each lower- tever-priority node approximately feeds into each higher higher-priority node. Gray dashed arrows denote arcs with approximately training, the smooth orientation matrix St.- converges to a discrete orientation and, consequently, the overall adjacency matrix W converges DAG.
Figure 1:
weight. By annealing its temperature during
only requires quadratic number of operations per optimization step, its constraint-fres scheme can be used as direct and faster replacement for the NOTEARS constrained optimization problem. Overall, we propose methodology perform constraint-free structure learning with smooth acyelic orientations, which name COSMO (Figure
Contributions. We summarize the key contributions of this paper as follows:
We introduce a differentiable relaxation of an acyclic orientation matrix, which we name smooth orientation matrix (Definition [T. The matrix depends on temperature value that controls the approximation of the discrete orientation matrix. We prove that we can represent all and only DAGs as the element-wise multiplication of weighted adjacency matrix and our novel smooth orientation matrix (Theorem)
We propose COSMO, an unconstrained optimization approach that learns a DAG entirely avoiding acyclicity constraints (Section 4.2). COSMO represents DAGs through smooth orientation matrix and requires solving unique optimization problem while annealing the temperature. Since reconstructing the DAG requires a number of operations quadratic on the number of nodes, COSMO is an order of magnitude faster than methods in literature.
constrained
We connect our proposed scheme to existing constrained approaches and prove that annealing the temperature during training effectively decreases upper bound on the acyclicity of the smooth orientation matrix (Theorem
We perform thorough experimental comparison of COSMO and competing structure learn ing approaches (Section 5). The empirical results report how COSMO achieves comparable structure recovery performances in significantly less time. Further, we highlight how COSMO consistently outperforms previous partially unconstrained structure learning proposals. For reproducibility purposes, we release the necessary code replicate experiments.
In the following, discuss related works in Section and report the necessary background on graph theory and structure learning in Section B3] Then, we introduce COSMO and our original contributions in Section 4 Finally, we report and discuss empirical analysis in Section
2

---

Published as a conference paper at ICLR 2024
Table 1: Summary comparison of our proposal, COSMO, with competing approaches. We propose parameterization that enables unconstrained learning of an acyelic graph without trading off on the adjacency matrix rank or the exactness of the acyclicity constraint. To express computational complexity, we define d as the number of nodes and k as the maximum length of iterative approaches. [i]: NOCURL requires a preliminary solution obtained by partially solving cubic-expensive con- strained optimization problem.
Method
Complexity
Constraint
NOTEARS Zheng et al. 2018
(
Exact
DAGMA Bello ct al. 2022
(
Exact
NOBEARS (Lee et al. 2019
(kk))
Approximated
TMPI Zhang et al. 2022
(kk))
Approximated
NOCURL Yu et al.
2021
( f
Partial
COSMO
(dd2)
None
RELATED WORKS
Both combinatorial and score-based structure learning approaches must deal with the number of possible DAGs, which grows exponentially with the number variables (Spirtes & Zhang, 2016). To avoid this issue, more recent line of research treats the space of DAGs continuous space and addresses the structure learning problem as an optimization task (Vowels et al., 2022). In this section, we report related works aiming to improve, approximate. or avoid altogether the constrained optimization scheme firstly introduced by NOTEARS (Zheng et al. (018)) as we summarize in Table Given the shared intuition, also briefly diseuss classical and continuous order learning approaches.
Low-Rank Approximation. Several works extended NOTEARS by assuming that the adjacency matrix of the underlying graph does not have full rank either to reduce the number of trainable parameters (Fang et al., 2023) or to improve computational complexity (Lopez et al. 2022). this work, we deal with possibly full-rank matrices and do not directly compare with low-rank solutions.
Constraint Reformulation. NOBEARS (Lee al. 2019) proposes estimate the acyclicity con- straint by approximating the spectral radius of the adjacency matrix. Given a maximum number k of iterations, the constraint can then be evaluated
graph with d nodes in O(kd²) time. Similarly, TMPI (Zhang et al. 2022) proposes an iterative approximation of the constraint that also results in O(kd)) computational complexity. Finally, DAGMA (Bello et al. (2022) reformulates the acyclicity constraint as the log degeemterminant of a linear transformation of the adjacency matrix. While still asymptotically cubic in complexity, the use of -og-determinant significantly faster in practice because of widespread optimizations in common linear algebra libraries (Bello et al. 2022 19)..
Unconstrained Methods. To tackle the exponential size of the space of DAGs, several works follow
the intuition of separately fitting the graph orientation and its adjacencies (Friedman & Koller 2003
Teyssier & Koller, 2012) Bernstein et al. 2020 Deng al.) (023). In the context of differentiable
approaches, of particular interest are causal discovery methods that avoid acyclicity constraints by
fitting distribution over permutations parameterized either explicitly (Cundy al.) [022]] 2 Charpentier et al. 2022) or through their topological ordering (Zantedeschi et al. (022). While sharing the overall intuition of fitting DAGs through their orientation, we propose a parameterization and an optimization scheme to fastly optimize acyclic adjacencies that could be then easily integrated in causal discovery methods. For this reason, we focus our analysis on competing acyclic optimization methods. However, we discuss the comparison with more complete approach, namely DAGUERREOTYES (Zantedeschi ctt al., 2022), in Appendix 19 Always in the context of causal discovery, ENCO (Lippe et al. 2022) decouples DAG in an adjacency matrix and edge orientation matrix. The authors explicitly parameterize the orientation matrix and prove that converges to an acyclic orientation whenever the training dataset contains sufficient number of interventions. Our structure learning proposal tackles instead non-intervened datasets and ensures acyclicity by construction. Similarly to us, NOCURL (Yu et al. 2021) proposes model that decouples the topological ordering from the adjacencies of an
3

---

Published as a conference paper at ICLR 2024
acyclic graph. However, the proposed optimization schemes are significantly different. Firstly, their approach extracts the nodes ordering from a preliminary solution obtained by partially solving the NOTEARS constrained optimization problem. Then, they fix the ordering and unconstrainedly optimize only the direct adjacency matrix. On the other hand, COSMO jointly learns priorities and adjacencies avoiding entirely acyelicity evaluations. Finally, GOLEM proposes unconstrained optimization problem regularized with the NOTEARS acyclicity constraint, which still needs to be evaluated on each optimization step and might nonetheless lead to cyclic DAGs (Ng et al. We report further discussion on the theoretical comparison with ENCO and NOCURL in Appendix and we carefully empirically compare with NOCURL in Section 5
BACKGROUND
Graph Theory. A directed graph is a pair D = (V, of vertices V (1,....,d d) and ares between them A C V x V. A directed acyclic graph (DAG) is a directed graph whose ares follow a strict partial order on the vertices. In a DAG, the parents of a vertex E V are the set of incoming nodes such that pa(v) = {u E V (u. v) E A} (Bondy & Murty
2008 We represent a directed graph as a define
binary adjacency matrix A E (0, 1) ]ddd where Aur 7 weighted adjacency matrix as the real matrix W E R Raxd where W + 0 + (u,v) E A.
(u, A. Similarly,
Structure Learning. A Structural Equation Model (SEM) models data-generating process as set of functions/ / - where fe= R Ripa(e)] - R for each variable v E V in the DAG (Pearl 2009). Given a class of functions F and a loss L, NOTEARS (Zheng et al. 2020) formalizes non-linear acyclic structure learning through the following constrained optimization problem
min L(F) s.t. tr(e (r(eW(/)w)())) -d= 0, FEF
where Wff) R is the adjacency matrix representing parent relations between variables in f. In particular, the constraint equals zero if and only if the adjacency matrix W() is acyelic. The authors propose to solve the problem using the Augmented Lagrangian method (Nocedal & Wright (999),, which in turn requires to solve multiple unconstrained problems and to compute the constraint value each optimization step. Notably, any causal interpretation of the identified ares depends on several assumptions on both the function class F and the loss function L (van de Geer & Bühlmann 2013 Loh & Bühlmann, 2014).
LEARNING ACYCLIC ORIENTATIONS WITH COSMO
In Subsection 4. we propose to parameterize weighted adjacency matrix function of direct matrix and a smooth orientation matrix. In this way, we effectively express the discrete space of DAGs in continuous and differentiable manner. Then, in Subsection 4.2 we introduce COSMO, an unconstrained optimization approach to learn acyelic DAGs. Furthermore, we prove an upper bound on the acyclicity of the smooth orientation matrix that connects our formulation to constrained approaches. To ease the presentation, we initially assume linear relations between variables. By doing so, the weighted adjacency matrix is the unique parameter of the problem. However, as with previous structure learning approaches, our proposal easily extends to non-linear models by jointly optimizing a non-linear model and an adjacency matrix either weighting or masking variables dependencies. We report one possible extension of COSMO to non-linear relations in Appendix
4.1 SMOOTH ACYCLIC ORIENTATIONS
To continuously represent the space of DAGs with d = |V| nodes, we introduce a priority vector Rd on its vertices. Consequently, given the priority vector p and strictly positive threshold
E >>00 0, we define the following strict partial order <p.e on the vertex set V
V(u, v) E V V: u "p.e"" U - P. - P 2..
(2)
In other terms, a vertex u precedes another vertex if and only if the priority of v is sufficiently larger than the priority of the vertex Notably, with a zero threshold E 0. the relation would be symmetric and thus not a strict order. On the other hand, whenever E is strictly positive, we can represent subset of all striet partial orders sufficient to express all possible DAGs.
4

---

Published as a conference paper at ICLR 2024
3
X3
3
0
Figure 2: (left) With infinite temperature, the sigmoid function is constant and connects all vertices. (center) Given two nodes, for positive temperatures the smooth orientation matrix has larger values on the ares respecting the priority ordering. (right) In the limit of the temperature to zero, the smooth orientation matrix contains non-zero entries if and only if the are respects the order, i.e., it directs node to another with sufficiently higher priority.
Lemma 1. Let W E Rddd be a real matrix. Then, for any & >> O, W is the weighted adjacency matrix of a DAG if and only if it exists priority vector p E Rd and real matrix H R such that
W H sp.c>
where Tp. E [0, 1]dxd is binary orientation matrix such that
{1 if <p,c v otherwise,
T sp. [we]
(4)
for any u. U E V.
Proof. We report the proof in Appendix A.1
While priority vectors enable the representation of strict partial orders in a continuous space, the construction of the orientation matrix still requires the con-differentiable evaluation of the inequality between priority differences from Equation 2 To this end, we approximate the comparison of the difference against the threshold e, using a tempered sigmoidal function. We refer to such approximation of the orientation matrix the smooth orientation matrix.
Definition 1 (Smooth Orientation Matrix). Let p R2 be a priority vector, E positive threshold, and l > 0 be a strictly positive lemperature. Then, the smooth orientation matrix of the strict partial order <p.s the real matrix St,s R such that, for any u, E it holds
be strictly
(p) = oz,e (p. - Pu),
(5)
where ot,s is the e-centered tempered sigmoid, defined as
(x)) = e-
Intuitively, the threshold E shifts the center of the sigmoid and breaks the symmetry whenever two variables approximately have the same priority. The temperature parameter > 0 regulates instead the steepness the sigmoid. Because the asymmetry introduced by the threshold, in the limit of the temperature zero, the zero-entries of smooth orientation matrix coincide with the zero-entries of the corresponding orientation matrix (Figure 2 Therefore, we prove that any directed acyclic graph can be represented as the element-wise product of a directed adjacency matrix and smooth orientation. Further, any directed graph resulting from this decomposition is acyclic.
Theorem 1. Let W E d be real matrix. Then, for any O, W is the weighted adjacency matrix of a DAG if and only if it exists a priority vector p E Rdd and real matrix H E xd such that
W H lim (p. (-0
where (p) the smooth orientation matrix of <p.--
Proof. We report the proof in Appendix A.2
5

---

Published as a conference paper at ICLR 2024
4.2 LEARNING ADIACENCIES AND ORIENTATIONS
Given our definition of smooth acyclic orientation (Definition D. we can effectively parameterize the space of DAGS continuous function of a weighted adjaceney matrix H Raxd and a priority vector P E IR. Therefore, the computational complexity of our solution reduces to the construction of the adjacency matrix W H St,e (P), which can be achieved in Ocd2) time and space per optimization step by computing each are as W uu = Hav o((p. P - =)/ In the literature, NOCURI. proposed similar model where each are has form W.. ur H.r ReLU(p, - P.). Despite the similarity, NOCURL is trained with significantly different optimization scheme, as we also discuss in Appendix D.2 In fact, to avoid a significant performance drop, their formulation requires preliminary solution from a constrained optimization problem and does not jointly learn the parameters corresponding to our adjacencies and priorities. In the following, we describe how COSMO effectively reduces to an unconstrained problem and avoids evaluating acyclicity altogether.
Temperature Annealing. The smooth orientation matrix St.e(p) represents an acyclic orientation only in the limit of the temperature - O. Nonetheless, the gradient loss vanishes whenever the temperature tends to zero. In fact, for an arbitrary loss function L, we can decompose the gradient of each component P, of the priority vector follows
ac(w)
OC(W)
Ow we
aC(W)
@we
(8)
0pu
OW
0p.
OW aM 0P.
DEV
aw
H
Tee (P. -P.](1 or. (P. -Pa)) -0.
(9)
0p. Heu
aw
Pp. 1 (P. - -Pa)).
(P,
(10)
Therefore, by property of the sigmoidal function ore it holds that both www./Op. and OW @Wow/Op /
tend to zero for - 0. To handle this issue, we tackle the optimization problem by progressively reducing the temperature during training. In practice, we perform cosine annealing from initial positive value fuan to a significantly lower target value fend R O. We further motivate our choice by showing the existence of an upper bound on the acyclicity of the orientation matrix that is monotone increasing function of the temperature. Therefore, temperature annealing effectively decreases the acyclicity upper bound during training of the smooth orientation and, consequently, of the adjacencies.
Theorem 2.
Let p E Rd be priority vector, e >> O be a strictly positive threshold, and t > 0 be strictly positive temperature. Then, given the smooth orientation matrix St,e (p) E RRddd it holds
h(St..c (p)) S do -1,
(11)
where (D)) = (r(eS...(P)) - d is the NOTEARS acyclicity constraint and a = o-E/())
Proof. We report the proof in Appendix B
Direct Matrix Regularization. To contrast the discovery of spurious ares we perform feature selection by applying LI regularization on the adjacency matrix H. Further, during the annealing procedure, even if a vertex u precedes v in the partial order <p.c. the weight of the opposite are " + u in the smooth orientation matrix will only approximately zero. Therefore, sufficiently large values of the weighted adjacency matrix H. might still lead to undesirable cyclic paths during training. To avoid this issue, we regularize the L2-norm of the non-oriented adjacency matrix.
Priority Vector Regularization. Other than for small temperature values, the partial derivatives in Equations 9 and 10 tend to zero whenever the priorities distances lp, p.l tend to infinity. Therefore, regularize the 1.2-norm of the priority vector. For the same reason, we initialize each component from the normal distribution P., ~ (O, e² 2/2,, so that each difference follows the normal distribution P. P,, ~ ~N(0 We provide further details on initialization in Appendix A.3
Optimization Problem.
We formalize COSMO as the differentiable and unconstrained problem
min HeRexd.ERe
L(H St.( A H| + X H 2 + Xp lllll
(12)

---

Published as a conference paper at ICLR 2024
where λ1, A2 the regularization coefficients for the adjacencies and the priorities. As the regularization coefficients À = {^1, 22, Ap). the initial temperature stat, the target temperature fend- and the shift E are hyperparameters of our proposal, whose choice discuss in Appendix C.4 Notably, Theorem 2 can also guide the choice of the final temperature value and the shift to guarantee maximum tolerance the acyclicity of the smooth orientation matrix.
EXPERIMENTS
We present an experimental comparison of COSMO against related acyclic structure learning ap- proaches. Our method operates possibly full-rank graphs and ensures the exact acyclicity of the solution. Therefore, we focus on algorithms providing the same guarantees and under the same conditions. Namely,
confront with the structure learning performance and execution time of NOTEARS (Zheng et al. 2018). NOCURL (Yu et al. 202T), and DAGMA (Bello et al.) 2022). As previously discussed, NOCURI. proposes similar model with substantially different optimization scheme. To highlight the importance of both our parameterization and optimization scheme, we also compare with an entirely unconstrained variant of the algorithm where we directly train the variables ordering without any preliminary solution. In the results, we refer to this variant as NOCURL-U.
We base our empirical analysis on the testbed originally introduced by Zheng et al. (018)) and then adopted as a benchmark by all followup methods. In particular, we test continuous approaches on randomly generated Rénös-Rényi (ER) and scale-free (SF) graphs of increasing size and for different exogenous noise types. For each method, we perform structure learning by minimizing the Mean Squared Error (MSH) of a model on a synthetic dataset using the Adam optimizer (Kingma & Ba 2015). In Appendix C. we report further details on the implementation of COSMO, the baselines, and the datasets. We include the code to reproduce our experiments in the Supplementary Materials !
5.1 EVALUATION OVERVIEW
In line with previous work, we retrieve the binary adjacency matrix by thresholding the learned weights against small threshold w = 0.3 (Zheng et al. 2018). While COSMO guarantees the solution to be acyelic, we maintain the thresholding step to prune correctly oriented but spurious ares. Then, we measure the Normalized Hamming Distance (NHD) between the solution and the ground-truth as the sum of missing, extra, or incorrect edges divided by the number of nodes. In general, testing weights against a fixed threshold might limit the retrieval of significant in the true model (Xu et al., 2022). For this reason, we also compute the Area under the ROC curve (AUC), which describes the trade-off between the True Positive Rate (TPR) and the False Positive Rate (FPR) for increasing values of the weight threshold (Heinze- -Dezze et al. 20188. Due to space limitations, we only report the main body the AUC results, which is the most comprehensive score. We provide detailed results for other metrics, including NHD, in Appendix
with small coefficients
5.2 RESULTS DISCUSSION
By looking at the AUC of the learned graphs, we observe that COSMO consistently achieves results that are comparable and competitive with those from constrained-oppiimzzzi DAGMA or NOTEARS across different graph sizes and noise types (Table[). This empirically confirms the approximation properties of COSMO, which can reliably discover DAGs without resorting to explicit acyclicity constraints.
solutions such as
Furthermore, COSMO performs better than NOCURL on most datasets. We recall that the latter the only existing structure learning approach combining constrained and unconstrained optimization. As pointed out in Yu et al. (2021), ) we also observe that the discovery performance of NOCURL drops when optimizing the variable ordering instead of inferring it from preliminary solution. The fact that COSMO outperforms NOCURL-U on all datasets highlights the substantial role and effect of our smooth orientation formulation and our optimization scheme for learning the topological ordering of variables from data in an unconstrained way. Overall, our proposal achieves, on average, the best or the second-best result for the AUC metric across all the analyzed datasets and correctly classifies ares also for large graphs (Figure D. Further, as we extensively report in Appendix E for different graph
https:: thub. rma dda ccamo

---

Published as a conference paper at ICLR 2024
Table 2: Experimental results on linear ER-4 acyclic graphs with different noise terms and sizes. For each algorithm, we report mean and standard deviation over five independent runs of the AUC metric and the time in seconds. We highlight in bold the best result and in italic bold the second best result. The reported duration of NOCURL includes the time to retrieve the necessary preliminary solution using an acyclicity constraint. We denote as NOCURL-U the quadratic version of NOCURL. Complete results on additional metrics and graph types are in Appendix
Gauss
Exp
Gumbel
d
Algorithm
AUC
Time
AUC
Time
AUC
Time
COSMO
0.984 0.02
88 1 2
0.989
0.01
89
0.914 0.10
87 1
DAGMA
0.985 + 0.01
781 192
0.986
0.02
744
75
0.973 = 0.02
787 + 86
30
NOCURI.
0.967 0.01
822 + 15
0.956
0.02
826
24
0.915 + 0.04
826 + 17
NOCURL-U
0,694 0.06
226 = 5
0.694 0.05
212
0.678 = 0.05
212 + 5
NOTEARS
0.973 1: 0.02
5193 170
0,966
0.03
284 5579
0.981 0.01
5229 1: 338
0,961 0.03
0.985
0.973
0.01
98
DAGMA
0,982 0.01
660
141
0.986
0.01
733 109
0.986
0.01
858
101
100
NOCURL
0.962 0.01
1664
14
0.950 0.02
1655
28
0.962
0.01
1675
34
NOCURL-U
0.682 0.05
267
10
0.693 1 0.05
242
0.663
0.04
247
NOTEARS
0.963 0.01
11000
339
0.972 0.01
10880
366
0.00 0.969
11889
343
COSMO
0.933 0.01
436 81
0.986
390
0,982
106
500
DAGMA
0,980 0.00
2575
2853
NOCURL-U
0.683
1546
1488
0.05
1342
209
The Wegnts
Predicted Weights
Weights Difference
100
1.0
10
05
05
0.5
90
00
00
-05
Figure 3: Visualization of the weighted adjacency matrix learned by COSMO (ER4, Gaussian noise, 100 nodes) against the ground truth. We also report the difference between the ground-truth and the learned weights. By thresholding the learned weighted adjacency matrix, COSMO correctly classifies most true (TPR = 0.96) and non-existing ares (FPR = 0.01), resulting in a limited number of errors (NHD 0.93) due to the narrow difference in the retrieved weights.
size and noise terms, our non-linear extension obtains comparable performances with DAG in significantly less time.
Unsurprisingly, due its quadratic computational complexity, COSMO significantly faster than constrained methods all datasets, especially for increasing graph sizes. Notably, despite employing carly stopping conditions for convergence, all competing methods incur the cost of solving multiple optimization problems with higher computational cost per step (Figure Z. In particular, while the unconstrained variant NOCURL-U has a comparable per-epoch average time cost, for a substantially worse graph recovery performance, NOCURL overall pays the need for preliminary solution computed with an acyclicity constraint. Therefore, already for graphs with 500 nodes, only COSMO, DAGMA, and NOCURL-U return solution before hitting our wall time limit. Finally, observe that the cubic computational complexity of DAGMA significantly emerges when dealing with large graphs. Therefore, despite the effective underlying optimizations the log-determinant computation, DAGMA'S 's acyclicity constraint still affects scalability.
Given the proposed parameterization, COSMO requires particular care in the choice of the regulariza- tion hyperparameters. In particular, we carefully considered the importance of regularizing the priority vector, which constitutes one of our main differences with previous structure learning approaches.
8

---

Published as a conference paper at ICLR 2024
Table Ablation test the priority regularization term Xp on DAGS with different noise terms and sizes. We name the configuration without priority regularization COSMO-NP and report mean and standard deviation over five independent runs. For each configuration, the best result in bold.
Graph
d
Algorithm
NHD
TPR
AUC
COSMO
0.867 + 1.01
0.953 + 0.04
0.984 + 0.02
30
COSMO-NP
1.893 + 0.92
0.870 + 0.07
0.937 + 0.05
ER4
COSMO
1,388 + 0.69
0.917 + 0.04
0.961 + 0.03
100
COSMO-NP
1.570 0.56
0.935 0.04
0.974 + 0.02
COSMO
4.087 + 1.12
0.838 0.06
0.921 + 0.04
30
COSMO-NP
4.153 + 2.38
0.819 = 0.12
0.885 - 0.09
ER6
COSMO
9.476 3.01
0.771 0.08
0.911 0.05
100
COSMO-NP
9.804 + 2.90
0.848 - 0.05
0.941 + 0.03
We found that our hyperparameter search procedure consistently returned relatively low priority regular- ization values (Ap ~ le-3). However, while might benefit structure learning for larger graphs, ablating priority regularization results a non-negligible per- formance drop for smaller graphs (Table 3).
Alywthin
COSMO
3P
Tire 04
Another important aspect of our proposal the tem- perature annealing procedure. While proved that in the limit of the temperature to zero the resulting model will be acyclic, the same procedure might lead to zero- gradients effectively impeding the optimiza- tion. In Appendix G show that COSMO consis- tently reaches acyclic solutions before completely annealing the temperature and thus that optimization continues even in the last epochs. Finally, we also find empirical evidence that our L2 regularization scheme on the directed adjacency matrix effectively contains acyclicity during training (Appendix F).
1011
107
Figure 4: Average duration of a training epoch for an increasing number of nodes on five independent runs on random ER-4 DAGs.
6 CONCLUSION
We introduced COSMO, unconstrained and continuous approach for recovering DAGs from data. Our novel definition of smooth orientation matrix ensures the acyclicity of the solution without requiring the evaluation of computationally expensive constraints. Furthermore, we prove that annealing the temperature of our smooth acyclic orientation corresponds to decreasing an upper bound on the widely adopted acyclicity relaxation from NOTEARS. Overall, our empirical analysis showed that COSMO performs comparably to constrained methods in significantly less time. Notably, our proposal significantly outperforms the only existing work partially optimizing the space of DAGs, NOCURL, and its completely unconstrained variant NOCURL-U, Experimental results highlight the role of our parameterization, which does not incur the necessity of preliminary solutions and provably returns a DAG without ever evaluating acyclicity.
In recent years, several authors debated using continuous acyclic learners as full-fledger causal discovery algorithms (Reisach et 2021 Kaiser & Sipos, 2022 Ng et b). In this context, our empirical analysis of COSMO shares the same limitations of existing baselines and, exactly like them, might not be significant in the causal discovery scenario. However, acyclic optimization techniques are a fundamental component of continuous discovery approaches (Brouillar et al.) 2020 Lorch et al. 2022). By reducing by order of magnitude the necessary time to optimize acyclic causal graph, COSMO opens up more scalable continuous causal discovery strategies without sacrificing demonstrated in this work the theoretical guarantees DAGS approximation capabilities.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
The work has been partially supported by the EU H2020 TAILOR project (n.952215).
REFERENCES
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: Learning DAGs via M-matrices and Bog-Deetrrmnnnt Acyclicity Characterization. In Advances in Neural Information Processing Systems, 2022.
Daniel Bernstein, Basil Saeed, Chandler Squires, and Caroline Uhler. Ordering- causal structure learning in the presence of latent variables. In International Conference Artificial Intelligence and Statistics, pp. 40980099900011 4108. PMLR, 2020.
Adrian Bondy and M. Ram Murty. Graph Theory. Graduate Texts in Mathematics. Springer- Verlag, 2008. ISBN 978-1-84628-969-,, 9781846289699
URL fttts::/www.springeee. com/ gp/book/
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexan- dre Drouin. Differentiable causal discovery from interventional data. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.). Advances in Neural In- formation Processing Systems, volume 33, pp. 21865-21877. Curran Associates, Inc., 2020. URL https://proceeedinggg neurips cwwpppprrrcppppprreeeeepep files/ files//apeer2220/ File/ FLLefFrnaaaaa0110555222211101122
Bertrand Charpentier, Simon Kibler, and Stephan Günnemann. Differentiable dag sampling. arXiv preprint wxXiv:2203008509 2022.
Gonçalo Correia, Vlad Niculae, Wilker Aziz, and André Martins. Efficient marginalization of discrete and structured latent variables via sparsity. Advances in Neural Information Processing Systems, 1789-1 2020.
Chris Cundy, Aditya Grover, and Stefano Ermon. Bed nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095-7110, 2021.
Chang Deng, Kevin Bello, Bryon Aragam, and Pradeep Kumar Ravikumar. Optimizing notears objectives via topological swaps. International Conference on Machine Learning, pp. 7563- 7595. PMLR, 2023.
Zhuangyan Fang. Shengyu Zhu, Jiji Zhang, Yue Liu, Zhitang Chen, and Yangbo He. On low-rank directed acyclic graphs and causal structure learning. IEEE Transactions on Neural Networks and Learning Systems, pp. 1-14, 2023. doi: 10. ALLLLO/NNNLS.2222...22222.
Nir Friedman and Daphne Koller. Being bayesian about network structure. bayesian approach to structure discovery in bayesian networks. Machine learning. 50:95-125, 2003.
Alexander Hägele, Jonas Rothfuss, Lars Lorch, Vignesh Ram Somnath, Bernhard Schölkopf, and Andreas Krause. Bacadi: Bayesian causal discovery with unknown interventions. International Conference on Artificial Intelligence and Statistics, pp. 1411-1436. PMLR, 2023.
Christina Heinze-Deml, Marloes H Maathuis, and Nicolai Meinshausen. Causal structure learning.
Annual Review of Statistics and Its Application, 5:371-391, 2018.
William Vallance Douglas Hodge. The theory and applications of harmonic integrals. CUP Archive, 1989.
Mareus Kaiser and Maksim Sipos. Unsuitability of notears for causal graph discovery when dealing with dimensional quantities. Neural Processing Letters,
2022.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, ICLR 2015, San Diego CA. USA, May 7-9, 2015, Conference Track Proceedings, 2015. TRI.tttp::// http://pp//aaa...oo///a/////////
6980
10

---

Published as a conference paper at ICLR 2024
Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient- based neural DAG learning. In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, OpenReview.net, 2020. URL https://
openreview www.eeeeeedddeeiie......
Hao-Chih Lee, Matteo Danieletto, Riccardo Miotto, Sarah T Chern,, and Joel T Dudley. Scaling structural learning with no-bears to infer causal transcriptome networks. In PACIFIC SYMPOSIUM ON BIOCOMPUTING 2020, pp. 391-402. World Scientific, 2019.
Phillip Lippe. Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclicity
constraints. In International Conference Learning Representations, 2022. URL https
//openreview net. Lttp://wwwe...tterooommdddeedddd
Po-Lingg Loh and Peter Bühlmann. High-dimensional learning of linear causal networks via inverse covariance estimation. The Journal of Machine Learning Research, 15(1):3065-3105, 2014.
Romain Lopez, Jan- Christiitt Hütter, Jonathan Pritchard, and Aviv Regev. Large-scale differentiable causal discovery of factor graphs. Advances in Neural Information Processing Systems, 35: 19290-19303, 2022.
Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. Amortized inference for causal structure learning. In Advances in Neural Information Processing Systems, 2022. URL
-MMXX
Ignavier Ng. AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. In Advances in Neural Information Processing Systems, volume 33, pp. 17943-17954. Curran Associates, Inc., a. URL htttts::::prrcceeiiiiii/o/ neurisss cc/ Abstract .tt.attthtmllll
Ignavier Ng. Biwei Huang, and Kun Zhang. Structure learning with continuous optimization: A sober look and beyond. Number arXiv:2304.02146 arXiv, b. doi: 10445500arXxx2234401146
URL http:: wrw.vvorrrrooo ab35/ 02146
Ignavier Ng, Sébastien Lachapelle, Nan Rosemary Ke, Simon Lacoste-Julien, and Kun Zhang. On the convergence of continuous constrained optimization for structure learning. In International Conference on Artificial Intelligence and Statistics, pp. 8176-8198 8198. PMLR, 2022a.
Ignavier Ng, Shengyu Zhu, Zhuangyan Fang, Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient-based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pp. 424 432. SIAM, 2022b.
Vlad Niculae, Andre Martins, Mathieu Blondel, and Claire Cardie. Sparsemap: Differentiable sparse structured inference. In International Conference on Machine Learning, pp. 3799-3808. PMLR, 2018.
Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 1999.
Judea Pearl. Causality. Cambridge university press, 2009.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated dag! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34:27772-277 27784. 2021.
Peter Spirtes and Kun Zhang. Causal discovery and inference: concepts and recent methodological advances. In Applied informatics, volume 3, pp. 1-28. SpringerOpen. 2016.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search, Second Edition. Adaptive computation and machine learning. MIT Press, 2000.
Marc Teyssier and Daphne Koller. Ordering- based search: A simple and effective algorithm for learning bayesian networks. arXiv preprint m2xiiv/200. 1429, 2012.

---

Published as a conference paper at ICLR 2024
Sara van de Geer and Peter Bühlmann. 20-penalized maximum likelihood for sparse directed acyelic graphs. The Annals of Statistics, 41(2):536-561 2013.
Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden. D'ya like dags? survey on structure learning and causal discovery. ACM Computing Surveys, 15(4):1-36, 2022.
Danru Xu, Erdun Gao, Wei Huang, and Mingming Gong. the sparse dag structure learning based on adaptive lasso. arXiv preprint arXiv:2209.0222 2222229999026 2022.
Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji. Dags with curl: An efficient dag structure learning approach. In International Conference on Machine Learning, pp. 12156-12166. PMLR, 2021.
Valentina Zantedeschi, Luca Franceschi, Jean Kaddour, Matt Kusner, and Vlad Niculae. Dag learning the permutahedron. In The Eleventh International Conference Learning Representations, 2022.
Zhen Zhang, Ignavier Ng, Dong Gong, Yuhang Liu, Ehsan Abbasnejad, Mingming Gong, Kun Zhang, and Javen Qinfeng Shi. Truncated matrix power iteration for differentiable dag learning. Advances in Neural Information Processing Systems, 35:1 83390-18402 2022.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric DAGs. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pp. 3414-3425, PMLR, 2020. URL httttp:://proceediggg.. mlr ISSN: 2640-3498.
press/v108// zheng20a. html
12

---

Published as a conference paper at ICLR 2024
A DEFERRED PROOFS
A.1 PROOF OF LEMMA I
Lemmal 1
Let W E Raxd be real matrix. Then, for any ε W the weighted adjacency
matrix of a DAG if and only if it exists priority vector p E RR and real matrix H E Rdxa such that
W = H T pp..
(13)
where T sp. E {0, 10,1ddxd is binary orientation matrix such that
{
if u <p.s v otherwise,
T Tp. [wv]
(14)
for any u,u V.
Proof. Firstly, we prove the existence of priority vector and an adjacency matrix H for each weighted acyclic matrix W of directed acyclic graph = (V, A). Being DAG, the arcs follow strict partial order < on the vertices V = (I..... Therefore, it holds that
A {u, v) v).
(15)
Consequently, for an arbitrary topological ordering of the variables и: V - fl..... d}, which always exists on DAGs, we define the vector p IR such that
P = ET)
(16)
Given the following implications
K T(b) > T(u)
(17)
P, - Pa = =(T(O) - =(u)) 2
(18)
= pp
(19)
holds that the order ip.s contains the order <. Finally, we can define the adjacency matrix as H = W, where W = H sp. holds since p.. [u, w] = 0 only if (u, v) & A.
To prove that any priority vector p and adjacency matrix H represent DAG, first notice that, since the ares follow a strict partial order, the orientation T <p. is acyelic. Then, by element-wise multiplying any matrix H we obtain sub- aub-graph of DAG, which acyelic by definition.
A.2 PROOF OF THEOREM
Let W Raxd be real matrix. Then, for any > O, W is the weighted adjacency matrix of a DAG if and only if it exists a priority vector p E Rd and real matrix H E Bdxd such that
Theorem I
W Ho lim (p),
(20)
(+0
where St,e (p) the smooth orientation matrix of <p.s-
know that for any acyclic weighted adjacency matrix W there exist priority vector p and a real matrix H such that W H o T Tep.. Further, by Definition the inner limit of
Proof. By Lemma
Equation 20 solves to
P. Pa E
lim (p) t00
1/2 P. P, = E
(21)
P. < E.
Therefore, define H' E Rdxd such that
2H
P, - P. =E otherwise.
HGp
H
(22)
H
13

---

Published as a conference paper at ICLR 2024
from which
W = Ho T Spc
HH' lim (p). t-30
(23)
Then, to prove the counter-implication of Theorem
we notice that
lim (p. = P. Pa <E
Ap. v.
(24)
Therefore, since the smooth orientation contains arc if and only if the vertex respect the strict partial order <p,s. acyclic. Consequently, as Lemma | the element-wise product with an acyclic matrix results in sub-graph of a DAG, which is also acyelic by definition.
A.3 PRIORITY VECTOR INITIALIZATION
By independently sampling each priority component from Normal distribution M(H, s² /2), each difference consequently sampled from distribution (0,0,²) s²). Therefore, we seek a value for the standard deviation that maximizes the partial derivative
ow
Hus
(P. o0. PP.).
(25)
0p.
for arbitrary vertices u, v. Given the definition of the campered-shifted sigmoid function, this object has maximum P, P, Therefore, by setting the variance as s² =², maximize the density function the point P., Р., & in the distribution s²).
B SMOOTH ACYCLIC ORIENTATIONS AND THE ACYCLICITY CONSTRAINT
In this section, we present the proof for the upper bound on the acyclicity of smooth acyclic orientation matrix. To this end, introduce two auxiliary and novel lemmas. Firstly, introduce lemma which binds the product of sigmoid on a sequence of values with zero sum (Lemma Then, we introduce another lemma on the sum of the priority differences in cyclic path (Lemma Finally, we are able to prove the acyclicity upper bound from Theorem 2
Lemma 2. (Sigmoid Product Upper Bound) Let {a4} be a sequence of n. real numbers such that
x =
Then, for any temperature > 0 and shift 2 0, holds that
II (a)) sa"
where a o(-s/t) the value of the tempered and shifted sigmoid in zero.
Proof. Before starting, invite the reader notice that, for any temperature / > 0, if the of the sequence {a;} } is zero, then also the sequence {æ./t} sums to zero. Therefore, we omit the temperature the following proof, and assume to divide beforehand all elements of the sequence by the temperature t. Further, we explicitly denote the shifted sigmoid by using the notation a c).
Firstly, we formulate the left-side of the inequality as
Een ai-e
II oo 5) II -
e
-
+ 6m--5"
Similarly, we rewrite the right side as
)
n
a" == ====(=-)"=
(1+e-s)*
14

---

Published as a conference paper at ICLR 2024
Therefore, proving the left-side smaller equal than the right-side, reduces to proving the left- denominator larger than the rightdenommnator. Formally,
+l+e c 2 (1+
or equivalently, by applying the logarithmic function,
log(1 + + 2 nlog(1 + +e=e).
(26)
To further ease the notation, we refer to the left side of inequality 26 as the target function
T(x) = Σ
log (1 + emi-O).
partieular, to prove 26 we show that
min T(x) n log(1 +e=s),
(27)
for = 0, which is the only stationary point due to the convexity of the target function.
Without loss of generality, we derive the partial derivative of the component oi on the target function T(x). To constraint the components sum to zero, we consider the components {z.} for i as free, and then = - 2ing 2g as a function of the remaining. The choice of æ1, is independent from the components ordering, and thus applies any possible pair. Consequently,
8T(x) O(log(1 + + log( -1 -ina #i 5) + Ex=3 log(1 + e((-()
(28)
8z
de
a(log(1 + + log(1 + Diss
(29)
dx
0 c)) o
x; E).
(30)
i=3
Since o--)) = (--)), the equation is satisfied, for any component 2,, by = Ö.
We finally prove Inequality 26] by showing that the value of the target function T(x), in its only stationary point æ equals the bound. Formally,
т(Ò) Σ log + 5)
(31)
==1
n log(( te
(32)
Lemma 3. (Sum of Differences in Cycle) Let {p:} be sequence of n. + 1 real numbers such that Pi = Pa41- Then, let { be sequence of n numbers such that 6, = Pi+1 Pi. Then,
Σ
=.
(33)
Proof. The proof immediate from the following sequence of equations:
8
Σ
pi+1 - Ps = -
Pi Pt Pn+ 0.
15

---

Published as a conference paper at ICLR 2024
Theorem 2 Let R be priority vector, > 0 be strictly positive threshold, and t be
strictly positive temperature. Then, given the smooth orientation matrix (t,( E IRdxd it holds
t( (p)) S do
(34)
where ((S),:(D)) tr(eSt.(()) - d the NOTEARS acyclicity constraint and a o(-e/t).
Proof. The left side of Inequality 34
corresponds to the following infinite series
2
Ir(eP)
hd! lr(P(k)
d
km0
tr(p(k))
k!
where (k) is the matrix power defined as p
p
and p I.
By definition of matrix power, the u-th element on the diagonal of p equals
((k)
Σ
VIEV
k-2 II
Σ
Σ VE-IEV
Pu
P P
Intuitively, the u-th element on the diagonal of p(k) amounts to the sum of all possible paths starting and ending in the variable X,, Therefore, being the same node, the priority of the first and the last node in the path priorities sums to zero. For this reason, given Lemma holds that the product of the corresponding sigmoids is smaller equal than a* Therefore,
equal by definition. Consequently, by Lemma [3] the difference between the
A-2
(
Σ
Σ
Pum
II P
P 15
VvEV
V--IEV
Σ
Σ
V
VE-
Consequently, upper bound the trace of the orientation matrix power as
tr (P
((k) uu
I d A
Finally, able to prove the Theorem
tr(eP) d=
tr(p(k)) -d
(P()
k! t
da
where the last passage is due the Taylor series of the exponential function.
16

---

Published as a conference paper at ICLR 2024
C IMPLEMENTATION DETAILS
In this section, we discuss the significant aspects of our implementation. We run all the experiments on our internal cluster of Intel(R) (eo Gold 5120 processors, totaling 56 CPUs per machine. The necessary code and instructions are in the Supplementary Materials. We report details on the evaluation (.II). the data generation procedure CC2 and the models (C.4).
C.I EVALUATION PROCEDURE
We ensure a fair comparison by selecting the best hyperparameters for each implemented method on each dataset. We describe the hyperparameter space for each algorithm the following subsections. Firstly, we sample fifty random configurations from the hyperparameter space. Since the hyperparam- eter space of COSMO also includes temperature and shift values, we extract more hyperparameters (200 800). Due to the significant speedup of COSMO, hyperparameter searches take comparable amount time, with NOTEARS being significantly longer on small graphs as well. Then, we test each configuration to the average AUC value. Finally, we perform a validation step by running the best configuration on five random graphs.
five randomly sampled DAGs. We select the best hyperparameters according
To report the duration of each method, we track the time difference from the start of the fitting procedure up the evaluation process, excluded. For constrained methods, this includes all the necessary adjustments between different problems. For COSMO, it contains the annealing of the temperature between training epochs.
Following previous work, we recover the binary adjacency matrix A of the retrieved graph by thresholding the learned weights W with small constant wi 0.3. Formally, A | w.
C.2 SYNTHETIC DATA
We include in our code the exact data generation process from the original implementation of NOTEARS Therefore, the dataset generation procedure firstly produces DAG with either the Rrdōs-Rényi (ER) or the scale-free Barabási- Albert (SF) models. Then, it samples 1000 independent observations. In the linear case, the generator simulates equations of the form
f((x) W; + 24
(35)
where we sample each weight Wi, from the uniform distribution 4 44-2, -055) (0.5, 2) and each noise term 2; from either the Normal, Exponential (A = 1), or Gumbel (μ = 0, 3 1) distributions. In the non-linear case, simulate additive noise model with form
[(x)) (x)) Z4,
(36)
where di a randomly initialized Multilayer Perceptron (MLP) with 100 hidden units and the noise term 2, is sampled from the Normal Distribution N (0, 1).
C.3 NON-LINEAR RELATIONS
We generalize COSMO represent non- non-linear relations between variables. To case the comparison. we follow the non-linear design of NOTEARS-MI.P (Zheng et al. 2020) and DAGMA-MI.P (Bello et al. 2022). However, it is worth mentioning that this only of the possible approaches for non-linear relations and that COSMO parameterization could be easily extended to mask the input of neural network instead of the weights, as done by Ng al. (2022b) or Brouillarr al. (2020). Similarly to NOTEARS (Zheng et al., 2020), we model the outcome of each variable X,, with neural network fu: Rd R. where we distinguish between the first-layer weights H E Rdx for & distinct neurons, and the remaining parameters &,, We ensure acyclicity by considering each weight matrix H as the w-th slice on the first dimension of tensor H E Rdxdxh Intuitively, each entry H represents the weight from the variable X,, to the i-th neuron in the first layer of the
NOTTARS implementation is published with Apache license at notears
17

---

Published as a conference paper at ICLR 2024
Table 4: Hyperparameter ranges and values for NOTEARS, NOCURL, and DAGMA.
Hyperparameter
Range/Value
Hyperparameter
Range/Value
Learning Rate
(Ie-4, le-3)
Learning Rate
(le-3, 1e-2)
Regularization
(Ie-4, le-1)
A
(le-4, le-3)
Initial Penalty
1.0
X2
(le-3, 5e-3)
Penalty Factor
10.0
Xp
(Ie-3, 30-3)
Max Penalty
le+16
start
0.45
Multiplier
1.0
lend
(5e-4, le-3)
DAGness Tolerance
le-8
(Se-3, 2c-2)
Progress Rate
0.25
Max Iterations
100
(a) COSMO
(b) NOTEARS
Hyperparameter
Range/Value
Hyperparameter
Range/Value
Learning Rate
(le-4, le-3)
Learning Rate
(1e-4, 1e-3)
Regularization
(le-4, le-1)
Regularization
(1e-4, le-1)
Initial Path Coefficient
1.0
Multiplier
(10.0, 1000.0)
Path Decay Factor
0.1
Prethreshold
0.3
Log-Det Parameter
0.90, 0.8, 0.7)
Progress Rate
0.25
Max Steps
(d) NOCURL
(c) DAGMA
MLP fe. Then, we broadcast the element-wise multiplication of smooth orientation matrix on the hidden dimensions. Formally, model the structural equation of each variable X,,
fu(x) = (x(((x) T [H (p)) Lu),
(37)
where yo is an activation function and 9v is Multilayer Perceptron (MLP) with weights d.. By
applying each MLP fu to each variable X,,, we define the overall SEM as the function f: Rd - Rd
which depends the parameters T- (..). the weight tensor H, and the priority vector p. Therefore, we formalize the non-linear extension of COSMO as the following problem
min HERdXa SERR4x"_pRd__0 pER
C(H St.- 4b + >1 H||
Ao H| + Apl|p||2.
(38)
C.4 MODELS
Since we focus on the role of acyclic learners as a building block within more comprehensive discovery solutions, slightly detach from experimental setups considering such algorithms as standalone structure learners. Therefore, instead of dealing with full-batch optimization, we perform mini-batch optimization with batch size B 64. Similarly, instead of explicitly computing the gradient of the loss function, we implement all methods in Py Torch to exploit automatic differentiation. By avoiding differentiation and other overhead sources, the time expenses results comparable between our implementations and the results reported in the original papers. However, our implementation choices are common to works that employed NOTEARS et similia to ensure the acyclicity of the solution (Lachapelle et al.) 2020 Brouillard et al. 2020
directly
Lopez et al.) 2022).
By checking the convergence of the model, both NOTEARS, DAGMA, and NOCURL can dynamically stop the optimization procedure. On the other hand, COSMO requires fixed number of epochs in which to anneal the temperature value. a fair comparison, while we stop optimization problems after maximum of 5000 training iterations, baselines. Therefore, when sufficiently large, the maximum number of epochs should not affect the overall execution time of the methods. For COSMO, interrupt the optimization after 2000 epochs. For the non-linear version of DAGMA, we increased the maximum epochs to 7000. Overall, we interrupt the execution of an algorithm whenever hits wall time limit of 20000 seconds.
do not disable carly-stopping conditions the
As previously discussed in Subsection C.I
we perform a hyperparameter search on each model for each dataset. In particular, sampled the learning rate from the range (1e-4, 1e-2) and the
18

---

Published as a conference paper at ICLR 2024
regularization coefficients from the interval (le-4, le-1). For the specific constrained optimization parameters, such as the number of problems or decay factors, we replicated the baseline parameters, for which point the reader the original papers our implementation. For COSMO, sample hyperparameters from the ranges Table given our theoretical findings on the relation between acyclicity and temperature (Theorem 2. ensure sufficiently small acyclicity values. In the non-linear variant, we employ Multilayer Perceptrons with h = 10 hidden units for each variable.
D DETAILED COMPARISON WITH RELATED WORKS
D.1 COMPARISON WITH ENCO
Lippe et al. 2022) propose to learn a directed acyelic graph by jointly learning the probability of an are being present and of the are direction. The overall method, named ENCO, defines the probability of direct edge X, -> X, as
W uv (PP))))
(39)
where H Raxd and P Raxd are free parameters and is the sigmoid function. Although similar to our formulation, ENCO defines arc orientations as matrix that might be cyclic. In fact, the matrix P does ensure the transitivity property that an orientation matrix grants instead. However, the authors proved that, in the limit of the number of samples from the interventional distribution, ENCO will converge to a directed acyclic graph. comparison, our formulation always converges to directed acyclic graph and is thus adapt to perform structure learning in the observational context.
D.2 COMPARISON WITH NOCURL
Motivated by Hodge theory (Hodge 1989) NOCURL (Yu et al. 2021) proposes to parameterize DAG function of d-dimensional vector and directed possibly eyclic graph. By adopting notation, introduced in Section 4 we could report their decomposition as
W uv H ReLU(P, -P,);
(40)
for an arbitrary weight from node X, to node X.. Compared to our defimition, the use ReLU does not correspond to the approximation of an orientation matrix. In fact, the distance between the priorities directly affects the weight. Instead, by employing the shifted-tempered sigmoid in the definition of smooth acyelic orientation, in COSMO the priorities only determine whether an are is present between two variables. Further, NOCURL requires preliminary solution from which to extract the topological ordering of the variables. In turn, such preliminary solution requires the use of an acyclicity constraint in multiple optimization problems. Therefore, in practice, NOCURL does not learn the variables ordering in an unconstrained way and adjusts adjacency weights in the last optimization problem.
E ADDITIONAL RESULTS
In this section, we report further results on simulated DAGs with different noise terms, graph types, and inereasing numbers of nodes. For each algorithm. we present the mean and standard deviation of each metric on five independent runs. We report the Area under the ROC Curve (AUC), the True Pos- itive Ratio (TPR), and the Structural Hamming Distance normalized by the number of nodes (NHD). The reported duration of NOCURL includes the time to retrieve the necessary preliminary solution through two optimization problems regularized with the NOTEARS acyclicity constraint. We denote as NOCURL-U the variation of NOCURL that solves unique unconstrained optimization problem without preliminary solution. When not immediate, we highlight in bold the best result and in italic bold the second best result. We do not report methods exceeding wall time limit of 20000 seconds. For completeness, we report and discuss results on variance normalized datasets (Section E.18) and further comparisons with GOLEM (Section E.17) and DAP (Section F.19
19

---

Published as a conference paper at ICLR 2024
E.1 ER4 GAUSSIAN NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
0.867 = 1.01
0.953 = 0.04
0.984 + 0.02
88 =
DAGMA
0.707 + 0.57
0.940 0.04
0.985 + 0.01
781 + 193
30
NOCURL
1.653 = 0.17
0.942 +
0.967 + 0.01
822 + 15
NOCURL-U
5.623 = 0.92
+ 0.492
0.694 0.06
227 -
NOTEARS
0.913 = 0.60
1.00
0.973 0.02
5193 170
COSMO
1.388 0.69
0.917 + 0.04
0.961
99 2
DAGMA
1.026 +: 0.40
0.876 0.02
0.982
661 + 142
100
NOCURL
5,226 1.34
0.921 = 0.02
0.962 0.01
1664 15
NOCURL-U
10.108 = 4.11
0.427 0.05
0.682
267 10
NOTEARS
2.380 + 2.10
0.898 + 0.03
0.963 + 0.01
11001 340
COSMO
4.149 + 1.14
0.819 + 0.02
0.933 + 0.01
437 81
500
DAGMA
2.246 + 0.40
0.882 + 0.01
0.980 + 0.00
2485 366
NOCURL-U
27.675 + 16.52
0.410 + 0.04
0.683 + 0.05
1546 - 304
E.2 ER4 EXPONENTIAL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
СОSMO
0.600 + 0.54
0.970 = 0.02
0.989 + 0.01
89 - 3
DAGMA
0.613 0.91
0.958
0.05
0.986 = 0.02
744 75
30
NOCURT.
2.300 + 0.97
0.918
0.04
0.956 0.02
826 + 24
NOCURL-U
5.313 0.17
0.423 + 0.05
+ 0.694 0.05
212 + 5
NOTEARS
1.320 + 0.67
0.880
0.10
0.966 + 0.03
5579 + 284
СОSMO
1.642 0.26
0.952 0.02
0.985 0.01
99 + 2
DAGMA
1.294 + 0.52
0.944 = 0.02
0.986 + 0.01
733 + 109
100
NOCURL
5.652 1.35
0.854 0.03
0.00
1655 28
NOCURL-U
11.642 4.34
0.478 0.05
0.05 0.05
242 = 4
NOTEARS
1.156 0.44
0.904 = 0.03
0.91
10880 + 366
COSMO
2.342 +: 0.86
0.944 + 0.02
0.986
390 - 102
500
DAGMA
2.147 1.08
0.902 0.04
0.984
2575 469
NOCURL-U
20.183 + 7.43
0.437 = 0.03
0.715 0.03
1488 + 249
E.3 ER4 GUMBEL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
СОSMО
2.220 + 1.65
0.862 1 0.14
0.914
0.10
87 + 2
DAGMA
1.680
0.73
0.937 = 0.03
0.973 0.02
787 86
30
NOCURL
3,873 = 1.26
0.853 0.08
0.915 0.04
826 + 17
NOCURL-U
5.260 + 0.57
0.475 + 0.08
0.678
0.05
212 + 5
NOTEARS
0.587 0.38
0.962 0.03
0.981 + 0.01
5229 + 338
CОSMO
2.398 0.70
0.936
0.02
0.973
0.01
98 +
DAGMA
1.132 + 0.79
0.921
0.04
0.986 + 0.01
858 + 101
100
NOCURL
4.714 1.77
0,905
0.03
0.962
1675 34
NOCURL-U
6.914 + 0.80
0.383 = 0.04
0.663
247 + 9
NOTEARS
1.402 0.40
0.869 - 0.04
0.90
11889 + 343
COSMO
3.574 +: 1.44
0.932 + 0.02
0.982 0.01
410 + 106
500
DAGMA
1.737 0.64
0.871 0.03
0.980 - 0.00
2853 + 218
NOCURL-U
18.182 + 9.28
0.462 + 0.06
0.728 0.05
1342 + 209
20

---

Published as a conference paper at ICLR 2024
E.4 SF4 GAUSSIAN NOISE
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
0.300 + 0.09
0.973 = 0.01
0.997 + 0.00
89 =
DAGMA
0.360 + 0.30
0.973 0.02
0.996 0.01
653 + 198
30
NOCURL
0.967 + 0.43
0.893 +
0.983 + 0.01
828 + 23
NOCURL-U
4.410 = 0.72
0.566 0.11
0.741 = 0.08
226 +
NOTEARS
0.553 + 0.54
0.944 0.06
0.984 + 0.02
5292 + 261
CОSMO
0.482 0.31
0.962 0.02
0.991 0.01
DAGMA
0.712 + 0.33
0.951 + 0.02
0.995 + 0.00
479 + 75
100
NOCURL
2,030 + 0.46
0.883 0.03
0.982 0.01
1667 + 25
NOCURL-U
5.521 + 0.61
0.596 = 0.09
0.788 = 0.06
269 9
NOTEARS
0.280 0.35
0.972 0.04
0.993 0.01
10112 + 492
COSMO
1.566 + 0.68
0.953 + 0.02
0.989 + 0.01
541 + 15
500
DAGMA
1.343 0.46
0.915 0.04
0.992 0.00
1345 33
NOCURL-U
7.146 + 3.19
0.504 + 0.08
0.780 + 0.07
1394 = 217
E.5 SF4 EXPONENTIAL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
СОSMO
0,613 + 0.39
0.965 = 0.02
0.985 + 0.02
87 - 2
DAGMA
0.127 0.20
0.991
0.01
0.999 = 0.00
592 - 200
30
NOCURT.
0.887 + 0.21
0.845
0.02
0.985 + 0.01
824 + 25
NOCURL-U
4.067 0.73
0.460
0.15
0.685 + 0.09
212 + 7
NOTEARS
0.513 + 0.30
0.962 0.03
0.984 +
5189 + 271
CОSMO
0.724 0.71
0.963
0.04
0.985 0.02
100 + 2
DAGMA
0.586 + 0.56
0.969 = 0.03
0.995 + 0.00
395 + 108
100
NOCURL
1,998 0.40
0.907 0.03
0.980
1670 28
NOCURL-U
5,912 1.54
0.575 = 0.06
0.00
245 -
NOTEARS
0.910 0.43
0.962 = 0.02
0.991
10243 + 723
COSMO
1.445 + 0.58
0.950 - 0.03
0.990 - 0.01
517 + 108
500
DAGMA
1.653 0.91
0.873 0.08
0.988 0.01
1466 1 247
NOCURL-U
12.140 + 7.84
0,482 + 0.08
0.727 + 0.06
1205 + 257
E.6 SF4 GUMBEL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
СОSMО
0.467 + 0.51
0.962 + 0.05
0.990 + 0.02
88 + 2
DAGMA
0.487 = 0.20
0.956
0.03
0.990
0.01
754 = 179
30
NOCURL
0.747 1 0.19
0.938 1 0.02
0.989 0.00
826 + 32
NOCURL-U
3.107 + 0.64
0.460 0.06
0.737 + 0.04
213 +: 5
NOTEARS
0.860 + 0.76
0.06 0.924
0.975 0.02
5199 + 130
CОSMO
0.864 0.24
0.968
0.01
0.992 0.01
98 + 2
DAGMA
0.388 + 0.30
0.975 0.02
0.997 + 0.00
422 + 103
100
NOCURL
1.806 + 0.40
0.898 0.03
0.982 0.01
1676 31
NOCURL-U
8.756 + 2.65
0.550 0.05
0.757 + 0.03
245 +
NOTEARS
1.134 + 0.81
0.894 + 0.08
0.989 0.01
11618 + 1309
COSMO
1.426 + 0.53
0.951 + 0.03
0.994 + 0.00
524 + 22
500
DAGMA
1.384 0.38
0.849 0.04
0.991 0.00
1359 34
NOCURL-U
8,931 + 7.05
0.430 0.10
0.741 + 0.08
1193 + 229
21

---

Published as a conference paper at ICLR 2024
E.7 ER6 GAUSSIAN NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
CОSMO
4.087 + 1.12
0.838 + 0.06
0.921
0.04
89 +
DAGMA
2.367 = 0.63
1047
0.958 0.01
665 249
30
NOCURL
4.480 + 0.92
0.869 + 0.03
0.908 = 0.03
909 + 18
NOCURL-U
7.490 + 1.18
0.459 + 0.08
0.672 0.06
226 +
NOTEARS
3.327 + 1.65
0.840 + 0.07
0.922 0.04
5239 427
COSMO
9.476 3.01
0.771 + 0.08
0.911 0.05
98 18:22
DAGMA
10.740 + 2.83
0.709 + 0.13
0.902 + 0.04
761 + 134
100
NOCURL
15.044 + 1.60
= 0.04
0.888 = 0.02
1687 = 26
NOCURL-U
30.719 5.20
0.43
0.580 0.04
268 9
NOTEARS
6.556 + 3.10
0.842 0.05
0.944 0.02
12053 940
CОSMO
25.443 + 4.47
0.736 + 0.01
0.937 = 0.01
526 + 100
500
DAGMA
15.952 + 1.67
0.553 + 0.05
0.925 0.01
3207 271
NOCURL-U
165.465 + 20.86
0.433 + 0.02
0.558 + 0.03
1226 + 293
E.8 ER6 EXPONENTIAL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
3.300 +: 0,95
0.897 + 0.05
0.947 + 0.03
89 + 2
DAGMA
3.480 = 1.42
0.861 0.06
0.945 + 0.03
672 + 177
30
NOCURL
4.573 = 0.78
0.846
0.902 0.03
897 = 13
NOCURL-U
8.700 + 0.89
0.426
0.615 0.06
226 9
NOTEARS
2.313 + 1.55
0.881 + 0.09
0.953 + 0.04
5516 + 652
COSMO
10.170 2.74
0.768 + 0.09
0.919 + 0.04
99 +
DAGMA
8.118 - 3.10
0.793 + 0.11
0.934 0.04
681 + 149
100
NOCURL
14,860
4.67
0.685
+ 0.10
0.863 + 0.06
1735 39
NOCURL-U
30.600 + 4.34
0,450 0.591 0.04
267 + 8
NOTEARS
5.208 = 2.54
0.796 0.09
0.943 0.03
12663 + 1555
COSMO
25.854 + 4.28
0.741 + 0.04
0.943 0.01
460 + 123
500
DAGMA
16.417 - 4.45
0.571
- 0.11
0.925 1 0.02
4069 + 580
NOCURL-U
152.336 + 31.97
0.425 + 0.02
0.567 + 0.03
1363 + 306
E.9
ER6 GUMBEL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
2.840 + 1.08
0.906 + 0.04
0.954 + 0.03
89 +
DAGMA
2.727 + 0.83
0.906 + 0.02
0.964 0.02
634 194
30
NOCURL
5.003 = 0.72
0.811 +.01 0.04
0.891 0.03
902 + 9
NOCURL-U
8.153 + 0.96
0,422 +: 0.07
0.629 + 0.04
226 +
NOTEARS
2.740 + 1.61
0.791 + 0.10
0.938 + 0.04
5416 446
COSMO
10.048 3.155
0.780 + 0.07
0.899 0.06
100 +
DAGMA
7.910 + 3.05
0.805 + 0.09
0.935 + 0.04
715 + 203
100
NOCURL
11.932 + 2.68
0.742 0.04
0.894 0.03
1688 + 34
NOCURL-U
27.401 + 4.42
0.431+0.00 0.05
0.600 1 0.04
266 +
NOTEARS
4.884 + 0.45
0.833 + 0.05
0.951 = 0.01
12634 + 639
COSMO
26.148 + 4.86
0.740 + 0.04
0.941 + 0.02
418 + 106
500
DAGMA
16.358 = 4.94
0.563 + 0.07
0.921 0.02
3527 241
NOCURL-U
125.858 + 36.61
0.367 + 0.06
0.571 + 0.02
1612 + 27
22

---

Published as a conference paper at ICLR 2024
E.10 SF6 GAUSSIAN NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
1.273 = 1.07
0.907 0.10
0.963 + 0.06
89 =
DAGMA
1.107 + 0.37
0.930
0.985 + 0.01
456 39
30
NOCURL
1.573 = 0.46
0 0.80
823 + 14
NOCURL-U
4.997 = 0.98
0.05
0.732 0.05
226 -
NOTEARS
0.933 + 0.71
0.919 1.05
0.984 + 0.02
5313 184
COSMO
4.478 2.88
0.776
0.874 + 0.11
99 +
DAGMA
2.024 + 0.71
0.914 + 0.02
0.987 + 0.00
396 53
100
NOCURL
2.824 = 0.39
0.980 0.00
1679 27
NOCURL-U
10.556 = 6.00
0.542 0.07
0.751 = 0.08
266 5
NOTEARS
1,412 + 0.59
0.939 + 0.03
0.990 + 0.01
11156 170
COSMO
4.670 + 1.99
0.912 + 0.02
0.984 + 0.00
460 + 70
500
DAGMA
3.825 0.19
0.746 + 0.03
0.982 + 0.00
1418 54
NOCURL-U
19.793 + 11.03
0.368 + 0.04
0.698 + 0.04
1137 + 231
E.11 SF6 EXPONENTIAL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
CОSMO
1.393 + 1.24
0.926
0.05
0.975 + 0.03
88 + 1
DAGMA
1.147
0.48
0.943
0.03
0.982 + 0.01
578 = 173
30
NOCURT.
1.987
0.54
0.757
0.08
0.967 0.01
820 + 8
NOCURL-U
4.787
0.99
0.534
0.07
0.761 + 0.06
227 + 7
NOTEARS
0.49 0.753
0.943 0.04
0.986 + 0.01
5312 + 258
CОSMO
3.836 2.75
0.864
0.09
0.944 + 0.05
98 + 2
DAGMA
1.532 + 0.61
0.887
0.04
0.988 + 0.00
373 + 88
100
NOCURL
2.890 0.61
0.910
0.02
0.977 + 0.00
1692 28
NOCURL-U
6.607 1.05
0.474
0.06
0.06
266 -
NOTEARS
1.784 0.52
0.939 = 0.02
0.988
11369 - 519
COSMO
3.144 + 0.47
0.919
0.02
0.989 - 0.00
457 + 81
500
DAGMA
3.854 0.34
0.750 0.01
0.977 0.01
1384 33
NOCURL-U
13.763 + 8.79
0.389 0.05
0.728 + 0.06
1436 + 230
E.12 SF6 GUMBEL NOISE
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
1.047 + 0.42
0.938 + 0.03
0.984 + 0.01
88 +
DAGMA
1.347 = 0.63
0.933 - 0.02
0.981 + 0.01
528 + 67
30
NOCURL
1.787 = 0.52
0.898 0.02
0.969 = 0.01
822 + 29
NOCURL-U
5.577 + 0..43
0.549 0.06
0.733 + 0.05
225 +
NOTEARS
1.053 0.59
0.911 0.04
0.978 + 0.02
5429 251
COSMO
3.486
2.62
0.879 + 0.10
0.947 0.06
99
DAGMA
1.418 + 0.34
0.910 + 0.03
0.990 + 0.00
424 + 90
100
NOCURL
3.074 + 0.50
0.893 + 0.02
0.976 0.00
1682 - 22
NOCURL-U
9,643 + 4.59
0.464 0.08
0.712 0.10
267 -
NOTEARS
1.586 + 1.39
0.913 + 0.06
0.987 + 0.01
11820 985
COSMO
3.288 + 0.50
0.931 + 0.01
0.992 + 0.00
429 + 87
500
DAGMA
4.055 0.88
0.802
- 0.03
0.981 + 0.00
1465 138
NOCURL-U
56.103 + 41.06
0.420 + 0.06
0.648 + 0.07
1201 + 253
23

---

Published as a conference paper at ICLR 2024
E.13
ER4 NON-LINEAR MLP
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
2.530 + 0.45
0.752
0.06
0.923
0.03
154 + 2
20
DAGMA-MLP
2,853 + 0.37
0.810 0.15
0.932 +9922+00 0.04
2252 + 40
NOTEARS-MLP
3.270 = 0.52
0.904
0.07
0.956 0.02
3554 = 154
COSMO
3.295 = 0.45
0.712 + 0.05
0.925 + 0.02
177 + 2
DAGMA-MLP
3.633 = 0.81
0.804 + 0.10
0.942 0.03
2622 + 28
CОSMO
5.164 + 1.38
0.582 + 0.04
0.895 + 0.02
311 10
100
DAGMA-MLP
2.337 + 0.18
0.506 + 0.05
0.912 + 0.01
4866 + 69
E.14
SF4 NON-LINEAR MIP
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
1.525 + 0.18
0.743 + 0.05
155 +
20
DAGMA-MLP
1,408 = 0.26
0.746 0.11
2270 + 27
NOTEARS-MLP
1.040 = 0.33
0.922 0.08
0.00
3400 184
COSMO
2.271 = 0.24
0.583 = 0.05
0.95
174 +
40
DAGMA-MLP
1,875 = 0.34
0.744
0.972
2588 + 45
CОSMO
3.144 = 0.23
0.455 = 0.04
313 + 6
100
DAGMA-MLP
3.051 = 0.14
0.260 + 0.04
0.961 0.01
4883 + 77
E.15
ER6 NON-LINEAR MLP
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
2.995 = 0.45
0.667 = 0.07
0.919
155 + 2
20
DAGMA-MLP
2.992 + 0.39
0.712 + 0.08
0.917 0.02
2252 + 35
NOTEARS-MLP
2.895 + 0.47
0.862 + 0.08
0.949 0.02
3557 + 130
CОSMO
4.837 = 0.66
0.598 + 0.07
0.02 0.891
174 +
40
DAGMA-MLP
3.732 + 0.71
0.633 + 0.13
0.919 + 0.02
2570 + 44
COSMO
6.049 + 0.91
0.501 + 0.04
0.875 + 0.02
308 + 4
100
DAGMA-MLP
3.790 + 0.24
0.535 + 0.05
0.908 + 0.01
4807 + 80
16
SF6 NON-LINEAR MIP
d
Algorithm
NHD
TPR
AUC
Time (s)
CОSMO
1,855 + 0.28
0.748 + 0.06
0.95
154 + 2
20
DAGMA-MLP
1.833 + 0.52
0.766+0.113
0.03
2249 + 40
NOTEARS-MLP
1.355 = 0.39
0.926 0.08
0.980
3512 = 193
CОSMO
3.394 = 0.25
0.501 + 0.06
0.946
174 +
40
DAGMA-MLP
2.720 = 0.38
0.724 + 0.08
0.967 0.01
2572 - 100
COSMO
4.744 = 0.21
0.342 0.04
0.931 0.01
308 +
100
DAGMA- -LPP
4.149 = 0.17
0.391 + 0.06
0.968 + 0.01
4699 = 338
24

---

Published as a conference paper at ICLR 2024
Graph
Algorithm
NHD
TPR
AUC
h(W)
Time (s)
COSMO
0.867 + 1.01
0.953 + 0.04
0.984 0.02
0.000 0.00
88 3
30
GOLEM
0.833 + 1.14
0.933 + 0.09
0.975 + 0.04
0.025 + 0.03
305 + 4
ER4
COSMO
1.388 + 0.69
0.917 = 0.04
0.961 0.03
0.000 + 0.00
99 2
100
GOLEM
4.874 + 1.02
0.802 + 0.04
0.936 + 0.01
0.548 + 0.22
660 = 11
COSMO
4.087 + 1.12
0.838 = 0.06
0.921 + 0.04
0.000 = 0.00
89 4
30
GOLEM
6,153 + 0.56
0.496 = 0.13
0.749 0.10
0.382 = 0.14
306 + 9
ER6
COSMO
9.476 + 3.01
0.771 = 0.08
0.911 = 0.05
0.000 + 0.00
98 =
100
GOLEM
9.572 +: 2.32
0.423 + 0.08
0.794 + 0.07
4.654 + 1.34
664 +
Table 5: Experimental comparison with GOLEM (gg al. ab.
E.17 GOLEM NG ET AL..
GOLEM (Ng et al. is unconstrained approach for continuous structure learning that employs the NOTEARS acylicity constraint as a regularization term. Given this strategy, GOLEM does not ensure DAG convergence and instead requires final pruning step to iteratively remove arcs until the resulting adjacency matrix is acyclic. As stated in the Related Works (Section 2). we focus our main empirical comparison methods ensuring the acyclicity of the solution. Therefore, in the following table, we only briefly report an empirical analysis of GOLEM. Notably, on small and sparse graphs (d 30, ER4), the resulting acyclicity might be tolerable, but the value is significant on both larger (d = 100) and denser (ER6) ones, We also report COSMO results on the same graph types to recall that our proposal achieves acyclicity zero by construction (Table
E.18 NORMALIZED VARIANCE RESULTS
As we remarked in the main body, the choice of optimizing an adjacency matrix by minimizing the MSE is known to exploit the variance of the variables (Loh & Bühlmann, 2014) and might not be significant outside of this assumption (Reisach et al., [021] Kaiser & Sipos [2022] Ng et b). In general, variance normalization is enough to induce a significant performance drop approaches adopting this loss function. However, it is worth remarking that structure learning formulations derived from NOTEARS, like ours, are loss-agnostic and have been used in several continuous causal discovery methods with different but appropriate loss functions (Brouillard al. 2020 Lorch et al. 2022. Hägele et al. 2023). Defining function classes, appropriate losses, and identifiailitty results are important requirements for causal discovery solutions, but are out of the scope of the current work. However, for clarity, we show that, when fitting a DAG through the MSE loss, COSMO is equally affected by variance normalization as competing approaches (Table
(b.
Reisach et al. 2021)) proposes a simple baseline named SORTNREGRESS to evaluate whether the variables order follow their variance in non-normalized datasets. The algorithm consists of two steps: (i.) sorting the nodes according to their variance and (ii.) regressing variables on their predecessors. In datasets where there is strict relation between the topological ordering and the nodes' variance, the SORTNREGRESS baseline is expected to perform consistently well. To better contextualize our work, we show that in the simulated test bed adopted by NOTEARS, DAGMA, NOCURL, and COSMO, the baseline performs well as expected (Table 7). Notably, it is evident how, for larger graphs, our parameterization is still faster the SORTNREGRESS baseline, despite the simplicity of the latter.
E.19 DAGUERREOTYPE ZANTEDESCHI ET AL. 2022)
DAGUERREOTYPE (Zantedeschi et al. 2022) proposes an end-to-end differentiable model for causal DAG discovery that parameterize the underlying causal ordering over the polytope of permutation vectors. In practice, the model fits a distribution over possible permutations parameterized by score vector analogous to the priority vector p used by COSMO. The differentiable transformation from the score vector to the permutation can be realized by either the SparseMAP (Niculae et al.) 2018) the Top-k sparsemax (Correia et al.) al.] (020) operators. Further, the scores and the functional relations adjacencies in the linear case
be either learned jointly with bi-level optimization scheme.
25

---

Published as a conference paper at ICLR 2024
Graph
Algorithm
NHD
TPR
AUC
k(W)
CОSMO
4.637 + 0.43
0.095 + 0.04
0.440 + 0.06
0.000 + 0.00
DAGMA
4.967 0.41
0.102 0.04
0.475
0.08
0.000 0.00
ER4
30
NOCURL
4.910 = 0.53
1.09==000
0.429 = 0.07
0.000 0.00
NOCURL-U
4.950 0.31
0.092 0.03
0.530 0.05
0.000 1.00
NOTEARS
4.920 + 0.36
0.098 + 0.04
0.473 + 0.06
0.000 +.000
COSMO
4.150 0.05
0.036 0.01
0.533 + 0.02
0.000 + 0.00
DAGMA
4.205 = 0.10
0.120 0.02
0.571 = 0.02
0.000 0.00
ER4
100
NOCURL
4.412 0.09
0.054 0.01
0.501 = 0.03
0.000 0.00
NOCURL-U
4.720 = 0.10
0.038 0.02
0.513 + 0.03
0.000 0.00
NOTEARS
4.221 0.10
0.121 0.02
0.613 0.02
0.000 0.00
COSMO
6.643 + 0.29
0.039 0.02
0.420 0.05
+.00
DAGMA
6,837 = 0.34
0.062 - 0.02
0.465 = 0.04
0.000 0.00
ER6
30
NOCURL
6.820 + 0.26
0.032 0.02
0.368 + 0.06
0.000 0.00
NOCURL-U
6.947 0.25
0.055 + 0.02
0.492 0.02
0.000
NOTEARS
6.917 = 0.37
0.056 + 0.03
0.497 = 0.04
0.000
COSMO
6.130 + 0.06
0.010 0.00
0.513 + 0.02
0.000 0.00
DAGMA
6.570 0.09
0.046 0.01
0.499 0.01
0.000 0.00
ER6
100
NOCURL
6.505 - 0.11
0.018
0.428 = 0.04
0.000
NOCURL-U
6.884 0.09
0.020 0.01
0.504 0.03
0.000 0.00
NOTEARS
6.533 = 0.14
0.049 0.01
0.532 = 0.02
0.000 0.00
Table 6: Experimental comparison on variance-normalized datasets.
Graph
d
Algorithm
NHD
TPR
AUC
Time (s)
COSMO
0.867
1.01
0.953
0.04
0.984 + 0.02
88 +
ER4
30
SORTNREGRESS
1.973 + 0.76
0.890 + 0.05
0.938
0.03
11 +
COSMO
1.388 + 0.69
0.917 + 0.04
0.961 + 0.03
99 +
ER4
100
SORTNREGRESS
1.747
0.85
0.937 0.03
0.972
0.01
98
COSMO
4.149 + 1.14
0.819 = 0.02
0.933 0.01
437 +
ER4
500
SORTNREGRESS
1.557
0.39
0.942 0.02
0.978 0.01
3825 171
COSMO
0.300 0.09
0.973 0.01
0.997 0.00
89
SF4
30
SORTNREGRESS
0.517 + 0.39
0.966 + 0.03
0.987 + 0.01
15 +
COSMO
0.482 0.31
0.962 0.02
0.991 0.01
99 + 3
SF4
100
SORTNREGRESS
0.765 + 0.33
0.958 + 0.02
0.989
0.01
99 +
COSMO
1.566 + 0.68
0.953 + 0.02
0.989 + 0.01
541 + 15
SF4
500
SORTNREGRESS
0.788 + 0.36
0.955 0.03
0.986 0.01
5373
101
Table 7: Experimental comparison with SORTNREGRESS
(Reisach al. 2021
26

---

Published as a conference paper at ICLR 2024
Graph
Algorithm
NHD
TPR
Time (s)
COSMO
0.867 + 1.01
0.953 0.04
88
DGT-SPMAP-JOINT
5.793 = 0.75
0.282
0.12
91 +
ER4
30
DGT-SPMAX-JOINT
6,353 0.80
0.398 0.09
140 +: 15
DGT-SPMAP
6.087 0.97
0.290 0.14
1829 896
DGT-SPMAX
5.727 = 0.82
0.198 0.10
2006 + 559
COSMO
1.388 = 0.69
0.917 = 0.04
99 =
DGT-SPMAP-JOINT
9.460 + 0.81
0.190 + 0.04
847 +: 297
ER4
100
DGT-SPMAX-JOINT
8.724 1.92
0.187 0.10
1261 429
DGT-SPMAP
8.506 + 0.91
0.149 + 0.07
20200
12514
DGT-SPMAX
8.152 0.75
0.128 0.05
24554
13113
COSMO
0.300 = 0.09
0.973 + 0.01
89 5
DGT-SPMAP-JOINT
3.813 = 0.35
0.382 + 0.08
72 + 17
SF4
30
DGT-SPMAX-JOINT
3.787 0.80
0.624 0.14
157 +
DGT-SPMAP
3.820 0.39
0.365 = 0.11
2246 + 701
DGT-SPMAX
3.740 + 0.24
0.333 0.12
2825 + 1218
COSMO
0,482 = 0.31
0.962 + 0.02
99 - 3
DGT-SPMAP-JOINT
4.686 0.73
0.317 0.16
946 249
SF4
100
DGT-SPMAX-JOINT
4.662 + 0.68
0.351 + 0.19
935 + 296
DGT-SPMAP
4.686 = 0.77
0.293 0.19
12688 11499
DGT-SPMAX
4.268 = 0.43
0.137 0.04
9812 + 5257
Table 8: Experimental comparison with DAGUERREOTYPE
Zantedeschi al. 2022
As further comparison, present empirical results testbed of DAGUERREOTYPE in its four main configurations: joint optimization with SparseMAP (DGT-SPMAP-JOINT) joint optimization with Top-k (DGT-SPMAX-IOINT), bi-level optimization with SparseMAP (DGT-SPMAP), and bi-level optimization with Top-k (DGT-SPMAX). In all experiments runs we assume Gaussian noise terms, employ default hyperparameters from the original implementationf] and fix the sparse operator parameter K = 10, as reported for the synthetic evaluation in Zantedeschi et al. 2022 Appendix D).
As highlighted by the original authors, joint optimization is significantly faster than bi-level opti- mization. which in turn takes significantly more time than COSMO on larger graphs (Table B) On all analyzed configurations, the result highlights the scalability issues of DAGUERREOTYPE when compared to simpler optimization schemes. Further, the empirical results confirm the expectation terms computational complexity; COSMO has quadratic complexity on the number of nodes, while DGT-SPMAP has complexity O(Kd²) and DGT-SPMAX has complexity O d(?²d)). Further, the higher variance of the computational time results of DAGUERREOT might be explained by the use of early stopping model convergence. As expected, bi-level optimization generally improves performance DAGUERREOTYES but we still notice a non-ne negligiiee performance degradation compared to sparser ER2 and SF2 graphs reported in Zantedeschi et al. (2022).
L2 REGULARIZATION
To ensure sparsity of the solution, most continuous structure learning approaches apply LI regulariza- tion to the weighted adjacency matrix (Ng et al. b). As reported Section 4 we adopt the same approach only on the direct matrix H avoid influencing the priority vector p. Apart from this, found beneficial to also apply L.2 regularization to the direct matrix H. Intuitively, this avoids that the directed matrix "constrast" small values of the smooth orientation matrix given by ares opposed the acyclicity during training against
to the partial order. In Figure 55] we offer a visualization the number of epochs. In particular. non-regularized matrices tend to be more acyclic (Figure 5a) and reach acyclicity later in the optimization process at parity of temperature (Figure 5b). In prac-
27

---

Published as a conference paper at ICLR 2024
1.73
14-
0.03
152
0.0
0.83
05
25
23
01
93
0.5
0.5
20
2.0
1.00
3.0
0.75
0.53
$20
1000
1750
1702
1900
1950
2002
Epoch
(a) All Epochs
(b) Last Epochs
Figure 5: We plot the acyclicity of the learned adjacency matrix W against the temperature / for different regularization terms A2 on the directed adjacency matrix H. Results are obtained over independent iterations of COSMO on an ER4 graph with Gaussian noise. In Subfigure (a) we plot the whole optimization process, while we zoom the last epochs in Subligure (b).
3.2
24-
3.5
0.5
3.5
250
756
1259
3508
2750
250
706
1259
1750
(a) Temperature
(b) Acyclicity
3.30
3.5
0.2
24
66
1007 1292
1928 1792
2003
240
14.
66
1993
1952
1988
1792 2002
NHD
(d) AUC
Figure 6: We report the values of (a) temperature, (b) acyclicity, (c) normalized Hamming distance, (d) and ROC-AUC against each optimization epoch of run COSMO on 1000 samples from an ER4 graph with 30 nodes and Gaussian noise. Notably, the model reaches an acyclic solution (As Epoch 1700) and keeps improving before the temperature reaches zero the last epoch.
tice, our hyperparameter search consistently awarded non-negligible L2 regularization coefficients (Qg Rs 2e-2).
28

---

Published as a conference paper at ICLR 2024
G TEMPERATURE ANNEALING AND MODEL OPTIMIZATION
In this section, we highlight the interaction between the temperature value and the acyclicity of the solution, together with the performance metries on graph reconstruction. We consistently observe that the solution passes across two distinct phases (Figure 6). In all the reported metrics, after an initial convergence, the model settles with a preliminary solution. Then, after the temperature hits sufficiently low value the acyclicity of the model rapidly decreases and the performance increase. In partieular, it is worth noticing that the model reaches an acyclic solution before hitting the minimum value of the temperature and that performance metrics still vary at acyclicity zero.
In Figure 7 we propose a more detailed visualization of the evolution of the parameters during the training of COSMO, where the same phenomenom can be observed. In fact, after a first phase (Epochs 1-500) where the parameters substantially vary, the optimization is substantially stable (Epochs 500-1500) until sufficient decrease in temperature (Epochs 1500-1700). As before, since the model converges to an acyclic solution before the temperature reaching zero (Epochs 1800-1999), we can notice that optimization continues even in this final phase, mostly increasing sparsity of the graph.
29

---

Published as a conference paper at ICLR 2024
Orientation (p)
Directed H
Weighted W
Binary [W] > 0.3
Epoch
Epoch 500
Epoch 1500
Epoch 1700
Epoch 1800
Epoch 1999
Ground Truth
Figure 7: We report the adjacency matrices learned by COSMO at significant epochs of 1000 samples from an ER4 graph with 30 nodes and Gaussian noise. As in the main body, we indicate positive values in blue and negative values in red.
30


Reviews:
[{"review_id": "9S29TPA68g", "reviewer": "Reviewer_Rsrs", "summary": "This paper introduces a constraint-free continuous optimization scheme (COSMO) to the structure learning problem that focuses on acyclic graph reconstruction. COSMO employs a differentiable relaxation of the acyclic orientation matrix, allowing it to learn a directed acyclic graph without explicitly enforcing acyclicity constraints. The key contributions include the introduction of a smooth orientation matrix, the development of COSMO, which is significantly faster than existing constrained methods, and a thorough experimental comparison highlighting COSMO's superior performance in terms of structure recovery and computational efficiency. The paper provides valuable insights and tools for more efficient and unconstrained acyclic structure learning.", "questions": "See above \"weeknesses\" section.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- A primary strength of the method is its provision of an unconstrained continuous optimization approach for acyclic structure learning.\n- The method introduces a novel, differentiable approximation known as the \"smooth orientation matrix,\" which depends on a temperature parameter.\n-  The method outperforms certain constrained approaches in terms of speed due to the quadratic number of operations needed to reconstruct the DAG.", "weaknesses": "-  In Figure 1, are you learning a directed graph or a partially directed graph with undirected edges to be identified later?\nPlease formally define the priority vector in the notation section.\n- While Figure 1 appears to be an illustration of how COSMO works, it's not very clearly explained. It would be helpful to provide a proof sketch or illustration with an example in the introduction to briefly explain how your algorithm works (given this example) and why it's an unconstrained optimization method.\n- Can you provide information on the computational time and sample complexity of COSMO?\n- How do you guarantee the identifiability of the recovered graph, or do the results lead to some equivalence class?\n- In the related work section, consider covering some classic DAG structure learning methods.\n- The method relies on the temperature parameter for the sigmoid function. While annealing the temperature is shown to be effective, selecting an appropriate temperature schedule could be a challenging task, and the method's performance may be sensitive to the choice of temperature.\n- The effectiveness of this approach might depend on the initial conditions. Finding suitable initializations for the priority vector and other parameters could impact convergence and the quality of the learned DAG.\n- Figure 4 does not seem very necessary, and it's not very clear. Consider moving Figure 4 to the appendix and using the space to discuss the sensitivity analysis of your method's parameter choices.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZfzDuBJxCn", "reviewer": "Reviewer_Ar2D", "summary": "the paper proposes a constrain-free DAG learning formulation, based on existing notears and nocurl framework. The main innovation is to use a smoothed matrix to indicate the topological order, which enables the joint optimization of the smoothed topological matrix and a weighted adjacency matrix. Results show improved accuracy, esp when the number of dimension is high.", "questions": "see above", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- the paper proposes an improved DAG learning approach without constraint, which is novel \n\n- empirical performance on the large number of variables is solid", "weaknesses": "- Experiment: When the number of variable is small or medium, the accuracy seems to be worse than some existing approaches. In addition, some experiment results are missing and some elaboration is needed.\n\n- Presentation could be improved to be more scientifically sound.  \n\n\nDetails and questions:\n\n- Figure 3: since you showed the results with 1e4 nodes, what is the accuracy performance of these methods with large d?\n- page 13 in appendix, last line: how did the denominator become $(1 + e^{-\\epsilon})^n$ in the last step?\n- authors seem to over claim the contributions and often make some imprecise statements. The formulation is based on the existing nocurl formulation, with $W = H \\cdot f(p)$. Nocurl also learn a topological order with p along with the typical weighted adjacency matrix H. The difference is that nocurl learns topology from a preliminary solution, while the proposed COSMO optimizes them jointly with a smoothed formulation. This in itself is a worthy contribution, no need to make exaggeration statements such as claiming \"the first unconstrained optimization approach\" (which would be nocurl-u by authors' definition) and nocurl is \"a slightly similar model\". Authors should discuss the differences. \n- learning topology seems to be a hard task by itself. Gradient-based approach may not do too well on it, which explains the poorer results when the number of variable is smaller. What are authors observation on the change in topology over time (by the orientation matrix)? Does it have significant change in the beginning and then barely change toward the end of optimization steps? How does the temperature annealing change the topology over time?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "CMGl2zmjbl", "reviewer": "Reviewer_tdcC", "summary": "The paper studies the problem of DAG structure learning from a continuous score-based perspective. In order to optimize a given score over the space of DAGs, the authors propose to decompose the parameterization by a real matrix $H$ and an orientation matrix induced by a priority vector $p$, a threshold $\\epsilon$ and a temperature $t$. Then, it is shown that asymptotically as $t\\to 0$ the orientation matrix corresponds to a adjacency matrix of a DAG. The authors optimize $H$ and $p$ jointly and observe speed ups over SOTA methods.", "questions": "* Consider a $p$ vector with all entries being different. Then, for some arbitrarily small $\\epsilon$, it seems to me that the orientation matrix is the adjacency matrix of a fully-connected DAG. In some sense, $\\epsilon$ controls the sparsity level of the output, however, this is problematic because $\\epsilon$ is not something being optimized but instead fixed a-priori. Please correct me if I am wrong.\n\n* In the paragraph about \"Direct Matrix Regularization\". For the values found for $H$ and $p$ after optimizing, why not simply threshold the smooth orientation matrix by a small amount so that the \"small entries\" are set to zero and are not affected by the possibly large values of $H$. Given that linear models with equal variances are identifiable from the score function, L2 regularization of the weights should not be necessary and I wonder if doing this affects the performance of COSMO.\n\n* I think more details about the optimization procedure should go to the main text. In particular, the scheduling for $t$ seems very important  and it was not clear how it was set in the main text.\n\n* Related to my point above, I might have missed this since I only briefly looked at the appendix. But the NOTEARS $h$ is not a good proxy for acyclicity as discussed by Bello et al. (2022). In fact, given that your method does not rely on $h$, why not simply threshold the smooth orientation matrix by a very small amount and check in $O(|V|+|E|)$ that the adjacency corresponds to a DAG?\n\n* In the experiments, the runtime for DAGMA looks different to what was reported by Bello et al. (2022). These might be because of the available implementation of DAGMA is using different hyperparameters for the optimization, but I am curious what values were used in the comparison for the main text.\n\n* I should also note a related work on learning DAGs by Deng et al. (2023) \"Optimizing NOTEARS objectives via topological swaps\", which obtains better performances in structure estimation although somewhat slower.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is overall clear and makes a fair contribution to structure learning, which is a significant problem. \nThe paper is also of good quality, the authors presented a comprehensive set of experiments, showcasing the competitive performance of their method.", "weaknesses": "The main weakness of the paper is its novelty. The re-parametrization itself is not new, that is, the idea of decomposing a weighted adjacency matrix by a real matrix and orientation matrix. Comparing to NOCurl, the change is basically from ReLU to Sigmoid. However, I should give credits to the authors for making this arguably small change work well in practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ryy7V36UWj", "reviewer": "Reviewer_LvdG", "summary": "The paper proposes a new continuous relaxation of the acyclicity constraint for learning DAGs that is based on learning a priority score per node (that implies a partial ordering of the nodes) and relaxing its corresponding binary directed adjacency matrix to obtain non-zero gradients.\nSimilar to NoCurl (Yu et al. 2021), a node is an ancestor of another iff its priority score is smaller than the priority score of the other node minus a margin $\\epsilon$ (that allows to represent also partial orderings). This Heaviside function is then relaxed by means of a tempered sigmoid, allowing for optimization by gradient descent.", "questions": "(minor) In page 15, the wrong theorem is reported.", "limitations": "", "rating": 5, "confidence": 5, "soundness": 1, "presentation": 4, "contribution": 3, "strengths": "1. The paper tackles an important problem of interest to the ICLR community by proposing a novel method for learning DAGs from observational data.\n\n2. The overall complexity of the optimization is quadratic in the number of nodes (which is one of the lowest among DAG learning methods) allowing for scaling to graphs of $> 10^4$ nodes in the experimental analysis. The empirical analysis shows that this lower complexity translates in practice to competitive training times.\n\n3. The paper is overall well presented. In particular, the proofs of appendix A and B are sound and easy to follow.", "weaknesses": "In order of gravity:\n\n1. (major) The paper neglects the order-based literature for DAG learning which is closely related to this work. Only VP-DAG (Charpentier et al. 2021) is mentioned, but a considerable line of research exists, starting from the seminal works [1,2], to the most recent ones [e.g. 3,4] based on permutation matrices, and in particular [5] that also makes use of a priority vector. The paper does not highlight the pros and cons wrt this literature nor compare with any of the order-based methods empirically.\n\n2. (major) Experiments are carried out on synthetic data with weight settings that are known to have trivial solutions [Reisach et al. 2021]. Hence they are insufficient to validate the method. The experimental analysis should be carried out on more challenging settings (e.g., smaller edge weights) where the trivial baseline sortnregress [Reisach et al. 2021] is not state-of-the-art.\n\n3. By using the tempered sigmoid at training, the adjacency matrix is never acyclic. Acyclicity is only obtained in the limit of the temperature tending to 0, that however turns the gradients to 0. This implies that careful temperature scheduling is required in order to converge to a good solution.\n\n4. The method relies on edge thresholding at test time to remove spurious relationships. However a method for choosing such threshold is not presented.\n\n[1] Nir Friedman and Daphne Koller. Being bayesian about network structure. A bayesian approach to structure discovery in bayesian networks. Machine learning, 50, 2003\n\n[2] Marc Teyssier and Daphne Koller. Ordering-based search: A simple and effective algorithm for learning bayesian networks. In UAI, pp. 548–549. AUAI Press, 2005. \n\n[3] Chris Cundy, Aditya Grover, and Stefano Ermon. Bcd nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34, 2021.\n\n[4] Ming Gao, Yi Ding, and Bryon Aragam. A polynomial-time algorithm for learning nonparametric causal graphs. In Advances in Neural Information Processing Systems, 2020\n\n[5] Zantedeschi, V., Franceschi, L., Kaddour, J., Kusner, M. J., & Niculae, V. (2023). DAG Learning on the Permutahedron. In ICLR, 2023\n\n### Update after discussion\nI increased my score to 5, to acknowledge the improvements made to the paper during the discussion. Still I cannot recommend acceptance as is for the following reasons:\n1. The results in the main text are still exclusively on possibly flawed settings. Because of this, these results show that the proposed method is much faster than the considered baselines, but they are not conclusive regarding its ability to perform structure learning. The authors added results on normalized setting in the appendix, however without reporting running times and only for a subset of the noise and graph configurations. The additional comparison with sortnregress, that would highlight the triviality of certain settings, is also reported separately from the main results and it again covers only a subset of the noise and graph configurations.\n2. After discussion, the authors added the literature on order-based methods to their related work and reported a comparison with one of them (DAGuerreotype) in the appendix. However, this comparison is again carried out on trivial settings. Moreover the authors justified not thoroughly comparing with this line of works because they are supposedly not scalable. However one of the VP-DAG variants is one of the few methods with quadratic complexity (like the proposed method) so it is not clear how the authors can make this claim without empirical evidence.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "WhK69CAB0P", "reviewer": "Reviewer_t9PA", "summary": "The paper aims to improve the computational complexity and scalability of differentiable structure learning method. A smooth acyclicity formulation based on priority vector and tempered sigmoid are presented. With this formulation, the final approach solves an unconstrained optimization problem, requiring quadratic number of operations per optimization step. The paper presents empirical analysis to study the performance and execution time.", "questions": "For some questions, see weaknesses section.\n\nOther questions:\n- Similar to NOTEARS and DAGMA, one would expect an $L_1$ regularization $\\|W\\|_1=\\|H\\circ S(p)\\|_1$. Instead, Eq. (12) involves three different regularizations $\\|H\\|_1$, $\\|H\\|_2$, and $\\|p\\|_2$, two of which are $L_2$ regularizations.  This seems counter-intuitive and involves many hyperparameters. Can the authors explain why $\\|W\\|_1=\\|H\\circ S(p)\\|_1$ would not work?\n- In practice, since the temperature cannot be annealed to exactly 0, will the solution be exactly an acyclic graph before thresholding?\n- Does the characterization in Eq. (7) cover the entire space of weighted DAGs $W$? \n- It seems that the references provided for Section 2 \"NOBEARS\" and Section 1 \"does not sacrifice theoretical guarantees of asymptotic convergence to acyclic solution ...\" are not correct?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "-   The goal of improving the computational complexity and scalability of differentiable structure learning method is well motivated.\n-   The smooth acyclicity formulation and optimization scheme are sound.\n-   The execution time and scalability are improved.", "weaknesses": "The novelty of the paper may be limited because:\n- Although the paper claimed several times that it provides \"the first unconstrained optimization approach that learns DAG from observational data\", several existing works have done so, such as https://arxiv.org/abs/2006.10201 and DP-DAG (https://arxiv.org/abs/2203.08509). DP-DAG has been discussed in Section 2, but it may not be clear why the paper still claims being \"first\". Also, it would be great to compare these methods in the experiments.\n- A key motivation of the work is to avoid computing the acyclicity constraint with cubic cost, which has been achieved by DP-DAG.\n\nFurthermore, the overall idea, to my knowledge, share similarities NOCURL and is not entirely new, which has  been well acknowledged and discussed in Section 4.2 and Appendix D.2.\n\nThere is quite some performance drop for the proposed method in many cases (especially for large $d$), although the improvement in running time is notable:\n- It may not be clear how to interpret AUC as compared with NHD. For instance, for the results of 500 nodes in Appendix E.3, the NHD of COSMO and DAGMA are 3.57 and 1.73 respectively, indicating that COSMO performs much worse than DAGMA (probably due to lots of false discoveries). However, the AUC of COSMO is slightly better than DAGMA.\n- Therefore, I would suggest reporting NHD in the main paper that might be more commonly used in structure learning (e.g. see DAGMA and NOTEARS papers).\n- It seems that there is much more performance drop for nonlinear case. Is there a reason for not including NOTEARS-MLP for 40 and 100 nodes?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)