Differentiable Causal Discovery with Smooth Acyclic Orientations

Riccardo Massidda 1 Francesco Landolfi 1 Martina Cinquini 1 Davide Bacciu 1

Abstractstrained optimization scheme to learn acyclic causal graphsfrom data in quadratic time. We represent Directed Acyclic
Graphs (DAGs) as a differentiable function of a directed
graph and a priority vector. By applying a tempered sigmoid
to the pair-wise priority differences, we define a smooth relaxation of an acyclic orientation matrix where each node
has an outgoing arc to all nodes with higher priority. In particular, our smooth orientation matrix is equal to the discreteorientation whenever the sigmoid temperature tends to zero.
Thus, by annealing the temperature during training, the solution of the optimization problem is acyclic by construction.
Further, since our approach only requires a quadratic number of operations per optimization step, its constraint-freescheme is significantly faster than existing methods. Overall,
we refer to our approach as COSMO, for Causal Ordering
Discovery with SMooth Acyclic Orientations.

Most differentiable causal discovery approaches
constrain or regularize an optimization problem
using a continuous relaxation of the acyclicity
property. The cost of computing the relaxation is
cubic on the number of nodes and thus affects thescalability of such techniques. In this work, we introduce COSMO, the first quadratic and constraintfree continuous optimization scheme. COSMO
represents a directed acyclic graph as a priority
vector on the nodes and an adjacency matrix. We
prove that the priority vector represents a differentiable approximation of the acyclic orientation of
the graph, and we demonstrate the existence of an
upper bound on the orientation acyclicity. In addition to being asymptotically faster, our empiricalanalysis highlights how COSMO performs comparably to constrained methods for graph discovery.

2. Background and Related Works
Graphical causal models represents causal relations between
d variables as a DAG. In the linear case, a Structural Causal
Model (SCM) consists of a weighted acyclic adjacency matrix W ∈ Rd×d such that, for any random vector x ∈ Rd ,it holds x = W⊤x + z, given an unobserved noise vector z.
NOTEARS (Zheng et al., 2018) formalizes causal discovery
as the constrained optimization problem

1. Introduction
Graphical approaches, such as Structural Causal Models(SCMs), emerged as the dominant framework to represent
causal information about the world (Pearl, 2009). A fundamental problem in this context concerns the discovery of
causal relations between a set of variables, i.e., the problemof identifying which arcs exist between nodes associated to
the variables of interest (Spirtes et al., 2000). Continuous
causal discovery techniques approach the problem by optimizing an acyclic causal graph (Vowels et al., 2022). In
this context, a well-established methodology (Zheng et al.,2018) defines a smooth relaxation of the acyclicity property
and frames causal discovery as a constrained optimization
problem. Despite their widespread adoption, exact acyclicity constraints impose a cubic number of operations in thenumber of nodes that severely affects scalability.

min L(X, W⊤X) + λ∥W∥1(1)s.t. tr(eW◦W ) − d = 0(2)

W∈Rd×dwhere L is the Mean Squared Loss (MSE). In particular, theconstraint equals zero if and only if the weighted adjacency
matrix W is acyclic. The authors propose to solve theproblem using the Augmented Lagrangian method (Nocedal& Wright, 1999), which in turn requires to solve multiple
unconstrained problems and to compute the constraint value
at each optimization step. Because of the matrix exponential,
computing the constraint requires O(d3 ) operations.

Motivated by such limitation, we propose the first uncon-

Subsequent works mostly tried to tackle the computational
cost of NOTEARS by either: (i) replacing the matrix exponential with theoretically cubic but practically faster operations,(ii) approximating the constraint,or (iii) defining an
unconstrained problem that starts from the solution of ashorter constrained problem.We report the relevant related
works along with their computational complexity in Table 1.

1

Department of Computer Science, University of
Pisa, Italy. Correspondence to: Riccardo Massidda <riccardo.massidda@phd.unipi.it>.
Published at the Differentiable Almost Everything Workshop of the
40 th International Conference on Machine Learning, Honolulu,
Hawaii, USA. July 2023. Copyright 2023 by the author(s).

1Table 1. Summary comparison of our proposal with related works.
To express computational complexity, we define d as the numberof nodes and k as the maximum length of iterative approaches.
[†]: NOCURL requires a preliminary solution obtained from a cubicexpensive constrained optimization problem.
MethodComplexity

NOTEARS (Zheng et al., 2018)
DAGMA (Bello et al., 2022)
NOBEARS (Lee et al., 2019)
TMPI (Zhang et al., 2022)
NOCURL (Yu et al., 2021)
COSMO (Ours)

3

O(d )
O(d3 )
O(kd2 )
O(kd2 )O(d2 )†O(d2 )

X1

X2

X1

X3

t→∞

Exact
Exact
ApproximatedApproximatedPartial
None

X2

X3

t→0

Figure 1. (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 arcs respecting the priority ordering. (right) In thelimit of the temperature to zero, the smooth orientation matrix
contains non-zero entries if and only if the arc respects the order,i.e., it directs a node to another with sufficiently higher priority.

Notably, while recent literature disputes the relevance ofplainly minimizing MSE to uncover causal relations in realworld scenarios (Reisach et al., 2021; Kaiser & Sipos, 2022).
NOTEARS -like formulations still have a relevant role asa building block in more complex continuous discovery
approaches (Lachapelle et al., 2020; Brouillard et al., 2020).

While priority vectors enable the representation of strict
partial orders in a continuous space, the construction of
the orientation matrix still requires the non-differentiableevaluation of the inequality in Equation 3. To this end,
we approximate the comparison of the difference against
the threshold ε, using a tempered sigmoidal function. We
refer to such approximation of the orientation matrix as thesmooth orientation matrix.
Definition 3.2 (Smooth Orientation Matrix). Let p ∈ Rd bea priority vector, ε > 0 be a strictly positive threshold, andt > 0 be a strictly positive temperature. Then, the smoothorientation matrix of the strict partial order ≺p,ε is the realmatrix St,ε (p) ∈ Rd×d such that, for any u, v ∈ V , it holds

3. Learning Acyclic Orientations with COSMO
To continuously represent the space of d-dimensional DAGs,
we introduce a priority vector p ∈ Rd . on its vertices V .
Consequently, given a strictly positive threshold ε > 0, we
define the strict partial order ≺p,ε as(3)St,ε (p)uv = σt,ε (pv − pu ),for any u, v ∈ V . In other terms, a vertex u precedes
another vertex v if and only if the priority of v is sufficiently
larger than the priority of the vertex u. Notably, with a zerothreshold ε = 0, the relation would be symmetric and thus
not a strict order. On the other hand, whenever ε is strictly
positive, we can represent a subset of all strict partial orders
sufficient to express all possible DAGs.

(6)

where σt,ε is the ε-centered tempered sigmoid, defined asσt,ε (x) =1.
1 + e−(x−ε)/t(7)Intuitively, the threshold ε shifts the center of the sigmoid
and breaks the symmetry whenever two variables approximately have the same priority. The temperature parametert > 0 regulates instead the steepness of the sigmoid. Because of the asymmetry introduced by the threshold, in thelimit of the temperature to zero, the zero-entries of a smoothorientation matrix coincide with the zero-entries of the corresponding orientation matrix (Figure 1). Therefore, we
prove that any directed acyclic graph can be represented as
the element-wise product of a directed adjacency matrix anda smooth orientation. Further, any directed graph resulting
from this decomposition is acyclic.
Theorem 3.3. Let W ∈ Rd×d be a real matrix. Then, forany ε > 0, W is the weighted adjacency matrix of a DAG
if and only if it exists a priority vector p ∈ Rd and a realmatrix H ∈ Rd×d such that

Lemma 3.1. Let W ∈ Rd×d be a real matrix. Then, forany ε > 0, W is the weighted adjacency matrix of a DAG
if and only if it exists a priority vector p ∈ Rd and a realmatrix H ∈ Rd×d such that
W = H ◦ T≺p,ε ,X1

X3

Constraintu ≺p,ε v ⇐⇒ pv − pu ≥ ε

X2

(4)

where T≺p,ε ∈ {0, 1}d×d is a binary orientation matrix
such that(
1 if u ≺p,ε v
T≺p,ε [uv] =(5)0 otherwise,
for any u, v ∈ V .

W = H ◦ lim St,ε (p),Proof. We report the proof in Appendix A.1.

t→0

2

(8)

where St,ε (p) is the smooth orientation matrix of ≺p,ε .

temperature. Then, given the smooth orientation matrix
St,ε (p) ∈ Rd×d , it holds

Proof. We report the proof in Appendix A.2.

h(St,ε (p)) ≤ edα − 1,

where h(St,ε (p)) = tr(eSt,ε (p) ) − d is the NOTEARS
acyclicity constraint and α = σt,ε (0) = σ(−ε/t).

Given our definition of smooth acyclic orientation, we caneffectively parameterize the space of DAGs as a continuousfunction of a weighted adjacency matrix and a priority vector. Formally, we propose to decompose the weight matrix
as the product
W = H ◦ St,ε (p),Proof. We report the proof in Appendix B.
To contrast the discovery of spurious causal relations, weapply L1 regularization on the adjacency matrix H in orderto perform feature selection. Further, during the annealing procedure, even if a vertex u precedes v in the partial
order ≺p,ε , the weight of the opposite arc v → u in thesmooth orientation matrix will only be approximately zero.
Therefore, sufficiently large values of the weighted adjacency matrix H, might still lead to cyclic paths during training. To avoid this issue, we regularize the L2-norm of the
non-oriented adjacency matrix.

(9)

where H ∈ Rd×d and p ∈ Rd . In line with previous work,
we propose to address causal discovery as a score-based
method minimizing a loss function on the observations. By
jointly learning adjacencies and priorities, we avoid theuse of acyclicity constraints and reduce to a unique unconstrained problem. Therefore, the computational complexityof our solution reduces to the construction of the weighted
matrix W, which can be achieved in O(d2 ) time and space
per optimization step.

Further, we also observe that the partial derivative∂Wuv /∂pu tends to zero whenever the priorities difference |pv − pu | tends to infinity. Therefore, we regularize the L2-norm of the priority vector. For the same reason, we initialize each component from the normal distribution pu ∼ N (0, ε2 /2), so that each difference follows thenormal distribution pv − pu ∼ N (0, ε2 ).

Notably, the smooth orientation matrix St,ε (p) represents anacyclic orientation only in the limit of the temperature t → 0.
Nonetheless, the gradient loss vanishes whenever the temperature tends to zero. In fact, for a given loss function L,
the gradient of the priority vector p has form
∂L(X; W⊤ X) X ∂L(X; W⊤ X) ∂Wuv
·
=
∂pu
∂Wuv
∂pu
v∈V

∂L(X; W⊤ X) ∂Wvu
+
·
,
∂Wvu
∂pu

Overall, we formalize COSMO as the differentiable and unconstrained problem
(10)

min

H∈Rd×d ,p∈Rd

lim

∂Wuv
Huv
= lim −
σt,ε (pv − pu )
t→0
∂pu
t
· (1 − σt,ε (pv − pu ))L(X, (H ◦ St,ε (p))⊤ X)(13)+ λ1 ∥H∥1 + λ2 ∥H∥2 + λp ∥p∥2 ,and, for each partial derivative ∂Wuv /∂pu , it holds that

t→0

(12)

where λ1 , λ2 , λp are the regularization coefficients for theadjacencies and the priorities. As the regularization coefficients λ = {λ1 , λ2 , λp }, the initial temperature tstart , thetarget temperature tend , and the shift ε are hyperparametersof our proposal. Nonetheless, Theorem 3.4 can guide thechoice of the final temperature value and the shift to guarantee a maximum tolerance on acyclicity. We delineate a possible strategy to model nonlinear relations in Appendix C.3.

(11)

= 0,
and, similarly, ∂Wvu /∂pu tends to zero.
To handle this issue, we tackle the optimization problem byprogressively reducing the temperature during training. In
practice, we perform cosine annealing from an initial positive value tstart to a significantly lower target value tend ≈ 0.
Further, we prove the existence of an upper bound on theacyclicity of the orientation matrix that is a monotone increasing function of the temperature t. Therefore, annealing
the temperature during training decreases the upper boundon the smooth orientation acyclicity.

4. Experiments
We discuss the experimental comparison of COSMO against
related state-of-the-art approaches that optimize a DAG
from data. Namely, we confront with the graph discovery
performance and execution time of NOTEARS (Zheng et al.,2020), NOCURL (Yu et al., 2021), and DAGMA (Bello et al.,2022). NOCURL consists of two phases: firstly it solves
a constrained and cubic-expensive optimization problem,
then it improves the solution with a further unconstrained
problem. We also compare with NOCURL - U, i.e., an enterely

Theorem 3.4. Let p ∈ Rd be a priority vector, ε > 0 be astrictly positive threshold, and t > 0 be a strictly positive3

Table 2. Experimental results on linear Erdős–Rényi simulated DAGs with different noise terms and sizes. For each algorithm, we reportmean and standard deviation over five independent runs. We highlight in bold the best result and in italic bold the second best result.
GaussdAlgorithm
COSMO
DAGMA

30

NOCURLNOCURL - U
NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS

500

COSMO
DAGMA
NOCURL - U

ExpGumbel

AUC

Time (s)AUCTime (s)AUCTime (s)0.984 ± 0.02
0.985 ± 0.01
0.967 ± 0.01
0.694 ± 0.06
0.973 ± 0.02

88 ± 2781 ± 192822 ± 15226 ± 55193 ± 170

0.989 ± 0.01
0.986 ± 0.02
0.956 ± 0.02
0.694 ± 0.05
0.966 ± 0.03

89 ± 3744 ± 75826 ± 24212 ± 55579 ± 284

0.914 ± 0.10
0.973 ± 0.02
0.915 ± 0.04
0.678 ± 0.05
0.981 ± 0.01

87 ± 2787 ± 86826 ± 17212 ± 55229 ± 338

0.961 ± 0.03
0.982 ± 0.01
0.962 ± 0.01
0.682 ± 0.05
0.963 ± 0.01

99 ± 2660 ± 1411664 ± 14
267 ± 1011000 ± 339

0.985 ± 0.01
0.986 ± 0.01
0.950 ± 0.02
0.693 ± 0.05
0.972 ± 0.01

99 ± 2733 ± 109
1655 ± 28
242 ± 410880 ± 366

0.973 ± 0.01
0.986 ± 0.01
0.962 ± 0.01
0.663 ± 0.04
0.969 ± 0.00

98 ± 1858 ± 101
1675 ± 34
247 ± 911889 ± 343

0.933 ± 0.01
0.980 ± 0.00
0.683 ± 0.05

436 ± 812485 ± 365
1546 ± 304

0.986 ± 0.00
0.984 ± 0.01
0.715 ± 0.03

390 ± 1022575 ± 469
1488 ± 249

0.982 ± 0.01
0.980 ± 0.00
0.728 ± 0.05

410 ± 1062853 ± 218
1342 ± 209

unconstrained variant of the latter that skips the preliminary
constrained phase.

Algorithm
COSMO
DAGMA
NOCURL-U
NOTEARS102

For each method, we perform causal discovery by minimizing the Mean Squared Error (MSE) of a model on observational data using the Adam optimizer (Kingma & Ba, 2015).
Then, we measure the Area under the ROC Curve (AUC) between the absolute value of the weighted adjacency matrix
and the ground truth graph (Heinze-Deml et al., 2018).

Time (s)101

100

By looking at the AUC of the identified causal graphs, weobserve that COSMO consistently achieves results comparable with constrained-optimization solutions such as DAGMA
or NOTEARS (Table 2). Furthermore, COSMO performs better than NOCURL on most datasets. As pointed out by its
authors, we also observe that the discovery performanceof NOCURL drops when optimizing the variable ordering
instead of inferring it from a preliminary solution. The
fact that COSMO outperforms NOCURL - U on all datasets
highlights the substantial role and effect of our smooth orientation formulation for learning the topological orderingof variables from data in an unconstrained way.

10−1
101

102

103

104

Nodes

Figure 2. Average duration of a training epoch for an increasing
number of nodes. We compute the optimization time on a random
ER-4 graph for five iterations over five independent repetitions.

vational data. Our novel definition of a smooth orientation
matrix ensures acyclicity of the solution without requiring
the evaluation of computationally expensive constraints. We
prove that annealing the temperature of our smooth acyclicorientation corresponds to decreasing an upper bound on
the widely adopted acyclicity relaxation from NOTEARS.
Our empirical analysis showed that COSMO performs comparably to constrained methods in significantly less time.
Furthermore, the analysis highlights the role of our parameterization, which does not incur the necessity of preliminary solutions and solves a unique unconstrained problem.
Overall, COSMO opens up more scalable and more causallygrounded continuous causal discovery strategies, without
sacrificing — as demonstrated in this work — the theoreticalguarantees on DAGs approximation capabilities.

Unsurprisingly, due to its quadratic computational complexity, COSMO is significantly faster than constrained methods on all datasets, especially for increasing graph sizes.
Other than requiring multiple optimization problems, allcompeting methods incur in an higher computational cost
per step (Figure 2). Therefore, already for graphs with
500 nodes, only COSMO, DAGMA, and NOCURL - U return asolution before hitting our wall time limit.

5. Conclusion
We introduced COSMO, an unconstrained and continuous
approach for learning an acyclic causal graph from obser4

References

Spirtes, P., Glymour, C., and Scheines, R. Causation, Prediction, and Search, Second Edition. Adaptive computationand machine learning. MIT Press, 2000.

Bello, K., Aragam, B., and Ravikumar, P. DAGMA:
Learning DAGs via M-matrices and a Log-Determinant
Acyclicity Characterization. In Advances in Neural Information Processing Systems, 2022.

Vowels, M. J., Camgoz, N. C., and Bowden, R. D’ya likedags? a survey on structure learning and causal discovery.
ACM Computing Surveys, 55(4):1–36, 2022.

Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien,S., and Drouin, A. Differentiable causal discovery
Yu, Y., Gao, T., Yin, N., and Ji, Q. Dags with no curl:from interventional data. In Larochelle, H., Ranzato,
An efficient dag structure learning approach. In InterM., Hadsell, R., Balcan, M., and Lin, H. (eds.),
national Conference on Machine Learning, pp. 12156–
Advances in Neural Information Processing Systems,
12166. PMLR, 2021.
volume 33, pp. 21865–21877. Curran Associates, Inc.,
2020. URL https://proceedings.neurips. Zhang, Z., Ng, I., Gong, D., Liu, Y., Abbasnejad, E., Gong,
M., Zhang, K., and Shi, J. Q. Truncated matrix power itercc/paper_files/paper/2020/file/ation for differentiable dag learning. Advances in Neural
f8b7aa3a0d349d9562b424160ad18612-Paper.
Information Processing Systems, 35:18390–18402, 2022.
pdf.

Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. P.
Dags with no tears: Continuous optimization for structurelearning. Advances in neural information processing
systems, 31, 2018.

Heinze-Deml, C., Maathuis, M. H., and Meinshausen, N.
Causal structure learning. Annual Review of Statistics
and Its Application, 5:371–391, 2018.
Kaiser, M. and Sipos, M. Unsuitability of notears for causal
graph discovery when dealing with dimensional quantities. Neural Processing Letters, 54(3):1587–1595, 2022.

Zheng, X., Dan, C., Aragam, B., Ravikumar, P., and Xing,
E. Learning sparse nonparametric DAGs. In Proceedingsof the Twenty Third International Conference on Artificial Intelligence and Statistics, pp. 3414–3425. PMLR,
2020. URL https://proceedings.mlr.press/
v108/zheng20a.html. ISSN: 2640-3498.

Kingma, D. P. and Ba, J. 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. URL
http://arxiv.org/abs/1412.6980.
Lachapelle, S., Brouillard, P., Deleu, T., and Lacoste-Julien,
S. 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.net/
forum?id=rklbKA4YDS.
Lee, H.-C., Danieletto, M., Miotto, R., Cherng, S. T., andDudley, J. T. Scaling structural learning with no-bears toinfer causal transcriptome networks. In PACIFIC SYMPOSIUM ON BIOCOMPUTING 2020, pp. 391–402. World
Scientific, 2019.
Lopez, R., Hütter, J.-C., Pritchard, J., and Regev, A. Largescale differentiable causal discovery of factor graphs. Advances in Neural Information Processing Systems, 35:
19290–19303, 2022.
Nocedal, J. and Wright, S. J. Numerical optimization.
Springer, 1999.
Pearl, J. Causality. Cambridge university press, 2009.
Reisach, A., Seiler, C., and Weichwald, S. Beware of thesimulated dag! causal discovery benchmarks may be easy
to game. Advances in Neural Information Processing
Systems, 34:27772–27784, 2021.
5

A. Deferred Proofs
A.1. Proof of Lemma 3.1
Lemma 3.1 Let W ∈ Rd×d be a real matrix. Then, for any ε > 0, W is the weighted adjacency matrix of a DAG if andonly if it exists a priority vector p ∈ Rd and a real matrix H ∈ Rd×d such that
W = H ◦ T≺p,ε ,(14)

where T≺p,ε ∈ {0, 1}d×d is a binary orientation matrix such thatT≺p,ε [uv] =(1
0

if u ≺p,ε v
otherwise,(15)

for any u, v ∈ V .
Proof. Firstly, we prove the existence of a priority vector p and an adjacency matrix H for each weighted acyclic matrix W
of a directed acyclic graph D = (V, A). Being a DAG, the arcs follow a strict partial order ≺ on the vertices V = {1, . . . , d}.
Therefore, it holds thatA ⊆ {(u, v) | u ≺ v}.

(16)

Consequently, for an arbitrary topological ordering of the variables π : V → {1, . . . , d}, which always exists on DAGs, wedefine the vector p ∈ Rd such thatpu = επ(u).

(17)

Given the following implicationsu ≺ v =⇒ π(v) > π(u)(18)=⇒ pv − pu = ε(π(v) − π(u)) ≥ ε⇐⇒ u ≺p,ε v,

(19)
(20)

it holds that the order ≺p,ε contains the order ≺. Finally, we can define the adjacency matrix as H = W, where
W = H ◦ T≺p,ε holds since T≺p,ε [u, v] = 0 only if (u, v) ̸∈ A.
To prove that any priority vector p and adjacency matrix H represent a DAG, we first notice that, since the arcs follow astrict partial order, the orientation T≺p,ε is acyclic. Then, by element-wise multiplying any matrix H we obtain a sub-graphof a DAG, which is acyclic by definition.
A.2. Proof of Theorem 3.3
Theorem 3.3 Let W ∈ Rd×d be a real matrix. Then, for any ε > 0, W is the weighted adjacency matrix of a DAG if andonly if it exists a priority vector p ∈ Rd and a real matrix H ∈ Rd×d such that
W = H ◦ lim St,ε (p),
t→0

(21)

where St,ε (p) is the smooth orientation matrix of ≺p,ε .
Proof. By Lemma 3.1, we know that for any acyclic weighted adjacency matrix W there exist a priority vector p and a realmatrix H such that W = H ◦ T≺p,ε . Further, by Definition 3.2, the inner limit of Equation 21 solves to
1
lim St,ε (p)uv = 1/2
t→0


0
6

pv − pu > ε
pv − pu = ε
pv − pu < ε.

(22)

Therefore, we can define H′ ∈ Rd×d such that
Huv =

(

2H′uv
H′uv

pv − pu = ε,,otherwise.

(23)

from which
W = H ◦ T≺p,ε = H′ ◦ lim St,ε (p).
t→0

(24)

Then, to prove the counter-implication of Theorem 3.3, we notice thatlim St,ε (p)uv = 0 ⇐⇒ pv − pu < ε ⇐⇒ u ̸≺p,ε v.

t→0

(25)

Therefore, since the smooth orientation contains an arc if and only if the vertex respect the strict partial order ≺p,ε , it isacyclic. Consequently, as in Lemma 3.1, the element-wise product with an acyclic matrix results in a sub-graph of a DAG,
which is also acyclic by definition.
A.3. Priority Vector Initialization
By independently sampling each priority component from a Normal distribution N (µ, s2 /2), each difference is consequently
sampled from the distribution N (0, s2 ). Therefore, we seek a value for the standard deviation s that maximizes the partial
derivative
Huv
∂Wuv
=
σt,ε (pv − pu )(1 − σt,ε (pv − pu )).
∂pu
t

(26)

for arbitrary vertices u, v. Given the definition of the tempered-shifted sigmoid function, this object has maximum
in pv − pu = ε. Therefore, by setting the variance as s2 = ε2 , we maximize the density function of the point pv − pu = ε
in the distribution N (0, s2 ).

B. Smooth Acyclic Orientations and the Acyclicity Constraint
In this section, we present the proof for the upper bound on the acyclicity of a smooth acyclic orientation matrix. To thisend, we introduce two auxiliary and novel lemmas. Firstly, we introduce a lemma which binds the product of a sigmoid on asequence of values with zero sum (Lemma B.1). Then, we introduce another lemma on the sum of the priority differences ina cyclic path (Lemma B.2). Finally, we are able to prove the acyclicity upper bound from Theorem 3.4.
Lemma B.1. (Sigmoid Product Upper Bound) Let {xi } be a sequence of n real numbers such thatnXxi = 0.

i=1Then, for any temperature t > 0 and shift ε ≥ 0, it holds thatnYσt,ε (xi ) ≤ αn ,i=1where α = σ(−ε/t) is the value of the tempered and shifted sigmoid in zero.
Proof. Before starting, we invite the reader to notice that, for any temperature t > 0, if the sum of the sequence {xi } iszero, then also the sequence {xi /t} sums to zero. Therefore, we omit the temperature in the following proof, and assume todivide beforehand all elements of the sequence by the temperature t. Further, we explicitly denote the shifted sigmoid byusing the notation σ(xi − ε).
Firstly, we formulate the left-side of the inequality as
Pn
Qn xi −εnnYYee i=1 xi −εe−nε
exi −ε
i=1
Q
Q
Q
=
=
=
.
σ(xi − ε) =
n
n
n
xi −ε
xi −ε
xi −ε
1 + exi −εi=1 1 + e
i=1 1 + e
i=1 1 + e
i=1
i=1
7

Similarly, we rewrite the right side as
αn = σ(−ε)n =



e−ε
−ε
i=1 1 + e

Qn

n

=

e−nε
.
(1 + e−ε )nTherefore, proving the left-side smaller or equal than the right-side, reduces to proving the left-denominator is larger than
the right-denominator. Formally,nY

1 + exi −ε ≥ 1 + e−ε

i=1

n

,

or equivalently, by applying the logarithmic function,nXlog(1 + exi −ε ) ≥ n log(1 + e−ε ).

(27)

i=1

To further ease the notation, we refer to the left side of inequality 27, as the target function
T (x) =nX

log(1 + exi −ε ).

i=1

In particular, to prove 27, we show thatmin T (x) = n log(1 + e−ε ),x(28)

for x = ⃗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 x1 on the target function TP(x). To constraint thencomponents sum to zero, we consider the components {xi } for i > 2 as free, and then x2 = −x1 − i=3 xi as a functionof the remaining. The choice of x1 , x2 is independent from the components ordering, and thus applies to any possible pair.
Consequently,PnPn∂(log(1 + ex1 −ε ) + log(1 + e−x1 − i=3 xi −ε ) + i=3 log(1 + exi −ε ))∂T (x)
=
(29)
∂x1
∂x1
Pn
∂(log(1 + ex1 −ε ) + log(1 + e−x1 − i=3 xi −ε )
(30)
=
∂x1
n
X
xi − ε).
(31)= σ(x1 − ε) − σ(−x1 −i=3

Since σ(−ε) = σ(−ε), the equation is satisfied, for any component xi , by x = ⃗0,
We finally prove Inequality 27, by showing that the value of the target function T (x), in its only stationary point x = ⃗0,equals the bound. Formally,T (⃗0) =nXlog(1 + e−ε )(32)= n log(1 + e−ε ).

(33)

i=1

Lemma B.2. (Sum of Differences in Cycle) Let {pi } be a sequence of n + 1 real numbers such that p1 = pn+1 . Then, let{δi } be a sequence of n numbers such that δi = pi+1 − pi . Then,nX

δi = 0.

i=1

8

(34)

Proof. The proof is immediate from the following sequence of equations:
n
X

δi =

i=1

n
X

pi+1 − pi = −p1 +

i=1

n
X

pi − pi + pn+1 = 0.

i=2Theorem B.3. (Orientation Acyclicity Upper Bound) Let W ∈ Rd×d be a real matrix. Then, for any ε > 0, W is theweighted adjacency matrix of a DAG if and only if it exists a priority vector p ∈ Rd and a real matrix H ∈ Rd×d such that
W = H ◦ lim St,ε (p),

(35)

t→0

where St,ε (p) is the smooth orientation matrix of ≺p,ε .

Proof. The left side of Inequality 35 corresponds to the following infinite seriestr(eP ) − d =
=

∞
X
1

k=0
∞
X

k=1

k!

tr(P(k) ) − d1tr(P(k) )k!

where P(k) is the matrix power defined as P(k) = P(k−1) P and P0 = I.
By definition of matrix power, the u-th element on the diagonal of P(k) equals to
(k)
Puu
=

X

Pv(k−1)
Pu,v1
1 ,u

X

X

v1 ∈V

=

···

Pu,v1

Pvi ,vi+1

i=1

vk−1 ∈V

v1 ∈V

k−2
Y

!

Pvk−1 ,u .

Intuitively, the u-th element on the diagonal of P(k) amounts to the sum of all possible paths starting and ending in thevariable Xu . Therefore, being the same node, the priority of the first and the last node in the path are equal by definition.
Consequently, by Lemma B.2, the difference between the priorities sums to zero. For this reason, given Lemma B.1, it holds
that the product of the corresponding sigmoids is smaller or equal than αk . Therefore,
P(k)
uu =

X

···

X

···

X

k−2
Y

Pvi ,vi+1

i=1αk

vk−1 ∈V

v1 ∈V

=d

Pu,v1

vk−1 ∈V

v1 ∈V

≤

X

k−1 k

α .

Consequently, we upper bound the trace of the orientation matrix power astr(P(k) ) =

d
X

u=1

9

k k
P(k)
uu ≤ d α .

!

Pvk−1 ,uFinally, we are able to prove the Theorem astr(eP ) − d =
=
≤

∞
X
1

k=0
∞
X

k=1
∞
X

k=1

k!

tr(P(k) ) − d1tr(P(k) )
k!
1 k k
d α
k!

= −1 + edα ,where the last passage is due to the Taylor series of the exponential function.

C. Implementation Details
We run all the experiments on our internal cluster of Intel(R) Xeon(R) Gold 5120 processors, totaling 56 CPUs per machine.
We report details on the evaluation (C.1), the data generation procedure (C.2), and the models (C.3).
C.1. 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 in the following subsections. Firstly, we sample fifty random
configurations from the hyperparameter space. Since the hyperparameter space of COSMO also includes temperature andshift values, we extract more hyperparameters (200 – 800). Due to the significant speedup of COSMO, hyperparametersearches take a comparable amount of time, with NOTEARS being significantly longer on small graphs as well. Then, we test
each configuration on five randomly sampled causal DAGs. We select the best hyperparameters according to the average
AUC value. Finally, we perform a validation step by running the best configuration on five new random graphs.
Following previous work, we recover the binary adjacency matrix A of the retrieved causal graph by thresholding thelearned weights W with a small constant ω = 0.3. Formally, A = |W| > ω.
C.2. Synthetic Data
As we remarked in the main body, continuous approaches are particularly susceptible to data normalization and might
exploit variance ordering between variables (Reisach et al., 2021). Therefore, empirical results on simulated datasets that do
not explicitly control this condition might not generalize to real-world scenarios. Nonetheless, our proposal aims at defining
a faster parameterization that could replace existing continuous approaches as a building block in more complex discovery
solutions. Therefore, as initially done by Zheng et al. (2018) and subsequent work, we tested COSMO and the remaining
baselines in the usual synthetic testbed without normalizing the variance.
We include in our code the exact data generation process from the original implementation of NOTEARS.1 Therefore, thedataset generation procedure firstly produces a DAG with either the Erdős–Rényi (ER) or the scale-free Barabási-Albert (SF)models. Then, it samples 1000 independent observations. In the linear case, the generator simulates equations of the formfi (x) = W⊤i x + zi ,(36)

where we sample each weight Wij from the uniform distribution U (−2, −0.5) ∪ (0.5, 2) and each noise term zi from either
the Normal, Exponential (λ = 1), or Gumbel (µ = 0, β = 1) distributions. In the non-linear case, we simulate an additivenoise model with formfi (x) = gi (x) + zi ,(37)

where gi is a randomly initialized Multilayer Perceptron (MLP) with 100 hidden units and the noise term zi is sampled fromthe Normal Distribution N (0, 1).
1

NOTEARS implementation is published with Apache license at https://github.com/xunzheng/notears.

10Table 3. Hyperparameter ranges and values for COSMO.

Hyperparameter

Range/Value

Learning Rate
λ1
λ2
λp
tstarttend
ε(1e-3, 1e-2)(1e-4, 1e-3)(1e-3, 5e-3)(1e-3, 3e-3)
0.45(5e-4, 1e-3)(5e-3, 2e-2)

C.3. Models
Since we focus on the role of acyclic learners as a building block within more comprehensive discovery solutions, we 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 PyTorch to exploit automatic differentiation. By
avoiding differentiation and other overhead sources, the time expenses results are not directly 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;
Lopez et al., 2022).
For the non-linear setting, similarly to NOTEARS (Zheng et al., 2020), we model the outcome of each variable Xu with aneural network fu : Rd → R, where we represent first-layer weights as a tensor H ∈ Rd×d×h . Intuitively, each entry Huvi
stands for the weight from the variable Xu to the i-th first-layer neuron in the MLP fv . Then, the weighted adjacency matrix
Ph
results from the summation on the hidden dimension. Formally, Wuv = i=1 Huvi .

By checking the convergence of the model, both NOTEARS, DAGMA, and NOCURL can dynamically stop the optimization
procedure. On the other hand, COSMO requires a fixed number of epochs in which to anneal the temperature value. Fora fair comparison, while we stop optimization problems after a maximum of 5000 training iterations, we do not disable
early-stopping conditions on the baselines. Therefore, when sufficiently large, the maximum number of epochs should
not affect the overall execution time of the methods. For COSMO, we 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 analgorithm whenever it hits a wall time limit of 20000 seconds.

As previously discussed in Subsection C.1, we perform a hyperparameter search on each model for each dataset. In particular,
we sampled the learning rate from the range (1e-4, 1e-2) and the regularization coefficients from the interval (1e-4, 1e-1).
For the specific constrained optimization parameters, such as the number of problems or decay factors, we replicated thebaseline parameters, for which we point the reader to the original papers or our implementation. For COSMO, we sample
hyperparameters from the ranges in Table 3, given our theoretical findings on the relation between acyclicity and temperature(Theorem 3.4), we ensure sufficiently small acyclicity values. In the non-linear variant, we employ Multilayer Perceptrons
with h = 10 hidden units for each variable.

D. Additional Results
In this section, we report further results on simulated causal DAGs with different noise terms, graph types, and increasing
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 Positive Ratio (TPR), and the Structural Hamming Distancenormalized 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 a 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 our wall time limit of 20000 seconds.
11D.1. ER4 - Gaussian Noised30

Algorithm

NHD

TPR

AUC

Time (s)

COSMO
DAGMA
NOCURLNOCURL - U

0.867 ± 1.01
0.707 ± 0.57
1.653 ± 0.17
5.623 ± 0.92
0.913 ± 0.60

0.953 ± 0.04
0.940 ± 0.04
0.942 ± 0.02
0.492 ± 0.08
0.940 ± 0.05

0.984 ± 0.02
0.985 ± 0.01
0.967 ± 0.01
0.694 ± 0.06
0.973 ± 0.02

88 ± 3781 ± 193
822 ± 15227 ± 55193 ± 170

1.388 ± 0.69
1.026 ± 0.40
5.226 ± 1.34
10.108 ± 4.11
2.380 ± 2.10

0.917 ± 0.04
0.876 ± 0.02
0.921 ± 0.02
0.427 ± 0.05
0.898 ± 0.03

0.961 ± 0.03
0.982 ± 0.01
0.962 ± 0.01
0.682 ± 0.05
0.963 ± 0.01

99 ± 2661 ± 142
1664 ± 15267 ± 1011001 ± 340

4.149 ± 1.14
2.246 ± 0.40
27.675 ± 16.52

0.819 ± 0.02
0.882 ± 0.01
0.410 ± 0.04

0.933 ± 0.01
0.980 ± 0.00
0.683 ± 0.05

437 ± 812485 ± 366
1546 ± 304

NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

D.2. ER4 - Exponential Noised

Algorithm

NHD

TPR

AUC

Time (s)

30

COSMO
DAGMA
NOCURLNOCURL - U
NOTEARS

0.600 ± 0.54
0.613 ± 0.91
2.300 ± 0.97
5.313 ± 0.17
1.320 ± 0.67

0.970 ± 0.02
0.958 ± 0.05
0.918 ± 0.04
0.423 ± 0.05
0.880 ± 0.10

0.989 ± 0.01
0.986 ± 0.02
0.956 ± 0.02
0.694 ± 0.05
0.966 ± 0.03

89 ± 3744 ± 75826 ± 24212 ± 55579 ± 284

COSMO

1.642 ± 0.26
1.294 ± 0.52
5.652 ± 1.35
11.642 ± 4.34
1.156 ± 0.44

0.952 ± 0.02
0.944 ± 0.02
0.854 ± 0.03
0.478 ± 0.05
0.904 ± 0.03

0.985 ± 0.01
0.986 ± 0.01
0.950 ± 0.02
0.693 ± 0.05
0.972 ± 0.01

99 ± 2733 ± 109
1655 ± 28
242 ± 410880 ± 366

2.342 ± 0.86
2.147 ± 1.08
20.183 ± 7.43

0.944 ± 0.02
0.902 ± 0.04
0.437 ± 0.03

0.986 ± 0.00
0.984 ± 0.01
0.715 ± 0.03

390 ± 1022575 ± 469
1488 ± 249

NHD

TPR

AUC

Time (s)2.220 ± 1.65
1.680 ± 0.73
3.873 ± 1.26
5.260 ± 0.57
0.587 ± 0.38

0.862 ± 0.14
0.937 ± 0.03
0.853 ± 0.08
0.475 ± 0.08
0.962 ± 0.03

0.914 ± 0.10
0.973 ± 0.02
0.915 ± 0.04
0.678 ± 0.05
0.981 ± 0.01

87 ± 2787 ± 86826 ± 17212 ± 55229 ± 338

2.398 ± 0.70
1.132 ± 0.79
4.714 ± 1.77
6.914 ± 0.80
1.402 ± 0.40

0.936 ± 0.02
0.921 ± 0.04
0.905 ± 0.03
0.383 ± 0.04
0.869 ± 0.04

0.973 ± 0.01
0.986 ± 0.01
0.962 ± 0.01
0.663 ± 0.04
0.969 ± 0.00

98 ± 1858 ± 101
1675 ± 34
247 ± 911889 ± 343

3.574 ± 1.44
1.737 ± 0.64
18.182 ± 9.28

0.932 ± 0.02
0.871 ± 0.03
0.462 ± 0.06

0.982 ± 0.01
0.980 ± 0.00
0.728 ± 0.05

410 ± 1062853 ± 218
1342 ± 209

100

500

DAGMA
NOCURLNOCURL - U
NOTEARSCOSMODAGMANOCURL - U

D.3. ER4 - Gumbel Noised

Algorithm

30

DAGMA
NOCURLCOSMO

NOCURL - U
NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

12

D.4. SF4 - Gaussian Noised30

Algorithm

NHD

TPR

AUC

Time (s)

COSMO
DAGMA
NOCURLNOCURL - U

0.300 ± 0.09
0.360 ± 0.30
0.967 ± 0.43
4.410 ± 0.72
0.553 ± 0.54

0.973 ± 0.01
0.973 ± 0.02
0.893 ± 0.03
0.566 ± 0.11
0.944 ± 0.06

0.997 ± 0.00
0.996 ± 0.01
0.983 ± 0.01
0.741 ± 0.08
0.984 ± 0.02

89 ± 5653 ± 198
828 ± 23226 ± 75292 ± 261

0.482 ± 0.31
0.712 ± 0.33
2.030 ± 0.46
5.521 ± 0.61
0.280 ± 0.35

0.962 ± 0.02
0.951 ± 0.02
0.883 ± 0.03
0.596 ± 0.09
0.972 ± 0.04

0.991 ± 0.01
0.995 ± 0.00
0.982 ± 0.01
0.788 ± 0.06
0.993 ± 0.01

99 ± 3479 ± 75
1667 ± 25
269 ± 910112 ± 492

1.566 ± 0.68
1.343 ± 0.46
7.146 ± 3.19

0.953 ± 0.02
0.915 ± 0.04
0.504 ± 0.08

0.989 ± 0.01
0.992 ± 0.00
0.780 ± 0.07

541 ± 151345 ± 33
1394 ± 217

NOTEARSCOSMODAGMANOCURLNOCURL - U

100

NOTEARS
COSMO

500

DAGMA
NOCURL - U

D.5. SF4 - Exponential Noised

Algorithm

NHD

TPR

AUC

Time (s)

30

COSMO
DAGMA
NOCURLNOCURL - U
NOTEARS0.613 ± 0.39
0.127 ± 0.20
0.887 ± 0.21
4.067 ± 0.73
0.513 ± 0.30

0.965 ± 0.02
0.991 ± 0.01
0.845 ± 0.02
0.460 ± 0.15
0.962 ± 0.03

0.985 ± 0.02
0.999 ± 0.00
0.985 ± 0.01
0.685 ± 0.09
0.984 ± 0.01

87 ± 2592 ± 200824 ± 25
212 ± 75189 ± 271

COSMO

0.724 ± 0.71
0.586 ± 0.56
1.998 ± 0.40
5.912 ± 1.54
0.910 ± 0.43

0.963 ± 0.04
0.969 ± 0.03
0.907 ± 0.03
0.575 ± 0.06
0.962 ± 0.02

0.985 ± 0.02
0.995 ± 0.00
0.980 ± 0.00
0.783 ± 0.04
0.991 ± 0.01

100 ± 2395 ± 1081670 ± 28
245 ± 710243 ± 723

1.445 ± 0.58
1.653 ± 0.91
12.140 ± 7.84

0.950 ± 0.03
0.873 ± 0.08
0.482 ± 0.08

0.990 ± 0.01
0.988 ± 0.01
0.727 ± 0.06

517 ± 108
1466 ± 247
1205 ± 257

NHD

TPR

AUC

Time (s)0.467 ± 0.51
0.487 ± 0.20
0.747 ± 0.19
3.107 ± 0.64
0.860 ± 0.76

0.962 ± 0.05
0.956 ± 0.03
0.938 ± 0.02
0.460 ± 0.06
0.924 ± 0.06

0.990 ± 0.02
0.990 ± 0.01
0.989 ± 0.00
0.737 ± 0.04
0.975 ± 0.02

88 ± 2754 ± 179826 ± 32
213 ± 55199 ± 1300.864 ± 0.24
0.388 ± 0.30
1.806 ± 0.40
8.756 ± 2.65
1.134 ± 0.81

0.968 ± 0.01
0.975 ± 0.02
0.898 ± 0.03
0.550 ± 0.05
0.894 ± 0.08

0.992 ± 0.01
0.997 ± 0.00
0.982 ± 0.01
0.757 ± 0.03
0.989 ± 0.01

98 ± 2422 ± 103
1676 ± 31
245 ± 711618 ± 1309

1.426 ± 0.53
1.384 ± 0.38
8.931 ± 7.05

0.951 ± 0.03
0.849 ± 0.04
0.430 ± 0.10

0.994 ± 0.00
0.991 ± 0.00
0.741 ± 0.08

524 ± 221359 ± 34
1193 ± 229

100

500

DAGMA
NOCURLNOCURL - U
NOTEARSCOSMODAGMANOCURL - U

D.6. SF4 - Gumbel Noised

Algorithm
COSMO

30

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

13

D.7. ER6 - Gaussian Noised30

Algorithm

NHD

TPR

AUC

Time (s)

COSMO
DAGMA
NOCURLNOCURL - U

4.087 ± 1.12
2.367 ± 0.63
4.480 ± 0.92
7.490 ± 1.18
3.327 ± 1.65

0.838 ± 0.06
0.847 ± 0.03
0.869 ± 0.03
0.459 ± 0.08
0.840 ± 0.07

0.921 ± 0.04
0.958 ± 0.01
0.908 ± 0.03
0.672 ± 0.06
0.922 ± 0.04

89 ± 4665 ± 249909 ± 18226 ± 65239 ± 427

9.476 ± 3.01
10.740 ± 2.83
15.044 ± 1.60
30.719 ± 5.20
6.556 ± 3.10

0.771 ± 0.08
0.709 ± 0.13
0.785 ± 0.04
0.435 ± 0.03
0.842 ± 0.05

0.911 ± 0.05
0.902 ± 0.04
0.888 ± 0.02
0.580 ± 0.04
0.944 ± 0.02

98 ± 2761 ± 134
1687 ± 26
268 ± 912053 ± 940

25.443 ± 4.47
15.952 ± 1.67
165.465 ± 20.86

0.736 ± 0.01
0.553 ± 0.05
0.433 ± 0.02

0.937 ± 0.01
0.925 ± 0.01
0.558 ± 0.03

526 ± 100
3207 ± 271
1226 ± 293

NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

D.8. ER6 - Exponential Noised

Algorithm

NHD

TPR

AUC

Time (s)

30

COSMO
DAGMA
NOCURLNOCURL - U
NOTEARS

3.300 ± 0.95
3.480 ± 1.42
4.573 ± 0.78
8.700 ± 0.89
2.313 ± 1.55

0.897 ± 0.05
0.861 ± 0.06
0.846 ± 0.05
0.426 ± 0.07
0.881 ± 0.09

0.947 ± 0.03
0.945 ± 0.03
0.902 ± 0.03
0.615 ± 0.06
0.953 ± 0.04

89 ± 2672 ± 177897 ± 13226 ± 95516 ± 652

COSMO10.170 ± 2.74
8.118 ± 3.10
14.860 ± 4.67
30.600 ± 4.34
5.208 ± 2.54

0.768 ± 0.09
0.793 ± 0.11
0.685 ± 0.10
0.450 ± 0.04
0.796 ± 0.09

0.919 ± 0.04
0.934 ± 0.04
0.863 ± 0.06
0.591 ± 0.04
0.943 ± 0.03

99 ± 3681 ± 149
1735 ± 39
267 ± 812663 ± 1555

25.854 ± 4.28
16.417 ± 4.45
152.336 ± 31.97

0.741 ± 0.04
0.571 ± 0.11
0.425 ± 0.02

0.943 ± 0.01
0.925 ± 0.02
0.567 ± 0.03

460 ± 123
4069 ± 580
1363 ± 306

NHD

TPR

AUC

Time (s)2.840 ± 1.08
2.727 ± 0.83
5.003 ± 0.72
8.153 ± 0.96
2.740 ± 1.61

0.906 ± 0.04
0.906 ± 0.02
0.811 ± 0.04
0.422 ± 0.07
0.791 ± 0.10

0.954 ± 0.03
0.964 ± 0.02
0.891 ± 0.03
0.629 ± 0.04
0.938 ± 0.04

89 ± 3634 ± 194902 ± 9226 ± 65416 ± 446

10.048 ± 3.15
7.910 ± 3.05
11.932 ± 2.68
27.401 ± 4.42
4.884 ± 0.45

0.780 ± 0.07
0.805 ± 0.09
0.742 ± 0.04
0.431 ± 0.05
0.833 ± 0.05

0.899 ± 0.06
0.935 ± 0.04
0.894 ± 0.03
0.600 ± 0.04
0.951 ± 0.01

100 ± 3715 ± 203
1688 ± 34
266 ± 4
12634 ± 639

26.148 ± 4.86
16.358 ± 4.94
125.858 ± 36.61

0.740 ± 0.04
0.563 ± 0.07
0.367 ± 0.06

0.941 ± 0.02
0.921 ± 0.02
0.571 ± 0.02

418 ± 106
3527 ± 241
1612 ± 27

100

500

DAGMA
NOCURLNOCURL - U
NOTEARSCOSMODAGMANOCURL - U

D.9. ER6 - Gumbel Noised

Algorithm
COSMO

30

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

14

D.10. SF6 - Gaussian Noised30

Algorithm

NHD

TPR

AUC

Time (s)

COSMO
DAGMA
NOCURLNOCURL - U

1.273 ± 1.07
1.107 ± 0.37
1.573 ± 0.46
4.997 ± 0.98
0.933 ± 0.71

0.907 ± 0.10
0.930 ± 0.03
0.864 ± 0.04
0.506 ± 0.05
0.919 ± 0.05

0.963 ± 0.06
0.985 ± 0.01
0.973 ± 0.01
0.732 ± 0.05
0.984 ± 0.02

89 ± 2456 ± 39823 ± 14226 ± 85313 ± 184

4.478 ± 2.88
2.024 ± 0.71
2.824 ± 0.39
10.556 ± 6.00
1.412 ± 0.59

0.776 ± 0.15
0.914 ± 0.02
0.818 ± 0.02
0.542 ± 0.07
0.939 ± 0.03

0.874 ± 0.11
0.987 ± 0.00
0.980 ± 0.00
0.751 ± 0.08
0.990 ± 0.01

99 ± 2396 ± 53
1679 ± 27
266 ± 511156 ± 170

4.670 ± 1.99
3.825 ± 0.19
19.793 ± 11.03

0.912 ± 0.02
0.746 ± 0.03
0.368 ± 0.04

0.984 ± 0.00
0.982 ± 0.00
0.698 ± 0.04

460 ± 701418 ± 54
1137 ± 231

NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

D.11. SF6 - Exponential Noised

Algorithm

NHD

TPR

AUC

Time (s)

30

COSMO
DAGMA
NOCURLNOCURL - U
NOTEARS1.393 ± 1.24
1.147 ± 0.48
1.987 ± 0.54
4.787 ± 0.99
0.753 ± 0.49

0.926 ± 0.05
0.943 ± 0.03
0.757 ± 0.08
0.534 ± 0.07
0.943 ± 0.04

0.975 ± 0.03
0.982 ± 0.01
0.967 ± 0.01
0.761 ± 0.06
0.986 ± 0.01

88 ± 1578 ± 173820 ± 8227 ± 75312 ± 258

COSMO

3.836 ± 2.75
1.532 ± 0.61
2.890 ± 0.61
6.607 ± 1.05
1.784 ± 0.52

0.864 ± 0.09
0.887 ± 0.04
0.910 ± 0.02
0.474 ± 0.06
0.939 ± 0.02

0.944 ± 0.05
0.988 ± 0.00
0.977 ± 0.00
0.760 ± 0.06
0.988 ± 0.00

98 ± 2373 ± 881692 ± 28
266 ± 211369 ± 519

3.144 ± 0.47
3.854 ± 0.34
13.763 ± 8.79

0.919 ± 0.02
0.750 ± 0.01
0.389 ± 0.05

0.989 ± 0.00
0.977 ± 0.01
0.728 ± 0.06

457 ± 811384 ± 33
1436 ± 230

NHD

TPR

AUC

Time (s)1.047 ± 0.42
1.347 ± 0.63
1.787 ± 0.52
5.577 ± 0.43
1.053 ± 0.59

0.938 ± 0.03
0.933 ± 0.02
0.898 ± 0.02
0.549 ± 0.06
0.911 ± 0.04

0.984 ± 0.01
0.981 ± 0.01
0.969 ± 0.01
0.733 ± 0.05
0.978 ± 0.02

88 ± 1528 ± 67822 ± 29225 ± 45429 ± 251

3.486 ± 2.62
1.418 ± 0.34
3.074 ± 0.50
9.643 ± 4.59
1.586 ± 1.39

0.879 ± 0.10
0.910 ± 0.03
0.893 ± 0.02
0.464 ± 0.08
0.913 ± 0.06

0.947 ± 0.06
0.990 ± 0.00
0.976 ± 0.00
0.712 ± 0.10
0.987 ± 0.01

99 ± 2424 ± 901682 ± 22267 ± 911820 ± 985

3.288 ± 0.50
4.055 ± 0.88
56.103 ± 41.06

0.931 ± 0.01
0.802 ± 0.03
0.420 ± 0.06

0.992 ± 0.00
0.981 ± 0.00
0.648 ± 0.07429 ± 87
1465 ± 138
1201 ± 253

100

500

DAGMA
NOCURLNOCURL - U
NOTEARSCOSMODAGMANOCURL - U

D.12. SF6 - Gumbel Noised

Algorithm

30

DAGMA
NOCURLCOSMO

NOCURL - U
NOTEARS
COSMO

100

DAGMA
NOCURLNOCURL - U
NOTEARS
COSMO

500

DAGMA
NOCURL - U

15

D.13. ER4 - Non-linear MLPdAlgorithm

20
40
100

NHD

TPR

AUC

Time (s)

COSMO
DAGMA

2.110 ± 0.37
3.540 ± 0.45

0.726 ± 0.05
0.323 ± 0.06

0.939 ± 0.02
0.740 ± 0.05

125 ± 41834 ± 34

COSMO
DAGMA

2.420 ± 0.48
3.902 ± 0.07

0.668 ± 0.04
0.024 ± 0.02

0.930 ± 0.01
0.767 ± 0.03

139 ± 32054 ± 45

COSMO
DAGMA

4.864 ± 0.95
3.609 ± 0.10

0.627 ± 0.02
0.244 ± 0.02

0.918 ± 0.01
0.848 ± 0.02

228 ± 63537 ± 32

NHD

TPR

AUC

Time (s)DAGMA

1.725 ± 0.34
3.290 ± 0.15

0.701 ± 0.06
0.069 ± 0.04

0.949 ± 0.02
0.604 ± 0.11

125 ± 51838 ± 48

COSMO
DAGMA

2.065 ± 0.18
3.750 ± 0.00

0.573 ± 0.05
0.000 ± 0.00

0.960 ± 0.01
0.706 ± 0.09

139 ± 42069 ± 26

COSMO

3.019 ± 0.14
3.835 ± 0.04

0.374 ± 0.04
0.017 ± 0.01

0.945 ± 0.01
0.721 ± 0.05

229 ± 83551 ± 34

NHD

TPR

AUC

Time (s)D.14. SF4 - Non-linear MLP
d
20
40
100

Algorithm
COSMO

DAGMA

D.15. ER6 - Non-linear MLPdAlgorithm

20

COSMO
DAGMA

2.945 ± 0.23
5.430 ± 0.24

0.666 ± 0.07
0.132 ± 0.04

0.938 ± 0.01
0.607 ± 0.07

124 ± 51826 ± 45

40

COSMO
DAGMA

4.345 ± 0.64
6.000 ± 0.00

0.594 ± 0.06
0.000 ± 0.00

0.907 ± 0.01
0.648 ± 0.04

138 ± 52056 ± 40

COSMO4.390 ± 0.31
6.355 ± 0.22

0.421 ± 0.02
0.161 ± 0.02

0.900 ± 0.01
0.806 ± 0.01

228 ± 63534 ± 29

NHD

TPR

AUC

Time (s)2.110 ± 0.29
4.915 ± 0.05

0.673 ± 0.08
0.007 ± 0.01

0.965 ± 0.01
0.610 ± 0.13

124 ± 51846 ± 34

DAGMA3.255 ± 0.32
5.460 ± 0.02

0.541 ± 0.04
0.003 ± 0.00

0.950 ± 0.01
0.590 ± 0.08

138 ± 22056 ± 42

COSMO
DAGMA

4.625 ± 0.20
5.654 ± 0.05

0.305 ± 0.04
0.023 ± 0.01

0.936 ± 0.01
0.667 ± 0.07

229 ± 63480 ± 60

100

DAGMA

D.16. SF6 - Non-linear MLP
d
20
40
100

Algorithm
COSMO
DAGMA
COSMO

16