Analytic DAG Constraints for Differentiable DAG Learning
Authors: Zhen Zhang
Ignavier Ng
Dong Gong
Yuhang Liu
Mingming Gong
Biwei Huang
Kun Zhang
Anton van den Hengel
Javen Qinfeng Shi
Venue: ICLR 2025 Poster (2025)
Keywords: DAG
Causal Discovery
Structural Learning
Primary area: causal reasoning

Abstract: Recovering the underlying Directed Acyclic Graph (DAG)
  structures from observational data presents a formidable challenge, partly due
  to the combinatorial nature of the DAG-constrained optimization
  problem. Recently, researchers have identified gradient vanishing as
  one of the primary obstacles in differentiable DAG learning and have
  proposed several DAG constraints to mitigate this issue. By developing
  the necessary theory to establish a connection between analytic
  functions and DAG constraints, we demonstrate that analytic functions
  from the set $\\{f(x) = c_0 + \\sum_{i=1}^{\infty}c_ix^i | \\forall i > 0, c_i > 0; r = \\lim_{i\\rightarrow \\infty}c_{i}/c_{i+1} > 0\\}$ can be employed to
  formulate effective DAG constraints. Furthermore, we establish that
  this set of functions is closed under several functional operators,
  including differentiation, summation, and
  multiplication. Consequently, these operators can be leveraged to
  create novel DAG constraints based on existing ones. Using these
  properties, we design a series of DAG constraints and develop an
  efficient algorithm to evaluate them. Experiments
   in various settings demonstrate that our DAG constraints
  outperform previous state-of-the-art comparators. Our implementation is available at https://github.com/zzhang1987/AnalyticDAGLearning.

Published as a conference paper at ICLR 2025
ANALYTIC DAG CONSTRAINTS FOR DIFFERENTIABLE DAG LEARNING
Zhen Zhang'
Ignavier Ng²
Dong Gong?
Yuhang Liu'
Mingming Gong 4.5
Biwei Huang
Kun Zhang2.
Anton van den Hengel!
Javen Qinfeng Shi!
1 Australian Institute for Machine Learning, The University of Adelaide
2 Department of Philosophy, Carnegie Mellon University
3 School of Computer Science and Engineering, The University New South Wales
4 School of Mathematics and Statistics, The University of Melbourne
5 Department of Machine Learning, Mohamed bin Zayed University of Artificial Intelligence
6 Halicioglu Data Science Institute (HDSI), UC San Diego
ABSTRACT
Recovering the underlying Directed Acyclic Graph (DAG) structures from ob- servational data presents formidable challenge, partly due to the combinatorial nature of the DAG-constrained optimization problem. Recently, researchers have identified gradient vanishing as one of the primary obstacles in differentiable DAG learning and have proposed several DAG constraints to mitigate this issue. By developing the necessary theory to establish a connection between analytic fune- tions and DAG constraints, we demonstrate that analytic functions from the set (f(x) = co + Eie qa"llVi > 0, G > 0r lim;- Ga/ //441 > 0} can be employed to formulate effective DAG constraints. Furthermore, establish that this set of functions closed under several functional operators, including differen- tiation, summation, and multiplication. Consequently, these operators can be lever- aged to create novel DAG constraints based on existing ones. Using these properties, design series of DAG constraints and develop efficient algorithm evalu- ate them. Experiments in various settings demonstrate that our DAG constraints outperform previous state-of-the-art comparators. Our implementation available at htttts::// com/
INTRODUCTION
DAG learning aims to recover Directed Acyclic Graphs (DAGs) from observational data, which is core problem in many fields, including bioinformatics (Sachs et al. 2005 Zhang al. (013). machine learning (Koller and Friedman, 2009), and causal inference (Spirtes et al.) 2000). Under certain assumptions (Pearl (2000) Spirtes et al. 2000), the recovered DAGs could be interpreted causally (Koller and Friedman 2009).
There are two main categories of DAG learning approaches: constraint-based. and score-based methods. Most constraint based approaches, e.g., PC (Spirtes and Glymour 1991), FCI (Spirtes et al.
al.) 2012).. rely on conditional independence tests, which typically necessitate large sample size (Shah and Peters 2020 Vowels et al. (021)) The score-based approaches, including exact methods based on dynamic programming (Koivisto and Sood [2004|/iingh
(995)Colomm
[Singh and Moore 2005)
Silander and Myllymaki 2006), A* search (Yuan et al. 2011 Yuan and Malone| 2013). and integer programming (Cussens) 201 and greedy methods like GES (Chickering] 2002). model the validity of a graph according to some score function and are often formulated discrete optimization problems. A key challenge for acore-based methods is the super-exponential combinatorial search space of DAGs w.r.t number of nodes (Chickering) 1996; Chickering et al. 2004)..
Recently, Zheng et al. (2018) developed a continuous DAG learning approach using Langrange Multiplier methods and a differentiable DAG constraint based on the trace of the matrix exponential

---

Published as a conference paper at ICLR 2025
of the weighted adjacency matrix. The resulting method, named NOTEARS, demonstrated superior performance in estimating linear DAGs with equal noise variances. Very recently, [hhan et al.) (2022) and Bello al. (2022) suggest that main issue for NOTEARS and its derivatives, such Yu et al. (2019). is gradient vanishing for linear DAG models with equal variance. They have thus proposed new continuous DAG constraints by based geometric series of matrices well as log-determinant of matrices.
In fact, many of the proposed Directed Acyclic Graph (DAG) constraints can be unified, as demon- strated in Wei et al. 2020)) Wei et al. (2020) reveals that, for a d x d adjacency matrix, an order-d polynomial of matrices is necessary and sufficient to enforce the DAG property. However, from a computational standpoint, computing general matrix polynomials can be challenging. Considering the fact that infinite-order polynomials that converge, i.e., power series, can give rise to analytic functions that are often simpler to evaluate than general polynomials, it prompts the question of whether analytic functions could be utilized in constructing DAG constraints. Furthermore, raises the possibility of employing techniques commonly used for analyzing analytic functions in the investigation of continuous DAG constraints.
The answer is yes. We demonstrate that any analytic function within the class of functions denoted as F = (xx = co + Ei= c,a"; C > 0, Vi > 0; lim;-+so C/ c/c++1 > 0} can be utilized to formulate Directed Acyelic Graph (DAG) constraints. In fact, the DAG constraints introduced in Zheng et al. (218)) Zhang et al. (2022), and Bello et al. (2022) can all be interpreted as being based analytic functions from F.. Furthermore, closed under various function operators, including differentiation, function addition, and function multiplication. Leveraging this insight, we can construct novel DAG constraints based on pre- existing ones. Additionally, we can analyze the performance of these derived DAG constraints using techniques rooted in analytic functions.
establish that the function class F remains
PRELIMINARIES
DAG and Linear SEM Given directed acyclic graph (DAG) defined over random vector x 12, 222,... r.d
F.] T, the corresponding distribution P(x) is assumed to satisfy the Markov assumption (Spiress et al. 2000 Pearl 2000). We consider x to follow linear Structural Equation Model (SEM):
x B x
Here, B E ddd represents the weighted adjacency matrix that characterizes the DAG g, and [e,, E2, .......,G. ea] represents the exogenous noise vector, comprising d independent random variables. To simplify notation, we use (B) to denote the graph induced by the weighted adjacency matrix B, and we interchangeably use the terms 'random variables' and 'vertices' or 'nodes'.
We aim to estimate the DAG G from n i.i.d. observational examples of x, denoted by X E Raxd Generally, the DAG G can be identified only up to its Markov equivalence class under the faithfulness (Spites et al. 2000)) or the sparsest Markov representation assumption (Raskutti and Uhler 2018). It has been demonstrated that for linear SEMs with homoscedastic errors, where the noise terms are specified up to constant (Loh and Bühlmann 2013). and for linear non- Gaussian SEMs, where no more than one of the noise terms Gaussian (Shimizu et al., (0066), the true DAG can be fully identified. In our study, we specifically focus on linear SEMs with equal noise variances (Peters and Bühlmann 2013), where the scale of the data may be either known or unknown. When the scale is known, it is possible to fully recover the DAG, However, in the case of an unknown scale, the DAG may only be identified up to its Markov equivalence class.
Continuous DAG learning In recent years, series of continuous Directed Acyclic Graph (DAG) learning algorithms Bello et al. (202)) Ng et al. )2020): Zhang et al. 2002 Yu et al.) (2021 2019); Zheng et al.] (2018) have been introduced, demonstrating superior performance when applied to linear Structural Equation Models (SEMs) with equal noise variances and known data scales. These
methods can be expressed as follows:
arg min S(B, X), s.t. (BB) 0. B
Here, S is a scoring function, which can take the form of mean squared error (Zheng et al. 2018) or negative likeliiod (Ng et al. 2020). The function h is continuous and equals to 0 if and only if
2

---

Published as a conference paper at ICLR 2025
the weighted adjacency matrix B defines valid DAG. Previous approaches have employed various techniques, such as matrix exponential (Zheng et al. 2018), log-determinants (Bello et al. 2022). and polynomials (Zhang et al. 2022). to construct the function h. However, these methods are known to perform poorly, when applied to normalized data since they rely on scale information across variables for complete DAG recovery (Reisach al. 2021
ANALYTIC DAG CONSTRAINTS
In this section, demonstrate that the diverse set of continuous DAG constraints proposed in previous work can be unified through the use of analytic functions. We begin by providing a brief introduction to analytie functions and then illustrate how they be used to establish DAG constraints.
3.1 ANALYTIC FUNCTIONS AS DAG CONSTRAINTS
In mathematics, power series
f(x) Co
+
G
(3)
which converges for |a]| < r
lim /c//c+1/, defines an analytic function f on the open interval <-000
(-r,r), and r is known as the radius of convergence. When we replace r with a square matrix A, we obtain an analytic function of matrix follows:
S(A) coT A
where I is the identity matrix. Equation would converge the largest absolute value of the eigenvalues of A, known as the spectral radius and denoted by p(A), is smaller than r.
We are particularly interested in the following specific class of analytic functions
F - = co Σ ca';
Vi > 0, cs 0; Jim a///1 > 0},
(5)
5-000
1=1
as any analytic function belonging to FF be applied to construct continuous DAG constraint.
Proposition 1. Let B Bax" with p(B) < y be the weighted adjacency matrix of a directed graph G, and let f be an analytic function in the form of equation f where we further assume > 0 we have > 0, then G is acyclic and only if
function
[[B)] cod.
(6)
interesting property the DAG constraint equation Glis that its gradients can also be represented as transpose of an analytic function follows, which allows us to use analytic functions as the gradients of DAG constraints.
Proposition 2. There exists some real number where for all (B E Rdxd r), the derivative of [B)] w.r.t. B is
VBt [B]]] = [V.f(r)
notable that for d weighted adjacency matrix B, an order-d polynomial of B sufficient and necessary to enforce DAGness (Wei et al. 2020 Ng et al. 2022 Mcanwhile, evaluating matrix polynomials efficiently highly nontrivial (Higham 2008). For matrix analytic functions such as exponentials or logarithms, however, efficient algorithms exist (Higham 2008).
The connection between matrix analytic functions and real analytic functions means that various properties of the matrix function be obtained from simple real-valued function. To pursue DAG constraints with better computational efficiency, we seek an analytic function whose derivative can be represented by itself to reduce the computation of different analytic functions. If function has
3

---

Published as a conference paper at ICLR 2025
such property, various intermediate results be saved for future computation of gradients. The exponential function exp(x) with 3 exp(2)) (app)/rrrrrr = exp(œ), is a natural contender, and this leads to the well-known exponential-basee DAG constraints (Zheng et al. 2018)
00
Constraints: tr (B))
Σ B' B'/il= /il= d, Gradient: VB exp(B) = exp(B)",
(8)
x=
which will converge for any B.
Recently Bello et al. (2022) and Zhang et al. (2022) have suggested that exponential-buse DAG constraints suffer from gradient vanishing. One cause of gradient vanishing arises from the small coefficients high-order terms. The convergence radius for the exponential is oc, that is lim |c/c++1| lim (i 1)!/El| co, which suggests that, compared the lower order terms, i->00
the higher order terms contribute almost nothing in the DAG constraints, which indicates that it would not be efficient to prohibit possible long loops in candidate adjacency matrices.
Due to the fact that the adjacency matrix of DAG must form nilpotent matrix, whose spectural radius are acutally 0, naturally the spectral radius candidate adjacency matrices would be close to O. As result, do not need function with infinite convergence radius. Instead, we analytic function with finite convergence radius r = lim /c/c+t| < co. Thus by using a sequence
t-
with geometric progression = 11/8-11 or harmonic-geometro progression = 1/(it-11) can obtain two analytic functions,
fine (x) (s x.)
x² 1 fing (2) - log(s x.)
x
s log s.
(9)
is
i=0
i=1
Then by our Proposition and Proposition two dag constraints can be obtained as follows:
Constraints: tr/inu (B) = d, Gradient: V fatrfine (
[Sine (B)
(10a)
Constraints:tr B 0, Gradient: V tr fr B [f [B)]
(10b)
where truncated version of fin is applied in Zhang et al. (2022). and the figm based constraints are equivalent those in Bello et al. (2022). One key difference between Zhang al. (2022): ; Bello et al.)
(2022) and the exponential- radius, which requires an additional constraints (BB) < s. Meanwhile, the adjancency matrix of DAG must be nilpotent, and thus its spectral radius must 0. this case, such additional constraints would not affect the feasible set.
DAG constraints (Zheng et al., 2018) is their finite convergence
3.2 CONSTRUCTING DAG CONSTRAINTS BY FUNCTIONAL OPERATOR
One can easily observe a coincidence between flog and finv as follows,
a (z) o (2), [iog(z) x I" (in +C,
(1I)
Ox
which suggests that it may be possible to derive group of DAG constraints from an analytic function by applying integration or differentiation. This because derivatives of any order of an analytic function is also analytic. More formally, if function is analytic at some point Zo. then its n derivative for any integer n exists and is also analytic at To. Thus we can derive DAG constraints from any f E F as follows.
Proposition 3. Let f(x) cp -i qaxt E F be analytic on (-r,r), and let n be arbitary integer larger than 1, then B Badd with spectral radius p(B) < forms DAG and only if
0 f(x))) Fr
tr
n!c..-
(12)
The above proposition suggests that the differential operator be applied to analytic function to form a new DAG constraints. Besides the differential operator, the addition and multiplication of analytic functions also be applied to generate DAG constraints. That
4

---

Published as a conference paper at ICLR 2025
Proposition 4. Let f1(2) c/-
c C/ F, and f2(x) cJ+ G F. Then for ad-
jancency matrix B 1220d with spectral radius p(B) & min C // +0 of / [2+1 the following three statements are equivalent: 1) B forms a DAG; 2) tr (B) + (2 (c) + +g) 3) tr ( f2 (B) ctead.
Particularly for (1) and (x), due to the specific property of fi (x), we have
gu+1 Flo (2)
2th S'no Oxn+1
x (s x) (n-1) = (z) (
(13)
Or"
which implies that the n derivative of function 1/(s x) is propositional to the the order-(n + 1) power of 1/(s z). Using this property. the value of (I- - B/s) can be cached and then used to generate a series DAG constraints as well as their gradients. Similarly, the value of matrix
exponential exp(( (B)/s) can also be cached during the evaluation of DAG constraints to accelerate the computation. Furthermore, the gradients of the DAG constraints will also increase increases.
Proposition 5. Let be any positive integer, the adjancency matrix B {B E Rdxd ppBB) < s}
forms DAG if and only if
tr (II-B/s) d.
Furthermore, the gradients of the DAG constraints satisfies that vB {B E Rdxd pBB) < <s}
B/)r < B B/s)--~R111
where k is an arbitrary positive integer, and || || denote an arbitrary matrix norm induced by vector p-norm.
Gradient Vanishing and Numeric Stability For the series of DAG constraints constructed from Equation (13. as gradient vanishing is one of the main challenges for differentiable DAG learn- ing, according Proposition [ we may prefer larger n achieve better performance in practice. Furthermore, choosing a smaller may also help to amplify the gradient of DAG constraints. There- fore, |Blll al. (2022) applied an annealing strategy on s improve performance, while Zhang et al. (2022) used a fixed 8 = 1.0 in their implementation. However, in practice, especially when incorporating the DAG constraints with first-order optimizers, the spectral radius of the candidate B can often be larger than s. Bello al. 2022) applied a simple heuristics search for the proper while Zhilg et al. (2022) truncated the power series to avoid numerical issues in higher-order terms. However, in practice, we observed that [hhng et al. (2022) encountered some numerical issues for large graphs, and the simple heuristics used by Bello et al. (2022) may result in sacrifice in performance. Based our analysis, if B goes out, it can be verified by checking if the power series E (B/s) converges to (I - B B/s)= and can be chosen based on the spectral radius of B.
Efficiently Computation
The specific structure of the power series Ei= (B/s)* allows for fast
evaluation. Let
L, = Σ (B/s),
(14)
1=
then it is evident that
L L, + (B/s)LL Lt:
(15)
which indicates that the term L, can be obtained with O t) time complexity. Furthermore, using Equation 13). the gradient of tr(1 -B/s) can also be easily derived from Lo Along with the strategy for searching 8, we can use Algorithm [ to efficiently compute the DAG constraints.
3.3 OVERALL OPTIMIZATION FRAMEWORK
The DAG constraints above are applicable only to positive adjancency matrices, we use the Hadamard product to map a real adjancency matrix to a positive one. Thus Equation ( becomes:
arg min SBB X), s.t. tB B cod, p( @ B) <r,
(16)
B
5

---

Published as a conference paper at ICLR 2025
where the analytic function f(x) co Ei= x² E F,, and @ denotes the Hadamard product.
Algorithm 1 Efficient Evaluation of Gradients
A Algorithm 2 Path following algorithm
Input: B, s, d, 0,
Input: X RRxdd S: F; Pot (0, 1);
Output: Vatrsiog (B)or Vgtr [B)]
Touter; Tinner; 7 >
Output: Estimated B
I: + - B/s, W B/s
i: 0, 14 f- No. Bg
=1
2: for i = 0,1 < T Touter i + + do
3:
while
D(I B) I| > c and < 2d do
3:
B4 < B,
Optimize over
4: W - W x W
[|S ([ A1 B [1] 1/2/(B@)
5:
D (W +1)
for = 0: Tinerij + + do
=-2k
5:
B +- B,41 ④ BB4
7: end while
B;41 < Ba (B. +
8: if |ID - B) I| > then
A Aisign [B] VV @ B.
8 pB) + g. goto line
end for
10: else
1A /
11: For Flog return D /s
9: B B
12: For (B) return n|D /s]
10: end for
13:
end if
11: Return B
In our work, we choose to use the path-following approach with an & regularizer, as in Bello et al. (2022). This because the Lagrange approaches applied in Zhhnn al. (2022): Yu et al. (2021): Zheng ct al. (018);; Yu et al. (2019), the Lagrangian multiplier must be set very large value to enforce DAGness, which may result in numerical instability. In the path-following approach, instead of using large Lagrangian multipliers, small coefficients are added to the score function $
as follows
arg min piS(B, + AII B [1] tB o B), s.t.
<T,
(17)
where X is the user-specified weight for the & regularizer. For the additional constraints p O B) < with properly chosen initial value and step-length, it usually be satisfied. Also it is notable that | [[ < r is a sufficient condition for p B @ B) < and thus the sparsity constraints also encourage this condition to be satisfied. Based on Bello et al. (2022), we implemented a path-following shown in Algorithm
The optimization model equation 17 observed to perform well for linear Gaussian SEMs with equal variance as well as other equal variance SEMs. Meanwhile, for unequal variance, or normalized data from linear Gaussian SEMs with equal variance where the scale information are missing, MSE score function is not consistent and often provides misleading information about the underlying DAG. Additionally, as observed by Ng et al.) (2023), the initialization of adjacency matrices cases of unequal variance can significantly affect performance, suggesting that non-convexity may pose serious challenge in such scenarios.
NON-CONVEXITY ANALYSIS OF ANALYTIC DAG CONSTRAINTS
The non-convexity of function be analyzed through the analysis of its Hessian. Particularly for our analytic DAG constraints, its Hessian can be obtained using the following proposition and then the non-convexity can be analyzed by analysis the spectral radius of the Hessian.
Proposition 6. The Hessian of DAG constraints equation 6 can be obtained as follows:
B
Lr/
K
B
(18)
f=0
where & denotes the Kronecker product, and Kad E {0, 1f x is the commutation matrix satisfies that for any d x d matrix A
Ka,avec(A A) = veo( A
'hhe constant cod can he dropped because fr (B:3 3 B) is bounded below by cod, detailed derivation is provided in the supplementary file.

---

Published as a conference paper at ICLR 2025
Obviously, the Hessian Equation (18) is symmetric and not positive semi-definite. One widely used way to convexify Hessian is to find a positive scalar such that
A V tr f + nI,
(19)
becomes positive semi-definite. require 17 to be no less than the absolute value of the most negative eigenvalue of V2 B (B). Here the Hessian are symmetric matrix with all non-negative entries. For this kind of matrices the absolute value of the most negative eigenvalue of tr fr (B) is upper bounded by the spectral radius of the Hessian, and the bound is tight under certain conditions (Spielman 2012). Thus it would be nature to use the spectral radius of the Hessian to measure the level of non-convexity of the analytic DAG constraints.
The Hessian V2 trr((B) can be viewed as linear combinations of series of symmetrie matrices
dK Σ ==0
B'
B 2- with all non-negative entries. The commutation matrix Kdd Magnus
and Neudecker 1979) would not have any effects on the spectral radius as it is orthonormal. Thus
larger 6 would result the spectral radius of a single term ie; dd Σj= B
B
to
increase, and finally lead the spectral radius of the Hessian to increase the following proposition.
Proposition 7. For two analytic function fi()) = Cu,1 +0,1+EE Eizi 01,10* and fz(x) co,2 Ei= G,22', if > 1 we have Ct 2 C,22 > 0, then
p B tDf A (2 (B).
(20)
where pl-) denotes the spectral radius of matrix.
Proposition suggests that the spectral radius of the Hessian would increase the coefficients o; in the analytic function increase. This implies that DAG constraints with larger c; may gain benefits from gradient vanishing, but suffer from non-convexity. In fact, using Proposition Z it would be straightforward to get the following corollary, which provides the level of non-convexity comparison for several DAG constraints.
Corollary 8.
PVV tr (B)) < ( tr (B)) & Arrfiin (B)). B
(21)
The optimization problem equation 17 can be viewed convex objective plus one non-convex constraint, and the convex mean square error (MSE) loss may play different roles in different scenarios. For data with known scale, the MSE loss is consistent and thus it provides enough information to identify the underlying model and thus the non-convexity may not be a serious issue. This is because in the path-following optimization framework (provided in Algorithm 2). the heginning the optimization direction are dominated by the MSE loss so that it will push the candidates to point that is not far from global optimal. Thus DAG constraints with finite convergence radius is preferred to escape from gradient vanishing. Meanwhile, for DAG learning problem with unknown scale, the MSE loss may not be very informative to the underlying graph structure. In this case, the highly non-convex DAG constraints may lead to the optimizer to get trapped into local minimum easily, and thus we may need additional constraints to reduce the search space, which may possibly make the objective flatter. In our experiments, we find that by allowing only edges to exist between nodes with strong correlation can significantly improve the performance.
EXPERIMENTS
In the experiment, we compared the performance of different analytic DAG constraints in the same fath-following optimization framework. We implemented the path-following algorithm (provided in Algorithm 2) using cupy and numpy based on the path-following optimizer in Bello et al. (2022). For analytic DAG constraints with infinite convergence radius, consider the exponential based DAG constraints. For analytic DAG constraints with finite convergence radius, we consider the following different DAG constraints generated by the differentiation operator multiply operator:
Order-1: tr tr/ing (B B)
7

---

Published as a conference paper at ICLR 2025
Graphs
#Nodes
DAGMA
Order
Order 2
Order
Order
500
44,90
32,95
33,40
23.46
31.70 + 19.47
30.60 :: 19.07
29.80
20.97
ER2-Gausiiar
1000
94,80 1. 35,80
69.60 1: 27.64
55,60 1: 19.13
52.40 1: 19.86
57.30 +: 21.39
62.76
176.00
153.30
35.92
135.60
38,91
1: 28.65
500
125.30
:: 44.55
101.10
39.03
90.30 1: 39.56
93.90
4: 31.26
92.90
45.56
ER3 Gaussinn
1000
339.60
67.80
242,80
72.21
210.30
60.98
184,90
1: 47.44
165.80
1: 35,95
2000
669.50
140.61
610.70 = 136.84
555.30 + 106.01
479.50
= 88.72
424.90
+ 64.39
500
307.60
116.53
261.40 1: 102.81
263.00 1: 122.34
246.70 1: 110.28
223.80 1: 97.46
ER4-Gaussian
1000
878.50
174.96
689.20 165.62
695.50 + 134.41
619.80
150.43
626.30 157,59
2000
1922.30
187.69
1785,40 184.47
1779.30
211.23
1655,40 = 181.75
1574.10 = 152.23
500
58.20
31.58
40.50 26.93
28.90 + 16.60
31.00 +: 25.67
35.20 34.32
ER2-Exp
1000
93,90 + 33.96
68,70 # 23.20
54,00 + 16.26
50,90 + 17.29
57.90 + 24.93
500
142,70 = 50.13
106.00 + 39.15
95.10 + 32.68
99.60 + 37.94
100.30 = 47.52
ER3-Esp
1000
321.10 1: 83.82
242.80 1. 68.51
212.40
67.87
187.60 1: 61.87
173.10 1: 49.03
500
336.00 + 124.19
292.70 +: 123.41 718.20
294.90 + 130.66
254.40 + 133.05
214.70 : 84.29
ER4-Esp
879.40
162.98
1:
133.44
ER2-Gumbel
500
45.10
33,28
22,60
20.04
21.30 + 18:41
19.80 +: 16.39
16.20
11.91
1000
80.50
42.65
49,90 24.54
39.90 1: 14.04
36.80 1. 15.18
45.90
23.18
ER3-Gumbel
500
54.19
76.60 +: 60.77
60.60 +: 31.34
85.30
1000
297.90 1: 72.40
215.40 1. 52.35
185.00 71.98
173.90 1: 57.09
147.70 1: 40.51
BR4-Gumbel
500
338.80
127.56
273.70 + 131.13
257.50 + 11::111:11 111.06
232.40 + 121.98
234.70 + 149,59
1000
919.90
182.38
722.80 1: 177.86
734.80 1: 177.78
620.70 d: 187.56
564.10 1: 170.46
Table 1: DAG learning performance (measured in structural hamming distance, the lower the better, best results in bold) of different algorithms (2,3,,)4 graphs with different noise distributions. All our algorithms performs better than the previous state-of-the-ee DAGMA (Bello al. (022). and higher order DAG constraints suffers less gradient vanishing, it tends have better performance.
Graphs
#Nodes
DAGMA Bello
Orler
Order
Order
Order
500
31.40 + 43.51
24.30 + 43,90
32,40
49.38
34,20 + 45.56
41.50
48.45
SF2
1000
44.90 1 34.38
41.20 1: 36.02
22.50
13.21
29.20
1. 20.07
58.10
27.58
2000
189.80 +: 99.47
162.90 73.30
172. + 74.35
152.20
+: 90.29
172.60
124.12
500
58.10 1: 33,90
$1.10 1: 32.10
41.10
17.91
49,80 1: 24.58
71.00
23.46
169
158.10
46.70
: 55.25
162.50
57.54
065.40:2
75.60
2000
928.70 + 148.70
896.00 = 101.85
897.10
146.78
891.50 = 143.40
999.70
206.36
500
131.20 42.63
136.80 41,71
134,40 +: 39.34
128.90
36.68
151.60 +:: 37.05
SF4
1000
431,70 + 119.22
404.00 88.89
400.30 = 76.49
386.90
93.16
394.50 111.57
1525.10
299.02
297.88
1444.70 ::
1395.60
264.90
1418.90
500
25.90 1: 44.45
23.40 1: 44.41
32,10
49.08
35,00 :: 48.84
37.20 45.21
SI22
1000
43.70 1: 34.48
41.20 1: 36.02
32.00
34.13
29.10
19.68
59.10 30.34
$7.70 d: 33.68
$7.70
33.64
43.20
15.75
66.70
SI3
1000
177.10 1: 67.53
165.40
57.70
171.60
1: 66.80
175.11
d:: 69.09
195.90 1: 80.37
500
127.50 1: 40.84
132.80
40.39
129.90
+: 42.07
135.50
44.21
152.40 +: 39.94
1000
408. 1: 119.71
419.70
1: 108.01
388.50
1: 53.01
394.30
1: 88.95
395.30
109.37
500
23.10 +:: 44.78
17.70 +: 44.51
16.70 :: 44.15
20.00 :: 45.75
33.40 =:: 49.30
SF2
29.20
24.77
24.85
16.20
13.96
500
33.50 4: 32.98
25,20 +: 27.85
19,40
12.37
19.00
7.44
50.00 22.41
S13
1000
107.50
50.50
114.50
59.80
106.60
64.77
103.79
58.15
133.60 1. 88.43
500
77.70 1: 41.43
76.20 +: 41.86
67.90
26.47
79.20 :: 23.87
101.40 4: 22.37
SF4
1000
333.10 + 118.06
348.80 80 + 110.93
309.20 + $1.86
321.50
83.13
339,70 111.17
Table 2: DAG learning performance (measured in structural hamming distance, the lower the better, best results bold) of different algorithms on SF(2,3,4) graphs with different noise distributions. Our algorithms usually performs better than the previous -tate-of-the-arts DAGMA [Bello et al.) 2022).
. Order-2: (r Ar/hnv( (B IB@B B B/s) d;
Order-3: tr(f [B B /s)] d;
Order-4: tr tr(ff (B 3 B/ B/s)] = d.
In our experiments, we use the same annealing strategy for as Bello et al. (2022). During the optimization, the spectral radius of B ) B may be larger than s, which make the DAG constraints
8

---

Published as a conference paper at ICLR 2025
PC
GES
DAGMA
Exponential
Order 1
Order 2
Order 3
Order 4
SHD 563.9 d 23.84 4490.2 62.52 588.8 18.33
488.61 24.29 429.6 +: 24.73 410.6 15:25 401 + 16.64 389.4 : 16.70
Table 3: DAG learning performance (measured in structural hamming distance, the lower the better, best results in bold) of different algorithms on 1000-node ERI graphs with Gaussian noise with observation data normalized. Our algorithms performs better than the previous approaches, and as higher order DAG constraints suffers less to gradient vanishing, it tends to have better performance. We compare differential DAG learning approaches with conditional independent test based PC (Spirtes and Glymour (99T) algorithm and score based GFS (Chickering 2002) algorithm. The result reported the format of average: standard derivation gathered from 10 different simulations.
invalid. In this case, we use the strategy provided in Algorithm I to reset s. We also tried the DAG constraints (Zhang et al.] 2022) with their code provided in their appendix, but it do have some numeric issues for large scale graphs.
We compare the performance of these DAG constraints using two different settings: linear SEM with known ground truth scale and with unknown ground truth scale. We also compare these methods with constraint based PC (Spirtes and Glymour 1991) algorithm and score based combinatorial search algorithm GES (Chikering 2002) implemented by Kalainathan et al. (2020)
5.1 LINEAR SEM WITH KNOWN GROUND TRUTH SCALE
For linear SEM with a known ground truth scale, our experimental setting is similar to Bello et al. (2022): ; Zhang al. (2022); Zheng al. (0018). We generated two different types of random graphs: ER (Erdós-Rényi) and SF (Scale-Free) graphs with different numbers of expected edges. We use ERn (SFn) to denote graphs with d nodes and nd expected edges. Edge weights generated from uniform distribution over the union of two intervals 2-2 0.5) J [0.5, 2.0] are assigned to each edge to form a weighted adjacency matrix B. Then, samples are generated from the linear SEM x B x + to form an n X d data matrix X, where the noise e is iid sampled from Gaussian. Exponential, or Gumbel distribution. As Bello et al. (022);; Zhang et al. 2022): p) Zheng al. (2018) achieved nearly perfect results on small and sparse graphs we considered more challenging large and denser graphs in our experiments. We set the sample size n = 1000 and consider 3 different numbers of nodes d = 500, 1000, 2000. For each setting, we conducted 10 random simulations to obtain an average performance. All these experiments were performed using an 100 GPU, and all computations were done in double precision. Our algorithms were compared with the previous Mlate-of-the-ar approach DAGMA Bello et al. (2022). The original version of DAGMA Bello et al. (2022) used numpy and on CPU; we replaced numpy with cupy to get GPU version of DAGMA, which performed identically to the CPU version.
The results on ER2, ER3, and ER4 graphs shown in Table In all cases, our algorithms outper- formed the previous state-of-the-art DAGMA. Our Order-1 algorithm is very similar to DAGMA, except for our annealing strategy of s derived from our theory, which indicates the efficiency of our theory. Furthermore, our theory shows that higher-order constraints suffer less from gradient vanishing, and in the experimental results, we observed that the performance of higher-order DAG constraints outperformed lower-order ones in most cases. The results of SF2, SF3, and SF4 graphs are shown in Table 2 On scale-free graphs, our algorithms usually performed better than DAGMA, and the higher-order constraints, Order-2 and Order-3, often outperformed Order-1. The performance of Order-4 constraints was not good, possibly due to stronger non-convexity.
The DAGMA algorithm actually employed the same DAG constraints our Order-1 method, but with a different strategy to search for 8. Our search strategy, derived from properties of analytic functions, provides tight bound for allowing smaller to be used than DAGMA without sacrificing the numeric stability. As a result, our algorithm suffers less from gradient vanishing and achieves better performance.
In terms of running time, all algorithms had similar rumning times, typically about minutes for a 500-node graph, 10-20 minutes for 1000-node graph, and around 2 hours for 2000-node graph. Due to limited time and resources, we only considered 2000 for Gaussian noises, and for other cases, we only considered d 500, 1000.

---

Published as a conference paper at ICLR 2025
Original GRAN
Order
Order 2
Order 3
Order 4
SHD
15
13
13
SHID CPDAG
10
9
9
SHD
13
13
13
SHD-CPDAG
10
Table 4: Nonlinear DAG learning performance (measured in structural hamming distance on DAG ad CPDAG, the lower the better) of different DAG constraints on Sachs et al. (2005)' ' dataset. Different DAG constraints was plugged into the GRAN (Lachapelle al. 2020) framework. Top two rows: results obtained from single- precision mode. Bottom two rows: results obtained from double-precision mode.
5.2 LINEAR SEM WITH UNKNOWN GROUND TRUTH SCALE
For linear SEM with an unknown ground truth scale, we applied the same data generation process as for the linear SEM with known ground truth scale and Gaussian noise, but normalized the generated data X to have zero mean and umit variance. In this normalization procedure, the scale information of the variables removed from the data. Particularly for Gaussian noise, this case, the true DAG is not identifiable, and we may only identify it to Markov Equivalent Class. Previously, it has been observed that direct optimization over equation 17 may result in poor performance, mainly due the non-convex nature. In our experiments, we added an additional constraint to only allow edges exist between highly correlated nodes. first computed the Pearson correlation coefficients between every pair of nodes, and if the absolute value of the coefficient between two nodes is larger than 0.1, then we allow the edge to exist. During optimization, at every gradient descent step, we removed the disallowed edges from the candidate graph. We generated 10 instances of 1000-node ERI graphs with Gaussian noise, and 1000 observational samples were generated for each instance.
The results shown in Table 3 Our algorithm outperforms PC (Spirtes and Glymour 1991) and GES (Chickering) 2002) in terms of SHD. Although higher-order DAG constraints may suffer more from non-convexity, by adding proper constraints satisfactory results. In the results, we can see that the performance of higher-order constraints is better than that of lower-order ones, and also better than the exponential- baseed DAG constraints, which suggests that gradient vanishing may still be one important reason for poor performance.
the candidate graphs, we can still achieve
5.3 EXPERIMENTAL RESULTS ON NONLINEAR CASES
Our DAG constraints can also be extended to continuous nonlinear DAG learning approaches by replacing their original DAG constraints. We incorporated our DAG constraints into Lachapelle et al.) (2020) to model nonlinear Structural Equation Models (SEMs) and conducted experiments using Sachs et al. 22005)ss is dataset pre-processed by [achapelle al. 20020 The GraN-DAG algorithm can operate in both single-precision and double-precision modes. The experimental results are shown in Table 4 The results suggest that DAG constraints with a finite spectral radius suffer less from gradient vanishing and, consequently, from numerical truncation errors. In contrast, the original GraN-DAG algorithm experiences gradient vanishing, particularly when running in single-precision mode, as higher-order constraints that prevent long loops truncated due limited machine precision.
CONCLUSION
The continuous differentiable DAG constraints play important role in the continuous DAG learning algorithms. We show that many of these DAG constraints can be formulated using analytic functions. Several functional operators, including differentiation, summation, and multiplication, can be leveraged to create novel DAG constraints based on existing ones. Using these properties, we designed a series of DAG constraints and designed an efficient algorithm to evaluate these DAG constraints. Experiments on various settings show that our DAG constraints outperform previous state-of-the-artss approaches.
Available http:////ggtht....////aaaaaaaaaa
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
This work was partially funded by the Centre for Augmented Reasoning (CAR), the Responsible Artificial Intelligence Research Centre (RAIR), the Australian Research Council (ARC) Discovery Early Career Researcher Award (DECRA) project DE210101624 to M. Gong and DE23010159) Gong, and the Discovery Projects DP240102088 and DP240103278. We would also like to acknowledge the support from NSF Award No. 2229881, Al Institute for Societal Decision Making (AI-SDM), the National Institutes Health (NIH) under Contract RO1HL and grants from Quris AI, Florin Court Capital, and MBZUAI-WIS Joint Program.
REFERENCES
Kevin Bello, Bryon Aragam, and Pradeep Kumar Ravikumar. Dagma: Learning dags via m-matrices and log-determinant acyclicity characterization. In Advances in Neural Information Processing Systems, 2022.
David Maxwell Chickering. Learning Bayesian networks is NP-complete. In Learning from Data: Artificial Intelligence and Statistics V. Springer. 1996.
David Maxwell Chickering, Optimal structure identification with greedy search. Journal of machine learning research, 30Nov):507-554, 2002.
Max Chickering. David Heckerman, and Chris Meek. Large-sample learning of bayesian networks is up- pp-hard. Journal of Machine Learning Research, 2004.
Diego Colombo, Marloes Maathuis, Markus Kaliseh, and Thomas S Richardson. Learning high dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pages 294-321, 2012.
James Cussens. Bayesian network learning with cutting planes. In Conference on Uncertainty in Artificial Intelligence, 2011.
Nicholas J Higham. Functions of matrices: theory and computation. SIAM, 2008.
Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. in python. The Journal of Machine Learning Research, (1(1): 100::140001110, 2020.
Causal discovery toolbox: Uncovering causal relationships
Mikko Koivisto and Kismat Sood. Exact Bayesian structure discovery in Bayesian networks. Journal of Machine Learning Research, 5(Dcc):549-57) 2004.
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-hased neural DAG learning. In International Conference on Learning Representations, 2020.
Po Ling Loh and Peter Bühlmann. High-dimensional learning linear causal networks inverse covariance estimation. Journal of Machine Learning Research, 2013.
Jan R Magmus and Heinz Neudecker. The commutation matrix: some properties and applications. The annals of stutistics. (2):33113944 1979.
Jan R Magnus and Heinz Neudecker. Matrix differential calculus with applications in statistics and econometrics. John Wiley & Sons, 2019.
Ignavier Ng. AmirEmad Chassami, and Kun Zhang, On the role of sparsity and DAG constraints for learning linear DAGs. Advances in Neural Information Processing Systems, 33, 2020.
Ignavier Ng. Sébastien Lachapelle, Nan Rosemary Ke, Simon Lacoste-Julien, and Kun Zhang. On the conver- genee of continuous constrained optimization for structure learning. In International Conference on Artificial Intelligence and Statistics, pages 8176 8198. PMLR, 2022.
Ignavier Ng, Biwei Huang, and Kun Zhang. Structure learning with continuous optimization: A soher look and beyond. arXiv preprint 2r: 046, 2023.
Judea Pearl. Models, reasoning and inference. Cambridge, UK: Cambridge University Press, 19, 2000.
Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal emorr variances. Biometrika, (1):11999) 2013.
11

---

Published as a conference paper at ICLR 2025
Garvesh Raskutti and Caroline Uhler. Learning directed acyclic graph models based on sparsest permutations. Stat, 111, 2018.
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-27766 2021.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein protein-signaline signaling networks derived from multiparameter single-cell data. Science, 308(5721):523 529, 2005.
Rajen Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure, The Annals of Statistics, (8(3):1514-155) 2020.
Shohei Shimizu, Patrik O Hoyer. Aapo Hyvärinen. and Antti Kerminen, linear non Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(Oct):2003-2) 2030, 2006.
Tomi Silander and Petri Myllymäki. A simple approach for finding the globally optimal Bayesian network structure. In Conference on Uncertainty in Artificial Intelligence, 2006.
Ajit Singh and Andrew W. Moore. Finding optimal Bayesian networks by dynamic programming. Technical report, Carnegie Mellon University, 2005.
Daniel A. Spielman. Lecture notes in spectral graph theory: The adjacency matrix and nth eigenvalue, Septem- ber 2012. http::ttpps:.......ee CB yale. edu/ wwwphomw..eeeeeeeeeeeeeee hwoome//speelii www.////eeemeeeeeeeeeeeeeeaaeees
Peter Spirtes and Clark Glymour. algorithm fast recovery of sparse causal graphs. Social science computer review, (1):6227 1991.
Peter Spirtes, Chris Meek, and Thomas Richardson. Causal inference in the presence latent variables and selection bias. In Conference on Uncertainty in Artificial Intelligence, 1995.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Matthew Vowels, Necati Cihan Camgoz, and Richard Bowden, D'ya like dags? a survey on structure learning and causal discovery. arXiv preprini arXiv:2103.0255: 2021.
Dennis Wei, Tian Gao, Yue Yu, Dags with No Fears: A closer look continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 33:3895-3906,, 2020.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: Dag structure learning with graph neural networks. In International Conference on Machine Learning, pages 71 7154-7163. PMLR, 2019.
Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji. DAGs with no curl: An efficient dag structure learming approach, In Proceedings of the 38th International Conference on Machine Learning, 2021.
Changhe Yuan and Brandon Malone, Learning optimal Bayesian networks: A shortest path perspective. Journal of Artificial Intelligence Research, 48(1):23365, 2013.
Changhe Yuan, Brandon Malone, and Xiaojian Wu. Learning optimal Bayesian networks using A search. In International Joint Conference on Artificial Intelligence, 2011.
Bin Zhang, Chris Gaiteri, Liviu-Gabriel Bodea, Zhi Wang. Joshua McElwee, Alexei Podtelezhnikov, Chun- sheng Zhang, Tao Xie, Linh Tran, Radu Dobrin, et al. Integrated systems approach identifies genetic nodes and networks in late-onset alzheimer's disease. Cell, 1530 2013.
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, 15:1829--8,40,, 2022.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. DAGs with No Tears: Continuous optimization for structure leamning Advances in neural information processing systems, 31, 2018.
12

---

Published as a conference paper at ICLR 2025
Appendices
A PROOF OF CONVERGENCE OF MATRIX ANALYTIC FUNCTION
To prove that fAA converges for any matrix A with spectral radius p(A) < r, we begin by considering the power series representation of the analytic function f(z). Since fzz converges for ax| < r, can be expressed as:
f(x)
c
A==
where the coefficients C₁ are such that the series converges absolutely for |x| < т. Now, A be an matrix with spectral radius p(A) < By the definition of the spectral radius, all eigenvalues A satisfy |^| < r.
We define the matrix function f(A) by substituting A into the power series:
S(A) Σ c
k=
To show that this series converges, sub-multiplicatio matrix norms of A. Specifically, for any e 0, there exists sub-multiplicative matrix norm such that ||||| S p(A) t.. Since p(A) < r. we can choose small enough that
the fact that the spectral radius p(A) the infimum of all
|||]| < r.
To justify this claim, we provide proof that the spectral radius p(A) is indeed the infimum of all sub- multiplicative norms of A. Let || || be any sub- cub-multiplicative matrix norm. For any eigenvalue À of A with corresponding eigenvector v, we have:
|||||||| || || s ||||||all
Since ||||| * 0, it follows that |\| < ||A|| Taking the supremum over all eigenvalues we obtain p(A) < A|] This shows that p(A) is a lower bound for all sub-multiplicative norms of A.
To show that p(A) is the infimum, we use the fact that for any >> there exists a sub-multiplicative ) + c. This can be achieved, for example, by using the Jordan canonical form of A and constructing an appropriate norm. Thus, p(A) is the infimum of all sub-multillicative norms of A.
norm |. |le such that || S p(A)
Now, consider the norm of the series:
DO Σ
S Σ
lc |I
Since
||||| and the original series kk=0 Car converges absolutely for |x| < k= Cx || also converges. Therefore, the matrix series
the series E Aa converges absolutely, and
((A)
dell-defined.
This completes the proof that f(A) converges for any matrix A with spectral radius p(A)
B PROOF OF PROPOSITIONS
Our proof are also based on the well-known propertics of analytic functions listed as follows:
1. Let fi(x), (2() be analytic functions on (-11,71) and (-r2, r2), then fi(x) + f2(x) and fi (2() are analytic functions on (- mim((1i, (2), min(F1.722)
2. Let f(x be analytic function (-r,r), then (f(x))OOr is analytic function (-r,r).
B. LEMMAS REQUIRED FOR PROOFS
Lemma 9. Let E RRdd be the weighted adjacency matrix of a graph with d vertices, G is
DAG and only B
13

---

Published as a conference paper at ICLR 2025
Proof. See Proposition 3.1 of Zhang al. (2022).
Lemma 10. Let B I220 be the weighted adjacency matrix of a graph G with vertices, 9 is DAG if and only
d
C3
0,
where 6 Ovi.
Proof. See Wei al. 2020).
B.2 PROOF OF PROPOSITION 0
Proposition 1 Let B B with p(B) be the weighted adjacency matrix of a directed graph G, and let f be analytic function in the form of equation 3 where we further assume Vi >> 0 we have 64 > 0, then G is acyclic if and only if
[(B)] cod.
Proof. Without loss of generality, assume that f can be formulated as:
f(x) = co + Σ Vi,c 0; lim c//4+11 > O.
(22)
60000
First if G is acyclic, by Lemma we must have
Ovk >
(23)
which also indicates that p(B) = 0. Thus we have
tr [(B)] =tr
co +
c
Σ
cBB'
i=1
d=
Equals 0, By
(24)
=tr[col] +
tr
B
Equals 0, By Lemma 10
=cd.
On the other hand, if tr [[B) = end, we must have that
tr
c
=
O.
By the fact all entries of are positive, we have that
0 < tr
B
=
0.
(25)
Then we must have
tr
Σ
B
0.
Finally by Lemma 10 we have that 9 is DAG.
14

---

Published as a conference paper at ICLR 2025
B.3 PROOF OF PROPOSITION
In all the paper, we consider analytic functions f from the functional class F defined in equation S Proposition 2 There exists some real number where for all {B E Rdxd (p < r. the derivative of [B) w.r.t. B
VBir [[B]] =
Proof. Without loss of generality, assume that can be formulated as:
f(x) co Σ cart 0; lim c//441 > O.
(26)
21000
For some i by basic matrix differentiation we have
dtr
(iBI
(27)
and then by the properties of power series we have
VBt
[(B)
-V
Co
I
CB
trre
(28)
x
Ci
cai
(vf((x)
where we can exchange v and Σ i=1 because after the exchanging the new power series will still converge (by properties of analytic functions).
B.4 PROOF PROPOSITION 3
Let f(z) = + Σ qrt E F be analytic function on (-r,r), and let be arbitary integer larger than then B Rdxd with spectral radius p(B) r forms DAG if and only if
Proposition 3
0" (xx
n.cn.
Ox
Proof. By properties of analytic functions, the th order derivative of an analytic function f(x on (-r,r) is still an analytic function on (-r,r). Particularly for fxx co Ei ca2 E F, we have
8r f(x)
0 ca
Oxn
Oxn
(29)
Gar
k
i=n
k=
1+1
12+
=nlc,
+
Ci+m x
15

---

Published as a conference paper at ICLR 2025
where by the fact > 0, Vi > 1, we have that proved the proposition.
8" f(t) Or
F. Then by Proposition
|
immediately
B.5 PROOF OF PROPOSITION
Proposition
Let f1(x) f + cf EF, and f2(x) - Got E Then for ad-
jancency matrix B R with spectral radius ((B) & min((imm->so c//c+1 / c//22+1))
the following three statements are equivalent:
1. B forms DAG;
2. tr(fi(B) + f2(B)] = (cb + cG)d:
3.
= cdr
Proof. By properties of analytic functions, we have
f1( + f2(x) = dj + Σ (c) + c))x
(30)
1=1
is an analytic function, and its convergence radius is given by lim (c) + (c) c.+1) min 1 21000
c//c/+1 lim 22/c2t1), i+:00
(31)
and thus by Proposition | the statement and are equivalent. Similarly by properties of analytic functions statement and 3 are equivalent. Thus the 3 statements are equivalent.
B.6 PROOF OF PROPOSITION 5
Proposition
Let n be any positive integer, the adjancency matrix B [B R R IP s)
and only if
tr trs)-
and the gradients of the DAG constraints satisfies that vB E {B E Bood pBB < 5}
VVgir(I - B/s) (Bs))
"Il < Brr - B/1
where l an arbitrary positive integer, and || || denote an arbitrary matrix norm induced by vector p-norm.
Proof. By Proposition 4 or Proposition 3 it would be straightforward that tr(I - B) 3 = d is
necessary and sufficient condition for an adjacency matrix B DAG.
{
E R
(( & $} to form
For the norm gradients, it straightforward that
0(1 x))
n21 x)) n-1
(32)
dx
For arbitrary we have
n-
(1 -1- (1-2)--
1+ 1=1
x
(33)
and obviously the coefficients is monotonic increasing w.r.t. n. Bdxdlp(B) 8) we have for any > 0, 0
Thus by the fact VB E
{ e
[I - B) < | B
-
(34)
As result, we have
V -I =I
B) - < (n+ B B)- -- '11
(35)
(n + k)]l - - Vgtr(I B) -n
16

---

Published as a conference paper at ICLR 2025
B.7 PROOF OF PROPOSITION 6
Proposition 6 Thus the Hessian of the DAG constraints equation 6
can be oblained as follows:
i-2
ttr = Kddd
ic
B
0
B
where & denotes the Kronecker product.
Proof. Firstly, the derivative of matrix power can be obtained using the following equation (Magnus
and Neudecker 2019).
2
tr
B
B
(36)
j=0
where denotes the Kronecker product. Thus the Hessian of analytic DAG constraints can be obtained as follows:
V
LrB(BB = Σ
Glr 2
B
B
1=
4-2
B 1-2-1
K
ic
(37)
2
j=0
B.8 PROOF OF PROPOSITION 7
Proposition 7 For two analytic function fi(x) = co.1 + E G,1x² E F and fa(œ) = C0.2 + 1 we have Ci Ct2, then
E 22,20 E FF
(Vattf > trf2
where p(-) denotes the spectral radius of matrix.
Proof. Obviously, each entries the Hessian of Irfi (B)) is larger than the corresponding ones in tr tB)). Thus for any unit length vector u with all positive entries we would have
V tr)i u tarfa((
and then it would be straightforward that
p (rB((
u T V tr (B)) u
max
muu30 1|==1
max u
u
V (B)) u ((
C ERROR ANALYSIS OF THE ANALYTIC FUNCTION BASED DAG CONSTRAINTS
Without loss of generalization, we analyze the truncated error for a DAG constraints
f(B) col + Σ B
which converges for p( < 1, and the results can be easily generalized to normal analytic function based DAG constraints.
17

---

Published as a conference paper at ICLR 2025
it would be obvious for truncated approximation
[x - col +
Σ G
==1
the error between the truncated approximation and the original analytic function can be obtained by
[ - (B)
BBI12
k+1
B 2
Let
sup
BBp(B)<1
B
where we define C(f. k)) (1)) EE= 12. It notable that |B*+1|| converges to exponentially w.r.t the increase of k due the fact (BB) 1. Meanwhile, it obvious that C(J,k) increases at most in polynomial rate due to that fact Di= è, 2,B' converges for ((B) 1. Finally, the residual of Algorithm 2 converges to 0 in rate of a polynomial of divided by the exponential of k.
D HYPER PARAMETERS
In terms of hyper-yarameeneeee parameters, our selection involves a = 0.1, A1 = 0.1, and T = For we use the same annealing approach Bello al. (2022), but with strategy reset when candidate graph goes out of the desired region.
In all experiments in this paper, for continuous based approaches parameter as [ello et al.] (2022). for conditional independent test and score based approaches we use the default parameter in Causal Discovery Toolbox
exactly the same hyper
E EXTRA EXPERIMENTAL RESULTS
In this section, provide additional experimental results, including true positive rate, false detection rate and running time for large scale graphs, as well as experimental results on small scale graphs.
1 scoveryToo aooverryloolbbo/httt /index.
html
18

---

Published as a conference paper at ICLR 2025
Table 5: DAG learning performance (measured in true positive rate, the higher the better, best results in bold) of different algorithms on large scale (500-2000 nodes) graphs with different noise distributions. Our algorithm performs better than previous approaches.
Graphs
Nodes
DAGMA
Order-1
Order-2
Order-3
Order-4
FR2-gaus
0.97 0.01
0.98
0.01
0.98
0.01
0.98
0.01
=:
ER2 guuss
1000
0.97 1: 0.01
0.98 1: 0.01
0.98
0.00
0.98
1: 0.00
0.98 0.00
ER2 gauss
2000
0.96
0.01
0.97 1 0.01
0.98
0.00
0.98
0.00
0.98 0.00
FR1-gauss
0.95
ERI 3R3ss
1000
0.94
0.01
0.96
0.01
0.96
0.01
0.97
0.01
0.97
0.01
ERI-gane
2000
0.94 1: 0.01
0.95 - 0.01
0.95
0.01
0.96 2: 0.01
0.96 2: 0.00
FR4-gan4s
500
0.92 - 0.02
0.94 : 0.02
0.94
0.02
0.94 =: 0.02
0.95
0.02
ER1
ER+4-guss
2000
0.89
0.01
0.90
0.01
0
0.01
0.91
0.01
0.91 : 0.01
ER2
500
0.96 H: 0.01
0.97
+ 0.01
0.98
0.01
0.98
0,01
0.98
# 0.01
BR2
1000
0.97 1: 0.01
0,98
0.00
0.98
0,00
0,98
0.00
0.98 0.00
2000
LRR exp
500
0.95
0.01
0.96
0.01
0.97
0.01
0.97
0.01
0.97 0.01
ER3
1000
0.94 1: 0.01
0.96
1: 0.01
0.96
0.01
0.97
1 0.01
0.97 1: 0.01
ERI-exp 2000
0.00 + 0.00
0.09
0.28
0.09
0.28
0,10
0.29
0.00 : 0.00
0.94
ERI-exp
1000
0.89 1: 0.01
0.92
0.01
0.92 + 0.01
0.93
0.01
0.94
1: 0.01
134 exp
2000
0.00 + 0.00
0.09 + 0.27
0.09
0.27
0.09
0,27
0.00 + 0.00
ER2 gumbel
0.95
0.99
0.01
0.01
0.99
0.01
0.99
0.00
LIK2 gumbel
2000
0.00
0.00
0.59
0.48
0.10
0.29
0.10
0,30
0.10 + 0.30
ER3 gumbel
500
0.95 1: 0.01
0.97
1: 0.01
0.98
0.01
0.98
1: 0,01
0.98
1: 0.01
ER3-g 1RI-gubbll
1000
0.95 1: 0.01
0,97
0.01
0.97
0.01
0.97
0.01
0.98 0.01
gumbel
2000
0.00
ER4 gumbel
500
0.93
0.02
0.95
0.01
0.96
0.01
0.96
0.01
0.96
0,02
ERl-gumber
1000
0,90
1: 0.01
0.94
1 0.01
0.94
1 0.01
0.95 - 0.01
0.95
1: 0.01
EX44 gumbel
2000
0.00
0,00
0,00
0.28
SF2-gmos
500
0.97 + 0.01
0.98 + 0.01
0:98
0.01
0.98 = 0.01
0.98 + 0.01
SF2-gones
SF2-gmms
2000
0.96
0.01
0.97
0.01
0.98
0.00
0.98
0.00
0.98
0.00
SF3-gmes
500
0.95 d 0.01
0.96 + 0.01
0.97
0.01
0.97
0.01
0.97
= 0.01
SF3 SF3-gmos
1000
0.94 2 0.01
0,96
0.01
0.96
0.01
0.97
0.01
0.97
0.01
2000
0.00
SF4-gams
0.92 0/02
0.94
0.94
0/02
0.94
0.02
SF4-gaoss
1000
0.89 1: 0.02
0.92 + 0.01
0.92
= 0.01
0.93
- 0.01
0,93
1: 0.01
SF4-gmos
2000
0.89 2 0.01
0.90
0.01
0:90
0.01
0.91
0.01
0,91
0.01
SF2-exp
1000
0.97
0.01
0.96
0.00
0.98
0.00
0.98
0.00
0.98
0.00
SF3-esp
500
0.95
2: 0.01
0.96
- 0.01
0.97
0.01
0.97
= 0.01
0.97 1 0.01
SF3-exp
1000
0.94
1: 0.01
0,96
0.01
0.96
0.01
0.97
0.01
0.97
0.01
SF4-esp
1000
0.89
0M01
0.92
8.01
0.92
0.01
0.93
= 6.01
0.94 -
8r2-gumbel
500
0.98
1: 0.01
0.99 1: 0.01
0.99
1 0.01
0.99
1 0.01
0.99
1: 0,00
8F2-gumbel
1000
0.98 2: 0.01
0.99
2 0.01
0.99
0.00
0.99
0,00
0.99
2 0.00
SFi-gunhel
303-pumbel
1000
0.95
0.01
0.97
1 0.01
0.97
0.01
0.97
0.01
0.98
0.01
SF4-gumbel
500
0.93
1: 0.02
0.95
- 0.01
0.96
2 0.01
0.96
2
0.01
0.96
2: 0.02
SF4 gumbel
1000
0.90
0,94
0.01
0.94
0.01
0.95
0.01
0.95
0.01
19

---

Published as a conference paper at ICLR 2025
Table 6: DAG learning performance (measured in false detection rate, the lower the better, best results in bold) of different algorithms on large scale (500-2000 nodes) graphs with different noise distributions. Our algorithm performs better than previous approaches.
Graphs
DAGMA
FR2-gauss
0.02
0.02
0.01
ER2 guuss
1000
0.02 0.01
0.02
0.01
0.01
0.01
0.01
0.01
0.01
0.01
ER2
2000
0.02
1: 0.01
0.02 + 0.01
0.02
0.01
0.02
2: 0.01
0.02
1 0,00
FR FR1-ganus
500
0.04 - 0.02
0.04 : 0.02
0.03
0.02
0.03
- 0.01
0.04
0.02
ER1 gauss
ERl-gaus
2000
0.06
0.01
0.06 - 0.02
0.05
0.01
0.04
0.01
ERI gauss
500
0.08
6: 0.03
0,08 +: 0,03
0.08
0.04
0.07
:: 0,04
0.07
:: 0.03
ERI
1000
0.12
0.03
0.11 +: 0,03
0,10
0.02
0.09
0,03
0,10 0.03
BR4-gauss
2000
1122 exp
500
0.03
0.02
0.02
0.01
0.01
0,01
0.01
0.02
0.02
1: 0.02
ER2-aNp
1000
0.02 1: 0.01
0.02
+ 0.01
0.01
0.01
0.01
:: 0,01
0.01
:: 0.01
ERI-eNp
500
0.05 +: 0.02
0.04
0.02
0,03
0.01
0.04
0.02
0.04
0.02
ER4-exp
0.09 0.04
0.09
0.04
1 0,05
0.07
0.04
0.06
1: 0.03
134 exp
1000
0.12 + 0.003
0.11 + 0,02
0.11
0.012
0.09
0,03
0.10 # 0.02
BR2 gumbel
0.03
0,01
0,01
0,01
0.01
0.01
L3k1 gumbel
500
0.06
0.002
0.04
0.02
0J03
0,03
0.03
0.01
0.04
0.02
BR3- gumbel
1000
0.06 à: 0.01
0,05 1: 0.01
0.04
1 0.02
0.04
1: 0.01
0.03
1: 0.01
ERi-gumbv
500
0,10
0.04
0.09
0,05
0,09
0.04
0,08
0,04
0.08 1: 0.05
184 gambel
1000
SF2-amss
0.02
SF2-gamss
0:02
0.02
0.01
SF2-gmms
2000
0.02 + 0.01
0.02 + 0.01
0.02
= 0.01
0.02 = 0.01
0.02 1 0.00
SF3-gmos
500
0.04 + 0,02
0.04 + 0,002
0.03
0,02
0.03
0,01
0.04
0.02
SF3-pmes
SF3-gmss
2000
0.06
0.01
0.06
- 0.02
0.05
0.01
0.04
0.01
SF4-gmmss
500
0.08
d: 0.03
0.08
+ 0.03
D
0.004
0.07
0.04
0J07 = 0,03
SF
1000
0.12
0.03
0.11
0.03
0.10
0.02
0.09
0,03
0.10
0.03
2000
SF2-cxp
0.02
0.02
0.01
0.01
0.02
SF2-cxp
1000
0.02
1: 0.01
0.02
+ 0.01
0J11
= 0.01
0.01
0.01
0.01
2: 0.01
SF3-csp
500
6.05
: 0.92
0.04
0,02
6.03
0.01
0.04
0J02
0.04
0.02
1000
SF4-exp
0.09
0.04
0.09
0.04
0.09
0.05
0.07
0.04
0
0
SF4-esp
1000
0.12
0.03
0.11 - 0.02
0111
0.02
0,09
0,03
0.10
0.02
SF2 gumbel
500
0.03
0.02
0,01
0,01
0,01
0.01
0,01
0.01
0,01
0.01
SF2-gumbel
1000
SFi-gunhel
0.06
0.02
0.04
0.03
0.03
0.03
0.01
0.02
SF3 sumbel
1000
0.06
1: 0.01
0.05
- 0.01
0.04
0.02
0.04
2 0.01
0 0.01
SF4-guneel
500
0.10
0.04
0.09
0.05
0.09
0.04
0,88
0,04
0.08 0.05
SF4-gunhel
1000
0.14
0.12
20

---

Published as a conference paper at ICLR 2025
Table 7: DAG learning performance (measured in running time (seconds), the lower the better, best results in bold) of different algorithms on large scale (500-2000 nodes) graphs with different noise distributions. Our algorithm performs better than previous approaches.
Graphs
DAGMA
Order-1
Order-2
Order3
Order-4
FRI-guus
204.19
62.47
403.73 + 78.79
482.01
17.50
512.86 +: 33.87
ER2 gauss
1000
364.95
38.47
687.38
69.45
1261.33
47.66
1256.89
33.76
1329.92 1: 315.91
ER2
2000
1187.11
1 108.06
3515.56 1 310.54
6063.65
2 554.57
6300199
1 197.69
6360.65 1 521.65
FR FR1-gmns
500
238.54
58.34
394.95
= 110.7
516.49 + 14.45
524.93
+ 48.79
528.08 + 49.59 50
1: 162.25
21,10
1355
ERI-guuss
2000
1636.23
101.64
4903.56
+: 782.06
6037.67
765.47
7072.04
2168.01
6969.54 : 724.52
ERI gauss
500
347.06
55.40
519.77
+: 56.88
532.34 +: 10.16
$17.46 +: 12.95
$34,98 +: 52.22
ERI ERI-gauss
1000
798.01
27.84
1328.52
1: 48.32
1318.45
62.71
1348.08
15.79
1312,53 1: 138.14
ER4-guus
6520.96
2151.48
ER2 exp
500
167.08
25.75
291.87 63,06
492.14 +: 30,89
196.83 20.55
504,17 1: 22.59
ER2-exp
1000
359.97
1: 28.62
708.23 1: 61.00
1245.71 1: 71.4
1279 58 1: 41.63
1277. 1 45.39
ER3-exp
500
235.71
56.44
440,70 = 189.91
564.88 145.09
$59.54 144.42
550.44 + 179.03
ER3 exp
1000
1100.15
1371.65
1401
321.35
ERi-exp
500
358.94
46.65
510.53 +: 46.91
513.70
19.38
$22.16
15.09
308.01 1: 34.87
ER4
1000
778.40
51.69
1344,07 +: 26.71
1324.00 + 111.40
1347.33 + 21.06
1298.91 + 100,30
BR2 gumbel
24.90
255,86 1: 33.56
501.55
85.14
11.06
301.44 1: 23.26
ER2-
1276.58
15.41
ER3 gumbel
300
232.45
50.71
381.033
85
521.98
8.54
514.79 12.92
206.95 + 18.01
ER3 gumbel
1000
525.92
1: 93.20
1013.67
1: 168.84
1331.88 1: 99.25
1302.90 1: 10L.5
1369.44 44 1: 39.47
ER ERM-gambel
500
366.06
33,40
514.95 + 57.57
$40.93 : 17.35
$30.18 1: 20.49
519.87 1: 13.95
ER4 gumbel
1000
1434.41 154.30
1335.20
SF2-gmss
294.19
493.73 78.79
17.59
$12.86
SF2-gaess
1000
687.38 69.45
1261.33 47.66
1256,39 33.76
1329.92 + 315.91
SF2-gmss
2000
1187.17 :: 108.06
3515.56 + 310.54
6063.65 = 554.57
6300 1 197.6
6360.65 1 521 65
SF3-gmss
500
238.54 + 58.34
394.95 + 11087
516,49 + 14.45
124 + 48.79
$28.08 + 49 50
SF3-gomss
1355
1349
SF3-gess
2000
1636.23
101.64
4903.56 783.06
6057.67
765.47
7072:044 2168.01
6969.54 +
SF4-2mss
500
347.06
55.40
$19.77 - 56.38
$32.34
10.16
517.46 12.95
514.98 + 52.22
SF4gomos
1000
798.01
27.84
1328.52 1: 48.32
62.71 1318.45
1348 15.79
1312,53 + 138.
2000
311.23
SF2-esp
25.75
496.83 20.55
504.17 +
SF2-esp
1000
359.97
1: 78.62
708.23 1: 61.00
1245.71 71.3
1279 38 41.63
1277 30 1 45.39
SF3-esp
500
235.71
1: 56.44
440.70 18999
564.88 +: 145.09
559.54 144.42
550.44 + 129.01
1000
226-9
321.35
SF4-exp
358.94
46.65
510.53 + 46.9
513.70 2 49.38
$22.16 15.09
508.01 + 34
SF4-esp
1000
778.40
1: 51.69
1344.07 + 26.1
1354.00
111.40
1347 26.06
1298.91 + 100.00
SF2-gumbel
161.36
24.90
255,86 1: 33.56
501.55
85.14
490.99 11.06
501.44 1: 23.26
SF2-gumbel
1000
1276
SFi-gunhel
50.71
381.03 67
8.54
514.79
506.95 +
SF3-gumbo
1000
525.92 1 93.20
1013.67 1: 168.84
1331.88 99.25
1302.90 101.85
1369.44 1: 39.47
SF4-gumbel
500
366.86
33.40
514.95 = 57.57
$40.93 : 17.35
530. 18 30.49
519.87 + 13.98
SF4-gumhll
1367.02
1360.93
1454-41
21

---

Published as a conference paper at ICLR 2025
Table 8: DAG learning performance (measured in structural hamming distance, the lower the better, best results in bold) of different algorithms on small scale (10-100 nodes) [2,3,4,4 graphs with different noise distributions. Our algorithm performs better than previous approaches.
Grupts
Nodes
MMPC
GES
NOTEARS
DAGMA
Order-1
Order-
Order-3
Order-4
ER2-goms
10
21:01 +: 3.77
14.66 + 9.25
301 2.7
2.01 +: 1.91
0.74
1.22
0.93 +: 1.49
0.97 1.65
1.27 +: 2.11
FR2-gooms
30
65.23
+ 9.26
$7.95 = 44.65
600 1
3.76 - 3.70
1.88 2 2.97
1.57 + 2.20
1.49 2 2.32
1.87 # 2,89
EK2-gamse
50
10K.74 + 11.35
90.49 = 58.75
11.59 + 10:25
4.21 + 3.66
2.28 3.06
2/05 +: 2.96
2:21 + 2.98
2.43 + 3.69
FR2-gomss
155.61
22.58 17.69
7.45
364
4.54
FR3-ga
$2.04 + 3.39
26.39 = 6.85
6.6 -: 3.69
1:6:=:2. 2.91
1,35 + 2.97
3:1
3.51 + 2.83
ER3 3R3-g
30
97.56 1 9.71
158.57 1: 49.92
17.89 12.29
9:22 1: 6.24
4,83 1: 4.56
5.33 5.52
S.54 5.71
4.96 1: 5.94
FR3 TR3-goms
50
160.21
+ 12.67
265.88 96
30,70 18.48
13.05 = 8.68
6.60
6.56
6,43 4: 8.72
6.09
6.44 4: 5.54
64-27
11.59
+: 9,92
10.66-99.02
11.78
ER4-pons
10
40.37 1: 2.24
29:07 5.77
13:71 1: 3.93
11.57 3.27
6.89 3.39
6,75 1: 3.15
7.08 3:40
7.04
FR4-goms
10
126
= 10,70
222.49 = 38,99
30.18 1: 21.79
20.67
9.26
11,12 7.76
10.54 1 7.36
12.58
8.81
12.81
= 10.74
FR4-gums
50
210.43 + 13.15
495,14 81.60
60.11 25.25
34
12.48
14.39
9.72
12,89 +: 9.53
15 + 10.31
16.22
+ 10.17
ER4 gans
1047.05 250.05
51.53
28.31
1:23.60
23,03
FR2-exp
28 + 3.70
15.23 = 9.49
2.64
2.19 - 1.87
0.80
1.34
0,90 + 1.57
1.74
1.35 + 7.30
ER2 exp
30
61.98 1: 9.21
39.37 1: d6.71
6.59 1: 7.08
3.83 1 4,10
2.22 1: 3.63
1,85 1: 2.91
1 1: 2.37
2.09 1: 3.28
ER2-exp
50
106.77 + 10.85
95:43 56
11.51 9.70
4.21 + 3:43
2.12
2.75
1,97 1: 2447
2662 3.01
2.80 4: 3.53
FR2-esp
215.70
17.48
+
1:54::44
5.15
ER3-exp
10
31 90 1: 3:41
26,29 6.80
9:32 4.42
6.98
1: 3.43
3.57 1: 2.91
3.43 1: 2,96
3.10 2.89
3,68 1: 2.82
ER3-esp
30
97. 1: 9.81
154.20 = 48.45
17.54 1: 11.51
9.11 + 5,51
5.07
2: 4,47
5.51 1 5.51
5.60 2 5:73
4.81 :: 5.76
FR3-exp
50
160 16 + 12.49
268.43 99.00
31.32
19.96
14.16 9.75
7.92
7.38
6,73 4 5.94
6.32
7.18 + 5.69
BR3-esp
494.33
30.84
12.26
11.52 9:337
ER4-exp
10
40.44 + 2.15
28.79 + 5.56
13.92 3.89
11,85
3,42
7.01
3,44
6,93 +: 3.36
7.07
+ 3.55
7.10 1: 3.32
BR4 cxp
30
125 98 + 10.62
221.86
d: 39.40
36.55
18.74
21,00
10.26
10.73
7,83
11.55 + 7.50
13.19 +: 9:95
12.87
:: 11.09
ER4-esp
50
210.42 +: 13:45
194,91
77.20
$0.90
25.93
27,10
12.94
14.76
9.50
13.83 1: 9.46
14.92 1: 10.98
16.44 1: 12.95
FR4-esp
230
48.63
27,02
24.07
22,89
22.40
ER2 gumbel
10
20.39 4: 3.72
16.07
8,92
1.69
2.04
1.11 + 1,52
0.67
1.36
0.82 1: 1,50
0.79 1.65
1.00 :: 2.24
ER2-gumbel
30
64.88 1: 9.55
61.33 1 43.70
5.65 1: 7.35
2.10 + 2.92
1.76 2.79
169 1: 2,90
2.13 1: 3.06
2.28 1: 3.74
ER2-gumbe)
50
106 + 10,76
92.35 57.81
10.24 10.44
2.81 + 3,30
2.43 3,03
2.64 1: 3.8
3.1 4:43
3,77 1 4.56
gumbel
215.52
161.75 89.91
22.39 19.17
4.89
4.78
6.00
ER3 gumbe]
10
31.95 1: 3.37
26,29
6.69
5:93 3.95
4.05 1: 2.96
1.90 1: 2.39
2.44 1: 2.71
2.09 2:31
1.75 1: 1.23
ER3 gumbel
30
97.72 + 9.72
158.60
29.62
13:44
11.65
4.97 + 4.07
3.72 4,11
4.13 + 4.29
5,99 + 6.27
5.96 + 6.92
ER3 gumbel
50
160.35 :: 12.66
281,98
102.11
27.28
20.76
9.43 1: 7,56
1.29 7,07
5.79 1: 5.76
7.89 6.610
8.17 8,16
ER3 gumbel
321.87
60.12
29.08
13.15
12.62
1:
12.45
9:99
ER4 gumbel
10
40.48 + 2.13
28.69
5,93
10.24 4.08
7.75 + 2.94
3,96
3.11
4.55 + 3:54
4.92 4-3.8
4.50 + 3.46
ERI gumbel
30
126.04 +: 10.54
219.39
:: 43.15
29.45
1: 19.08
12.35
= 7.64
8.37
8: 6.70
8.14 1: 7.53
11:29 1: 9.49
11.84
1: 9.91
ER4-zumbe)
50
210.16
:: 13.34
499.02
75.43
55.57
28.67
18.03 15:20
12.07
1: 9.36
10.92 1: 9M4
13:44 1: 10.94
15.06
:: 11.37
ER4 gumbel
100
20.06 424.19
1031.32
243.28
114.40
52.86
22.69
27,00
28:51
22.36
21.24 23.58
25.37
25.67
SP2-gunss
14.82 1.99
5.01 6.58
2.67
2.01 + 1.91
1,22
0:93 1.49
0.97 1.65
1.27 2.11
SF2-guss
$4.01 + 3.54
19:50
6:31
3.26 + 3.70
2.97
1.57 + 2.29
1.49
1.87 + 2,89
SP2-gauss
50
97 30 5.13
43.39 + 40
11:59 +: 10.25
4.21 + 3.66
2.28 2: 3.06
2,05 1: 2.96
2:21 + 2.98
2.43 1: 1.69
SF2-guiss
100
215
= 8.24
106.51 62.29
22.58 2 17.69
7:45 + 6.70
3.49 3.97
3.64 + 4.48
4.09 4.40
4.11
+ 4.54
SF3-gams
+: 2:77
7.78
1.433
2,91
3.35
2.83
SFi-ganss
30
65
7.29
31.51
32.61
17.09 + 12.29
9.32 2 6.24
4.33
4.56
5,33 +
044::577 5.71
4.96 + 5.94
SFi-guiss
50
119.92 =: 10.81
69,43 58.72
50.70 2: 18.48
13.05 = 8.68
6.60 2: 6.56
6.43 + 5.72
6.55 - 6.09
6.44 + 5.54
SF3-gaus
100
271 55 + 16.84
157,11 = 99.22
64.27 33.28
21.66
15,48
11.59
9.90
10.69
1: 9.92
10.86:1:90
11.28 9,99
SF4-ganss
+ 3.20
13.71
7
SF4-gusss
30
71 + 9.41
46.72
30.18 21.79
9.26
11.12
7.76
10.54 + 7.35
12.58
12.81 +
SRi-gass
50
137.4 = 13.92
90.76 = 69.01
60.61 1 25.25
24.96
12.48
14.39
9.72
12.89 1: 9.52
15.15 1 10.31
16.22 : 10.17
SF4-guiss
100
$15. 31 + 74.17
174.28 = 109 72
118 16 51.53
42.81
23.74
28.31
25.58
24.47 4 24.31
23.22 + 23.60
23,03
= 21.51
SF2-exp
1.88
2.64
2.19
0190 + 1.57
1.35
2.30
SP2-exp
30
54
3.69
23.86
24.42
6.59 7.08
3 1: 4.10
2.22
3.63
1,85 2:1
1
2-37
2.09 + 3.28
SP2-oxp
50
97 + 5.42
47.29 = 42.50
11.51 1: 9.70
4.1 - 3.43
2.12 2 2.75
1.97 1 2,47
2,6 2 3.01
2,80 : 3.53
ST22-ep
100
216
7,98
113.30 63.41
21.87 17.48
7.31 1 6.17
3.71
4.77
4.11 4.558
4.54 1:441:477 4.72
4.82 1 5.15
SF3-esp
698
1.82
SF3-cpp
30
65
7.38
32.67
36.15
11.51
9.11
- 5.51
5/07
5:51 + 5:51
- 7.15 6.332
4.11 5.76
SF3--xp
50
119.86 1: 10.63
64.53 + 55.26
31.32 1: 19.96
14.16
9.73
7.92 1: 7.38
6.73 1: 5.94
7.18 1: 5.69
SF3-exp
100
272. +: 16.44
157.53
109 02
66.39 2: 30.84
22.02
15,03
12.26
9.06
12.44
+ 11:03
1..3
11.14 8,36
SF4-2xp
+
13.92 1:
3.55
SRi-exp
30
73.49
1: 8.97
47.03
40.43
36.55
18.74
21,00
10.26
10.73
7,83
11.55 1 7.70
13.19 9:95
12.87 1: 11.09
SRi-exp
50
137.96 +: 14.04
84.43
63.15
60.90
25.93
20
12.94
14.76
=: 9.50
13,33 1 9.46
14:92
+ 10.98
16.44
= 12.95
SR4-esp
100
314.36 1. 24.83
176.43
114.42
116.91 18.63
12.64
23.78
27,02
24.07
24.32 1. 22.06
22,89 23.39
24.12
1: 22.60
ST22-gumbel
1.91
6.51
L.1
0.82 1:
1.00
2.24
SP2-gumbel
30
54.01
+ 3.73
23.38
5 7.35
1 :: 2.92
1.76 2 2.79
1.69 1: 2,90
2.13 1
2.28 2 1.74
SE2 gumbel
50
97.21
1. 5.45
43.96 1: 37,99
10.24 1: 10.44
2.81 1: 3.30
2.43
3,03
2.64 1: 3,81
3.41 1: 4.63
3.77 1: 4.56
SF2-gumbel
100
217.
+: 8.29
112,49 66.23
22.39 19.17
4.89
1: 5,94
4.21 5.17
4.78 1:
5.59 5:42
6.63
1: 6.00
SFi-gumbel
16
+
1
4.05
2.331
1.75
1.23
gum-eumbel
30
65.91
1: 7.36
33.87
35.72 65.40
13:44 11.65
1.97 1 4.07
3.72
1: 4.11
4.13 1: 4.29
5.99 019-1:60
5.56 1: 6.92
SF3-gumbel
50
120.
+ 11.19
55.66
27.25 1: 20.76
9.43 1 7.56
7.29
2 7.07
5.79 1: 5.76
7.59
6:641
8.17 1: 8.16
SF3 gumbel
100
272.00
16.44
149.45 100.23
60.12 29.08
13.15
12.62
12.83
12,46 4: 12.84
12.45 13.35-19..0
13.87 4: 10,60
ST4 gumbel
3.177
7.27
1024
1.15 1 294
4,55 3.45
SF4-gumbel
30
73.54
9.23
45.15
41.75
29:45 1: 19.08
12.35 = 7.64
8.37
6.70
8.14 7.53
11
9.44
11.84
9.91
SF4 gumbel
50
137.23
4: 13.67
81.88 4: 67.58
55.57 4 28.67
18.03
15.20
12.07
4: 9.36
10,92 1: 904
13.44 +: 10.94
15.06
4: 11.37
SF4 gumbel
100
315.13
1: 23.64
185.00 111.88
114.10 1: 52.86
29,94
22.69
27,00
28.51
22.36 1: 21.7
21.24
1: 23.58
25.37
1: 25.67
22

---

Published as a conference paper at ICLR 2025
Table 9: DAG learning performance (measured in true positive rate, the higher the better, best results in bold) of different algorithms on small scale (10-100 nodes) ER|,,3,4) graphs with different noise distributions. Our algorithm performs better than previous approaches.
Graphs
Nodes
MMPC
GES
NOTEARS
DAGMA
Order-1
Order-2
Order-3
Order-1
ER2-gaues
10
0.59 0.14
0.70
+: 0.23
0.8R
0.09
0.91
= 0.07
0.97
0.04
0.97
:: 0.04
0.97
0.04
0.96
: 0.05
F82-gams
30
0.53
0.11
0.80
+ 0.14
0.93
0.::0.0 0.05
0.95
= 0.04
0.98
0.02
0.99
= 0.02
0.99
= 0M02
0.98
= 0.00
BR2-gams
50
0.55
+ 0.00
0.83
+ 0.30
0.93
+ 0.05
0.97
= 0:02
0.99
0:01
0.99
= 0:01
0.99
0:01
0.99
= 0.01
ER2-gaus
100
0.07
0.87
0.06
0.04
0.02
0.99
0.01
0.92
0.01
10
0.54
+ 0.07
0.8
0.75
0.80
= 0.10
0.91
0:08
0.1
= 0:08
0.91
0.07
091
- 0.07
BR3-gan gauss
30
0:29
0.07
0.63
0.13
0.88
1: 0.06
0,92
1 0.04
0.97
0.02
0,97
:: 0.03
0.97
2 0.03
0.97 1 0.03
FR3-pui
50
0.50 0.06
0.69
2 0.10
0.89
0.06
0.93
= 0.03
0
0.02
0.98
0.02
0.98
0.02
0.99
0.02
FR1-
100
+
0.78
0.07
0.95
0.98
0.01
0.9M
0,01
ERt-games
10
0.28
0.05
0:43
0.14
0.68
0.09
0.73
- 0.07
0.84
0.07
0,85
0.07
0.84
0.07
0,4
: 0.07
10
0.18
+ 0.04
0.54
0.10
0.80
= 0.08
0.87
= 0.04
0.94
0.03
0,95
= MM3
0.94
- 0.03
0.4
= 0.03
50
0.18 0.03
0.57
+ 0.00
0.8
0.05
0.91
0:03
0.96
0.02
0.96
0.02
0.96
0.02
0.96 0.02
0.68
0.07
0.93
0.97
0.97
0,01
FR2-cxp
10
0.58
14
0.67
+ 0.34
0.7
0.09
0.1
= 0.07
0.97
0.04
097
- 0.04
0.97
0:04
0,96
= 0.05
BR2 exp
30
0.54
1: 0.11
0,80
1: 0.34
0.93
0.05
0,95
: 0.04
0.98
0.02
0,98
= 0.02
0.99
0.02
0.98
:: 0.02
BR2
50
0.56
0,10
0.82
1 0.09
0.93
0.05
0.97
0.02
0.99
0.99
0.01
0.99
0.01
0,99
0.01
FR2-cxp
100
0.97
0.02
0.99
0.01
0.0
= 0.01
ER3-oxp exp
10
0.34
0.08
0.48
1: 0.16
0.75
0.11
0,81
= 0.09
0.90
0.07
0,9
:: 0.08
0.92
0.07
0,00
:: 0.07
BR3-<ap
10
0.29
1 0.07
0.64
+ 0.13
0.88
+ 0.06
0,92
= 0.04
0.97
0.02
0,96
= 0.03
0.97
- 0.03
0,97
0.03
FR FR1-axp
50
0.51
0.06
0.69
+ 0.10
0.88
0.06
0.93
0:03
0.97
0.02
0.9M
0.02
0.97
0.02
0.97
0.01
BR3-oxp exp
0.78
0.07
0.95
0.02
0.96
0.01
0,96
0,01
ERI-exp
10
0.28
0.05
0.44
1 0.13
0.68
0.09
0.72
= 0.07
0.84
0.08
0.84
= 0.07
0.84
0.08
084
: 0.07
BRI exp
30
0.17
0.04
0.54
+ 0.11
0.81
0.07
0.87
= 0.05
0.94
0,03
0.94
= 0:03
0.94
0.03
Q.4
= 0.04
BRI-00
50
0.18
0.03
0,58
1: 0.09
0.83
0.05
0,90
0.03
0.02
0,96
0,02
0.02
0,96
0.02
0.68
0.07
0.93
0.03
0.97
0.02
0.97
= 0,0
BR2 gumbel
10
0.57
0.15
0.64
1: 0.24
0.93
0.0S
0,96
=: 0.05
0.98
0.04
0.8
=: 0.04
0.98
0.04
0,88
:: 0.04
BR2 gumbel
10
0,54
1 0.11
0.78
1: 0.13
0.95
0.05
0.96
0.02
0.99
0.02
0,99
=: 0.02
0.99
: 0.02
0.99
: 0.02
BR2-g plb
50
0.56
0.10
0.83
2 0.09
0.95
0.04
0,98
0.02
0.99
0.99
0.01
0.99
0.01
0.99
0.01
BR2 gumbel
0.87
0.05
0.99
0.02
0.01
0.99
0.01
0.99
0,99
0.01
ER3- gumbel
10
0,35
0.08
0.49
1: 0.16
0.84
0.10
0.89
:: 0.07
0.95
0.06
0.94
:: 0.06
0.95
0.06
0,95
: 0.06
BK3 gumbel
20
0.28
0.077
0.63
+ 0.13
0.92
0.06
0.96
0.02
0.98
0.02
0.98
+ 0.02
0.97
0:02
0197
+ 0.02
BR3- gumbel
50
0.30
0.69
1: 0.10
0,96
0.02
0,98
0,01
0.6
0.02
BR3-g gumbel
0.78
0.08
0.97
: 0.02
0.9K
0.99
0.01
0.9
1 0,01
B84 gumbel
10
0.28
0.05
0.64
+ 0.13
0.71
0.09
0.82
+ 0.07
0.91
0.07
0.90
+: 0.07
0.89
0:07
0,9
+ 6.08
BRI gumbel
30
0.17
0.04
0.54
1: 0.11
0,87
0.07
0.93
:: 0.03
0.97
0.02
0,97
:: 0.02
0.96
0.03
0.96
:: 0.03
BR pumbel
50
0,18 0.03
0.57
1: 0.09
0.87
0.05
0,95
0.03
0.98
0.02
0,98
0,01
0.97
0.02
0,97
1 0.02
BR4 gumbel
100
0.69
0.07
0.96
0.02
0.98
0.02
0.98
0.01
SF2-2mas
10
0.78
0.11
0:91 + 0.19
0.91 = 0:47
0.97
0:84
0.9
+ 0.05
St2-gamas
B&2 0.09
0:95 + 0.05
0.95
0.95 = 0.04
0.98
0:02
0.99
= 0.02
0.99
0.02
0.8
+ 0.02
SF2-gmss
50
0.57 + 0.08
0.92 : 0.11
0.93 + 0.05
0.97 = 0.02
0.99
0.01
0.99
= 0.01
0.99
- 0.01
0,99
: 6.11
SF2-2mas
100
0.54 0.06
0.92 + 0.08
0.94
0.04
0.97
0.02
0.99
0:01
0.99
0.01
0.99
0:01
009
- 0.01
SF3 SF3-pnas
0,86
0.19
0,80
= 0.10
091
= 0.08
0.91
0407
160
= 0.7
SF1-gres
0.53 0.10
0.92
+ 0.00
0
0.06
0.92
= 0.04
0.97
0:02
087
= 0.03
0.97
= 0:03
SF3-gmes
50
0:47 + 0.00
0.90
+ 0.10
0.89
+ 0.06
0.93
= 0.83
0.9%
0:02
0,88
= 0.02
0.9%
0:02
= 0.02
SF3
100
0.43 0.05
0.90
1: 0.09
0.89
0,95
0.02
0.98
0.9
0.01
0.98
6UW 0.01
SF4-2ms
0.15
0.07
0,85
=
0x4
- 0.07
SF4-gmas
0.87
+ 0.12
0.00
0.08
088
= 0.04
0.94
0:03
0,9
0.4
0:03
094
= 6.03
SF4-pmss
50
0.42 : 0.07
0.88
1: 0.10
0.83
0.05
0.91
0.05
0.96
0.02
0,96
0002
0.96
- 0.02
096
= 0.02
SF4-gmas
100
0.59 0.05
0.90
: 0.00
0.85
0.06
0.93
0.03
0.97
0.02
007
8:02
0.97
6:02
0.97
0.01
SF2-esp
0.7
0.97
0%
= 0.05
SF2-exp
30
0.62
0.09
0.93
: 0.11
0.93
0.05
0,9
= 0.04
0.98
0.02
088
- 0.02
0.99
0402
0,8
= 0.02
SF2-esp
50
0.57
+ 0.08
0.92
: 0.11
0.93
0.05
0.97
0.02
0.99
0:01
0.99
= 0.01
0.99
= 0:01
0.99
= 0.01
SF2-exp
100
0.54 0.06
0.91
1: 0.08
0.94
0.97
0.02
0.99
0.99
0.01
0.9
0.01
0.1
0.00
- 0.07
SF1-c1p
0.53
10
0.92
+ 0.10
038
0.96
0.9
= 0.04
0.97
0:02
0.9
= 0:03
0.7
0:03
6.7
= 0.03
SP3-exp
50
0.47
1 0.07
0.91
1: 0.09
0.88
0.06
0.93
= 0.03
0.97
0.02
0,98
0
0.97
1 0.02
059
0.01
SF3-exp
100
0.43 0.05
0.89
2 0.10
0.99
0.04
0.95
0.02
0.98
0.01
0.98
6:01
0.98
0:01
0.01
SF4-csp
0.88
0.16
0.07
0.84
0,4
= 0.07
SF4-cxp
30
0.49
0.08
0.88
1: 0.12
0.81
0.07
0.87
: 0.05
0.94
0403
0.4
: 0.03
0.94
0.03
094
: 0.04
SF4-exp
50
0.42
+ 0.07
0.89
1 0.10
0.83
0.05
0
0.03
0.96
0:02
0.96
= 0/22
0.96
- 0.02
096
= 0.02
SF-
100
0.39
0.05
0,90
1: 0.10
0.85
0.93
0,03
0.97
0.97
0.02
0.97
0.02
0.7
1 0.01
3F2-gumbel
0.91
0.18
0,96
0.96
0.04
0.98
0.98
1 0.04
SF2-gumbel
0.62 0.00
0.94
2: 0.08
0.95
0.05
0.98
= 0.02
0.99
8:02
0.99
= 0.02
0.99
0:02
0.99
: 0.02
SF2 sumbe
50
0.57 0.08
0.92
1: 0.10
0.95
0.04
0.96
1 0.02
0.99
0.01
0.99
:: 0.01
0.99
1 0.01
0.99
: 0.01
SF2-zumbel
100
0.53 0.06
0.92
1: 0.08
0.95
0.03
0.99
0.02
0.99
0.01
0.99
6.01
0.99
0.01
0.99
1 0.01
SF1-gumhel
0.17
0.94
0.95
= 0,06
SF3 sumbe
30
0.54
0.10
0.92
1: 0.10
0.92
0.06
0.96
: 0.02
0.98
0.02
0.98
:: 0.02
0.97
0.02
0,7
1 0.02
SF3-gumbel
50
0.48 1 0.07
0.91
1 0.09
0.92
2 0.05
0.96
- 0.02
0.98
0.01
0
= 0.01
0.98
- 0.01
0.8
: 0.02
SF3 gumbel
100
0,90
0.99
0,97
0.02
0.99
0.01
099
1 0.01
SF4 gumbel
0.13
0.87
0.16
0.82
1 0.07
0.90
0.07
059
1 0.08
SF4-gumbe)
30
0.48
0.08
0.88
1: 0.13
0.87
0.07
0.3
- 0.03
0.97
0.02
0.97
: 0.02
0.96
0.03
096
1 0.03
SP4 gumbel
50
0.42
4: 0.07
0.89
4: 0.11
0.57
0.05
0.95
4 0.03
0.98
0.02
0.98
0.01
0.97
0.02
097
: 0.02
SF4 gumbel
100
0.39
0.05
0.89
1: 0.10
0.88
0.06
0.96
2 0.02
0.98
0.02
0.96
0.02
0.98
0.01
0.98
1 0.01
23

---

Published as a conference paper at ICLR 2025
Table 10: DAG learning performance (measured in true positive rate, the higher the better, best results in bold) of different algorithms on small scale (10-100 nodes) ER graphs with different noise distributions. Our algorithm performs better than previous approaches.
Graphs
Nodes
MMPC
GES
NOTEARS
DAGMA
Order-1
Order-2
Order-3
Order-4
ER2-gaues
10
0.55
0.05
0:46
+: 0.23
0,05
0.07 0.02 = 0.05
0.02
0.04
0.02
:: 0.05
0.02
0.06
0.03
: 0.08
F82-gams
30
0.56 0.03
0.45
0.21
0.04
0.06
0.02
= 0.04
0.02
0.04
0.02
= 0.03
0.02
= 0M03
0.02
= 0.044
BR2-gams
50
0.56
+ 0.03
0.4
+ 0.16
0,05
0.00
0.1
= 0:02
0.02
0.02
0.01
= 0.02
0.02
0.02
0.02
= 0.03
ER2-gaus
100
0.00
0.44
0.11
0.05
0.02
0.02
0.02
0.02
10
0.57
+ 0.05
0:60
+ 0.14
0.10
0.08
0.06
= 0.96
0.04
8:05
0.4
= 0:05
0M05
004
- 0.05
BR3- gauss
30
0.60
1: 0.04
0.70
1: 0.13
0.08
1: 0.08
0.03
= 0.04
0.03
0.04
0.03
:: 0.04
0.03
1 0.04
0.03
: 0.04
FR3-pun
50
0.59 0.03
0.71
+ 0.10
0.10
0.07
0.03
0.03
0.03
0.03
0.03
= 0.03
0.03
0.03
0.03
: 0.03
100
+
0.12
0.03
0.02
0.02
0.03
ERt-games
10
0.53
0.04
0.57
1: 0.13
0.09
0.06
0.06
= 0.05
0.05
0.04
0.05
2: 0.04
0.04
0.04
0,05
:: 0.04
10
0.61
0.04
0.76
4: 0.07
0.14
- 0.10
0.06
0:05
0.08
0.04
0.04
= 0.44
0.08
- 0.05
0.06
= 0.06
50
D.62
0.03
0.80
+ 0.05
0.15
0.07
0.04
0.34
0.04
0:03
0.04
0,03
0.04
0:04
0,05
0.04
guuss
0.77
0.06
0.04
0.04
0.04
0.04
FR2-cxp
10
0.55
0.05
0.4
+ 0.34
0.05
0.07
0.03
= 0.05
0.02
0.04
0.02
= 0.05
0.03
0.06
0,04
= 0.08
BR2 exp
30
0.56
1: 0.03
0.47
1: 0.21
0,05
0.07
0.02
0.04
0.02
0.04
0.02
= 0.04
0.01
0.03
0.02
: 0.04
BR2
50
0.56 0.03
0.49
1 0.15
0,05
0.06
0.01
0.02
0.01
0.01
0.02
0.02
0.02
0,02
: 0.03
FR2-cxp
100
0.12
0,00
+ 0.02
ER3-oxp exp
10
0.57
0.06
0.60
1: 0.14
0.11
0.08
0,06
0.06
0.04
0.05
0,04
:: 0.05
0.03
0.04
004
: 6.05
BR3-<ap
10
0.60
1 0.04
0.69
+ 0.13
0.08
+ 0.07
0.03
0.04
0.03
0.04
0.03
= 0.04
0.03
- 0.04
0.03
= 0.04
FR FR1-axp
50
0.59
0.03
0.70
+ 0.12
0.10
0.08
0.03
0.04
0.03
0.04
0.03
0,03
0.03
0.03
0,03
0.03
BR3-oxp exp
0.64
0.11
0.03
0.03
0.03
0,03
0.02
ERI-cxp
10
0.53
0.03
0.57
+ 0.12
0.09
0.06
0.07
= 0.05
0.05
0.04
0.05
- 0.04
0.04
0.04
0,0
: 0.04
BRI exp
30
0.61
0.04
0.76
+: 0.07
0.13
0.09
0.06
= 0.05
0.04
0.04
0.05
=: 0.04
0.06
0:05
005
+ 0.06
BRI-08
50
0.62
0.03
0.60
1: 0,05
0.15
0,05
0.04
0.04
0,03
0.04
0,5
0.04
0.77
0.06
0.04
0.04
000
0.04
BR2 gumbel
10
0.55
0.04
0.51
1: 0.23
0,03
0.06
0.02
0.05
0.02
0.06
0.03
+: 0.06
0.03
0.06
0.0
: 0.08
BR2 gumbel
10
0,56
1 0.03
0.49
1: 0.19
0,05
0.07
0.02
0.03
0.02
0.04
0.02
= 0.04
0.03
0.04
0.03
: 0.05
BR2-g plmbel
50
0.56
0.03
0.48
2: 0.16
0,06
0.06
0.02
0,03
0.02
0.03
0.02
0.03
0.03
0.04
0.03
0.04
BR2 gumbel
0.44
0.12
0,01
0,02
0.02
0.02
0.03
0.03
ER3- gumbel
10
0.57
0.05
0.60
1: 0.15
0.07
0.07
0.04
:: 0.05
0.03
0.05
0.04
:: 0.06
0.03
0.05
0.02
: 0.04
BK3 gumbel
20
0.60
0.04
0.70
+ 0.13
0.08
0.08
0.02
0,03
0.03
0.04
0.03
+ 0.04
0.05
0:05
0.4
+ 0.05
BR3- gumbel
50
0.59
0.70
1: 0.12
0.10
0.03
0,03
0.03
0,03
0.04
001
0.04
BR3-g gumbel
0.64
0.12
0,03
0,03
0.03
0,03
000
1 0,03
B84 gumbel
10
0.53
0.04
0.56
+ 0.13
0,08
0.06
0.05
0.04
0.03
0.04
0.44
+ 0.04
0.04
0:04
004
+ 0.04
BRI gumbel
30
0.61
1: 0.04
0.75
1: 0.08
0.13
2: 0.09
0.04
:: 0.04
0.04
0.04
0.04
:: 0.04
0.06
1 0.06
0.6
: 0.06
BR1-g pumbel
50
0.61 0.03
0.80
1: 0.05
0.15
0,08
0,05
0.05
0.04
0.04
0.04
0,03
0.05
0.04
0.06
1 0.04
BR4 gumbel
100
0.77
0.05
0.04
0.04
0.04
0.04
0.05
0.04
SF2-2mas
10
0.78 +
0:91 + 0.19
0.91 = 0.37
0.97
0:84
0.04
0.
= 0.05
St2-gamas
B&2 0.00
0:95 + 0.05
0.95
0.95 = 0.04
0.98
0:02
0.99
= 0.02
0.99
0.02
0.8
+ 0.02
SF2-gmos
50
0.57 + 0.08
0.92 : 0.11
0.93 + 0.05
0.97 = 0.02
0.99
0.01
0.99
= 0.01
0.99
- 0.01
0,99
: 6.11
SF2-2mas
100
0.54 0.06
0.92 + 0.08
0.94
0.04
0.97
0.02
0.99
0:01
0.99
0.01
0.99
0:01
009
- 0.01
SF3 SF3-pnas
0,86
0.19
0,80
= 0.10
091
= 0.08
0.91
0407
160
= 0.0
SF1-gres
0.53 0.10
0.92
+ 0.00
0
0.06
0.92
= 0.04
0.97
0:02
087
= 0.03
0.97
= 0:03
SF3-gmes
50
0:47 + 0.00
0.90
+ 0.10
0.89
+ 0.06
0.93
= 0.83
0.9%
0:02
0,88
= 0.02
0.9%
0:02
= 0.02
SF3
100
0.43 0.05
0.90
1: 0.09
0.89
0,95
0.02
0.98
099
0.01
0.98
OUW 0.01
SF4-2ms
0.15
0.73
0.07
0,85
=
0x4
- 0.07
SF4-gmas
0.87
+ 0.12
0.00
0.08
088
= 0.04
0.94
D:03
0,9
0.4
0:03
094
= 6.03
SF4-pmss
50
0.42 : 0.07
0.88
1: 0.10
0.83
0.05
0.91
0.05
0.96
0.02
0,96
0002
0.96
- 0.02
096
= 0.02
SF4-gmas
100
0.59 0.05
0.90
: 0.00
0.85
0.06
0.93
0.03
0.97
0.02
007
8:02
0.97
6:02
0.97
0.01
SF2-esp
0.7
0.97
0%
= 0.05
SF2-exp
30
0.62
0.09
0.93
: 0.11
0.93
0.05
0,9
= 0.04
0.98
0.02
088
- 0.02
0.99
0402
0,8
= 0.02
SF2-esp
50
0.57
+ 0.08
0.92
- 0.11
0.93
0.05
0.97
0.02
0.99
0:01
0.99
= 0.01
0.99
= 0:01
0.99
= 0.01
SF2-exp
100
0.54 0.06
0.91
1: 0.08
0.94
0.97
0.02
0.99
0.99
0.01
0.9
0.01
0.1
0.00
- 0.07
SF1-c1p
0.53
10
0.92
+ 0.10
038
0.96
0.9
= 0.04
0.97
0:02
0.9
= 0:03
0.7
0:03
6.7
= 0.03
SP3-exp
50
0.47
1 0.07
0.91
1: 0.09
0.88
0.06
0.93
= 0.03
0.97
0.02
0,98
0
0.97
1 0.02
059
0.01
SF3-exp
100
0.43 0.05
0.89
2 0.10
0.99
0.04
0.95
0.02
0.98
0.01
0.98
6:01
0.98
0:01
0.01
SF4-csp
0.88
0.16
0.07
0.84
0,4
= 0.07
SF4-cxp
30
0.49
0.08
0.88
1: 0.12
0.81
0.07
0.87
: 0.05
0.94
0403
0.4
: 0.03
0.94
0.03
094
: 0.04
SF4-exp
50
0.42
+ 0.07
0.89
1 0.10
0.83
0.05
0
0.03
0.96
0:02
0.96
= 0/22
0.96
- 0.02
096
= 0.02
SF-
100
0.39
0.05
0,90
1: 0.10
0.85
0.93
0,03
0.97
0.97
0.02
0.97
0.02
0.7
1 0.01
3F2-2umbel
0.91
0.18
0,96
0.96
0.04
0.98
0.98
1 0.04
SF2-gumbel
0.62 0.00
0.94
2 0.08
0.95
0.05
0.98
= 0.02
0.99
8:02
0.99
= 0.02
0.99
0:02
0.99
: 0.02
SF2 sumbe
50
0.57 0.08
0.92
1: 0.10
0.95
0.04
0.96
1 0.02
0.99
0.01
0.99
:: 0.01
0.99
1 0.01
0.99
: 0.01
SF2-zumbel
100
0.53 0.06
0.92
1: 0.08
0.95
0.03
0.99
0.02
0.99
0.01
0.99
6.01
0.99
0.01
0.99
1 0.01
SF1-gumhel
0.17
0.94
0.95
= 0,06
SF3 sumbe
30
0.54
0.10
0.92
1: 0.10
0.92
0.06
0.96
: 0.02
0.98
0.02
0.98
:: 0.02
0.97
0.02
0,7
1 0.02
SF3-gumbel
50
0.48 1 0.07
0.91
1 0.09
0.92
2 0.05
0.96
- 0.02
0.98
0.01
0
= 0.01
0.98
- 0.01
0.8
: 0.02
SF3 gumbel
100
0,90
0.99
0,97
0.02
0.99
0.01
099
1 0.01
SF4 gumbel
0.13
0.87
0.16
0.82
1 0.07
0.90
0.07
059
1 0.08
SF4-gumbe)
30
0.48
0.08
0.88
1: 0.13
0.87
0.07
0.3
- 0.03
0.97
0.02
0.97
: 0.02
0.96
0.03
096
1 0.03
SP4 gumbel
50
0.42
4: 0.07
0.89
4: 0.11
0.57
0.05
0.95
4 0.03
0.98
0.02
0.98
0.01
0.97
0.02
097
: 0.02
SF4 gumbel
100
0.39
0.05
0.89
1: 0.10
0.88
0.06
0.96
2 0.02
0.98
0.02
0.96
0.02
0.98
0.01
0.98
1 0.01
24

---

Published as a conference paper at ICLR 2025
PC
GES
DAGMA
Exponential
Order
Order
Order 3
Order 4
SHD 563.9 23.84 4490.2 d: 62.52
588.8 18.33
488.61 24.29
429.6 24.73
410.6 15.25
401 0 16.64
389.4 +: 16.70
Exp MLE
Order MLE
Order MLE
Order MLE
Order MLE
SHD
518.00 23.02
453.70 + 42.12
447.30 + 51.85
409.50 + 31.02
433.00 68.98
PC Exp
PC Order-1
PC Order-2
PC Order-3
PC Order-4
SHD
275.40 16.01
274.40 # 15.44 273.10 15.72
271.80 + 14.75
276.00 + 14.66
PC Exp MLE
PC Order-I MLE PC Order-2 MLE PC Order-3 MLE PC Order-4 MLE
272.20
14.04
273.00
17.79 270.20
12.58
PC
GFS
DAGMA
Exponential
Order
Order
Order 3
Order 4
SHDD 321.30 10.30-1
27.77 4626.20 i: 69.05
674.00
466.50
458.40
447.20
Exp MLE
Order
MLE
Order
MLE
Order
MLE
Order
MLE
SHDD
$74.50 1. 42.84 490.80 1. 66.99
486.80 1: 76.47
444.30
42.22
479.50 di 100.71
PC Exp
PC Order-1
PC Order-2
PC Order-3
PC Order-4
SHDC
236.20 15.16 236.20
15.16
236.20
15.16
236.20
15.16
236.20 15.16
Exp
Order-3
SHDD
231.30 13.64 257.50 26.14 236.10 1: 16.23 236.80 23.77
231.60 13.17
Table 11: DAG learning performance (measured in structural hamming distance, the lower the better, best results in bold) different algorithms on 1000-node ERI graphs with Gaussian noise with observation data normalized. Our algorithms performs better than the previous approaches, and higher order DAG constraints suffers less to gradient vanishing, it tends to have better performance. We compare differential DAG learning approaches with conditional independent test based PC (Spirtes and Glymour 199T> algorithm and score based GES (Chickering (002) algorithm. The result is reported in the format of average standard derivation gathered from 10 different simulations. The results are reported as averages + standard deviations, calculated from 10 independent simulations. In addition to the MSE score function, we also applied the MLE score function described in Ng et al. (2020). Furthermore, rather than only considering edges between variables with correlation coefficients greater than 0. we also evaluated cases where edges are restricted to those in the PC-estimated CPDAG (algorithms whose names begin with 'PC").
F IMPLEMENTATION FOR ALGORITHM I
def gradd W, eps-le-20)
M. /
Iw self. Id - M_ # self.Id is identity matrix
icnt
Inv self. Id + M
while ient < 2 * self.d:
M_ M_ 8 M.
Inv Inv + Inv M_
icnt x- 2
if self.np.mn. (wwlf.np......... (M)) eps:
break
if anlf.np.anyyselff.. isnan (Inv)) break
if self.np... any (self. (Inv))
return sels.ap.
www.s_erossseeeee (Inv)
if self.np.anc (self.np.isna)) (Inv))
return self. zeros_ (Inv)
return Inv / 8
def
W, a):
25

---

Published as a conference paper at ICLR 2025
191
AMASE 018112
AAASEA ONRY
Frobenius norm of gradients vs. gradient descent steps for the first two iterations in Algorithm 2. Bottom: Frobenius norm of gradients vs. gradient descent steps for the last three iterations in Algorithm 2. In most cases, our higher-order DAG constraints exhibit larger gradient norms compared to DAGMA, enabling our algorithm to often converge to better solutions than DAGMA.
Figure 1: Gradient Norm v.s. Optimization Iterations. Top:
M self. self._h_grad (W, s)
if self. self.np.any(es (self. isnan (M) or self ap.linalg... norm self.Id, ord='fro')
M 0 (s self. W 1e-6:
if isinstance (W, supy.ndarres
-r s, v supy.linalq.ae. svd (W W) # cupy does not have eig 1ib, thus use spectral norm as an estimation
cs cwwy.max (s) + 0.1 self.h_order
else:
CS np.max (np. abs (np. www.iaaggeigvalle..
M * W)) 0.1 self.h_ordes
else:
return M, cs
G EMPIRICAL RESULTS ON GRADIENTS VANISHING
In this section, we present empirical results on gradient issues in DAG learning. According to Proposition 5, the gradients of higher-order DAG constraints should have larger norms than those of lower-order constraints, given the same candidate adjacency matrix. Additionally, the behavior of our Order-1 DAG constraints is expected to align closely with that of DAGMA.
In the path-following algorithm described in Algorithm 2, five iterations are used to tune the scale of the score function. During each iteration, tens of thousands of gradient descent steps are performed. We recorded the gradient norms for various DAG constraints over the five iterations using 1000-node ER3 DAG learning problem with Gaussian noise. The results are shown in Figure [ In most cases. our higher-order DAG constraints exhibit larger gradient norms compared to DAGMA, enabling our algorithm to often converge to better solutions than DAGMA.
26


Reviews:
[{"review_id": "dYbPHC0LKP", "reviewer": "Reviewer_xBhf", "summary": "The submission presents a mathematical derivation of a family of analytical DAG constraints that can be used in the recent proposals of smooth optimization for DAG search, \nTheoretical results are proposed to show that the derivation is valid and that the proposed family has favorable properties to be used in the optimization algorithms. \nAdditionally, a large set of experiments is performed.", "questions": "### suggestions\n\nFor general suggestions see Weaknesses, additionally:\n\n1. In the first paragraph of the Introduction: \"the recovered DAGs can constitute a causal graphical mode\"  --> \"the recovered DAGs could be interpreted causally\" or something similar, as the DAG by itself can't be a causal model since it lacks at least a compatible model for the probabilities or more generally a functional causal model.\n2. End of second paragraph in Introduction: \"w.r.t number of nodes(Chickering, 1996\" a space is missing after \"nodes\".\n3. In \"as demonstrated in Wei et al. (2020). This study reveals that\"  The  \"This\" in the beginning of the second sentence refer to Wei et al. (2020) but it could be mistakenly taken as \"This study == this paper\" I suggest to change the wording to clearer reflect that the following is work in Wei et al. (2020).\n4. Section 2 Preliminaries: Similar to suggestion 1. above,  I think $\\mathcal{G}$ should only denote a DAG (that is a graph). And not a DAG model (which is not even clear what a DAG model is). Then given a DAG the authors could define what a compatible statistical model or directly a compatible SEM is (that is when the probability factorizes according to the DAG or for the linear SEM with indep. noise that the matrix B encodes a weighted adjeceny matrix of the graph G). \n5. It would be nice to have an extended intuitive explanation for the following: \"Due to the fact that the adjacency matrix of a DAG must form a nilpotent matrix, we do not need a function with infinite convergence radius\"\n6. before section 4 starts, \"mainly due to the non-convexity of the overall objective functionNg et al. (2023).\" a space is missing after \"function\" and before \"Ng et al.\". \n7. In caption of Table 1: \"t state-of-the-arts DAGMA(Bello et al., 2022),\" a space is missing after DAGMA. \n8. Table 3, is the value for GES wrong? it seems to be one order of magnitude bigger ?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is generally well written (except for some minor details, see suggestion section in Questions) and well structured. \n- The experiments section is quite extensive (even if the number of replications could be bigger). \n-  The theoretical part seems sound and it is well written (I have not checked the full proofs on the supplementary). \n- The subject is interesting even if theoretically it is restricted to equal noise variance and linear SEM.", "weaknesses": "- In the simulation experiments, it would be nice to have a larger number of replicates instead of 10.  I understand the time complexity but \nmaybe lower dimensional systems could be explored with a higher number of repetitions to have a set of more robust results that could complement the proposed ones. \n- In the experiments section the CausalDiscoveryToolbox is used to run some of the state-of-the-art methods. I would discourage that. the CausalDiscoveryToolbox would be a nice tool, but unfortunately, it is no longer maintained, While probably there are no issues in running the wrapper of other algorithms provided there I would suggest the author of this submission use other implementations. For example, PC and GES in the CausalDiscoveryToolbox are wrappers of the pc-alg R implementations, I would suggest to use directly those for reproducibility and additional safety. \n- Is the code for the experiments and the implementation available somewhere? I did not find it in the submission", "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": "l6jnvi1Ttc", "reviewer": "Reviewer_abAq", "summary": "the paper proposes a framework on continuous DAG constraints, based on analytic functions. It shows many existing dag constraints can be unified under such a framework. Specific choices on the analytic functions can improve computation efficiency, and in experiments even accuracy due to vanish gradient problems.", "questions": "see above\n\n---------------------- after rebuttal\nthank the authors for the new updated results. It would be good to have a formal statement around error bounds, which is not yet in the draft. I have raised my score accordingly.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The paper investigate in depth some properties of analytic function based DAG constraints, which shows some insights \n\n- Empirically extensive experiments are performed. The results show the method can outperform existing baselines.", "weaknesses": "- Some convergence statements, such as convergence rate and eror bound, of using these analytic functions would be good to further improve the understanding of these constraints. \n\n- Experiments contain some mixed messages, and it seems performance depends on data property and choices of orders. These make the proposed constraints harder to use in practice.\n\nSpecific Comments:\n- L136: f(B) in the form of eq 3 or eq 4?\n- L152: how efficient of these algorithms?\n- L158: beside exponential, polynomial also holds the property. Any other function classes with such property? Could they be used as DAG constraints, any pro and cons?\n- L224: how much can this acutally save in practice?\n- What is the complexity of Algoritahm 1, compared to standard computations?\n- L305~307: missing sentences here?\n- Given the convergence property of analystical functions, one natural question to ask is: what are the convergence rates and error bounds of  the formulation using these function based dag constraints? how does different convergence radius affect these bounds?\n- Experiments shows different orders of the propsed constraint, but it is unclear which order should be used in practice. \n- Computation time comparison on constraints here would also be good to justify author's claim on efficiency. \n- What are the differences between table 5 and 6? Titles are the same but the results are very different. Table 6 and 7 also show DAGMA are better, instead of the claim in table titles? Is the message that DAGMA should be use for large scale graphs? I assume this is true because longer paths exists and hence needs more orders?", "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": "d3zgU0EqHq", "reviewer": "Reviewer_5eCV", "summary": "The paper tackles the gradient vanishing problem of the DAG constraints in differentiable DAG learning methods. The main contributions are:\n1. Showing that any analytic function can be used as a DAG constraint and unifying existing constraints into this framework. In addition, the closeness property of the function class allows for easy generation of constraints with large gradient magnitude.\n2. An algorithm for choosing the spectral radius to amplify the gradient during training.\n3. Demonstrating the outperformance of the proposed method through extensive experiments on linear simulated data and non-linear real data.", "questions": "1. For the experiments on linear SEM with unknown variance, which score function is used? MSE or the log-likelihood? Is the SHD on the Markov equivalence class?\n2. Does Table 7 describe the runtime? The caption seems to be incorrect.\n3. It would be nice if the author could compare the choice of $s$ (and/or gradient magnitude) for the proposed method and DAGMA on an example to visualize the effectiveness of Algorithm 1 and the DAG constraint on dealing with gradient vanishing.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The class of analytic DAG constraints unifies the existing DAG constraints. In addition, the author provides a simple solution to amplify the gradient by multiplying the constraint function.\n2. The discussion on the convexity of the DAG constraints provides intuition on the behavior of the high-order constraint, suggesting a tradeoff between the convexity and the gradient magnitude may be considered.\n3. The experiments are extensive. The proposed method is evaluated on different DAGs ranging from low-dimensional to high-dimensional and the outperformance is well-established.", "weaknesses": "1. Some statements could be more rigorous. For example, \n   * Lines 304-307, the discussion of the optimization objective (17). Why the statement is only restricted to Gaussian SEMs? According to Loh and Bühlmann, (2013), MSE is valid for the linear equal variance model and is proved to have good performance by Bello et al. (2022). For unequal noise variance, MSE can be inconsistent. The second sentence is also not well-structured. Why is this a problem with convexity rather than the score function?\n   * Lines 362-367, same as above, MSE is not a consistent object for unequal scale, and therefore, it may not be a matter of the convexity that leads to bad results.\n2. I think the evaluation of the method on the non-linear data could be even more sound if the author also considers simulations on synthetic data with equal noise variance as Bello et al. (2022).", "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": "6Utlu1Mb3F", "reviewer": "Reviewer_Khru", "summary": "The authors provide a theoretical analysis of analytic matrix functions serving as acyclicity constraints in the context of differentiable DAG learning.\nThey first show that analytic functions are valid DAG constraints, that the gradients of analytic functions are easily computable, and that analytic functions do not suffer from vanishing gradients (in contrast to other acyclicity constraints). Next, the authors show that the family of analytic functions is closed under differentiation, addition, and multiplication, thus introducing an entire family of new acyclicity constraints. Additionally, a theoretical analysis of the non-convexity of said constraints is provided, indicating that there is a trade-off between high/informative gradients and non-convexity that has to be considered when designing acyclicity constraints. Finally, an experimental evaluation shows that differentiable DAG learners equipped with acyclicity constraints based on analytic functions perform better than other constraint functions.", "questions": "- line 195-198: that's an interesting property but what is it actually good for? Because using the gradient/derivative of an analytic function essentially yields another analytic function with a shift by one in the polynomial degrees. As we only need a nilpotent matrix to ensure DAGness, this property may not necessarily yield a benefit? Could the authors please provide some more information on that?\n-  paragraph in line 356-367: Since most of your experiments are built upon this paragraph, could the authors please provide more intuition on why the known/unknown variance has an influence on MSE. What about other losses like LL? Is this related to the findings of [2]?\n\n\n# References\n[1] Reisach et al. Beware of the simulated dag! causal discovery benchmarks may be easy to game. NeurIPS 2021.\n\n[2] Seng et al. Learning Large DAGs is Harder Than You Think. ICLR 2024.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- the paper provides an interesting and important contribution to analyzing acyclicity constraints for differentiable DAG learners formally\n- the theoretical contributions (i.e., the propositions in 3 and 4) are mathematically sound, and the proofs are correct. Note that I rated soundness with \"fair\" due to major weaknesses in the experimental section. The theory section is strong.\n- the work provides a theory-backed \"out-of-the-shelf\" approach to design acyclicity constraints using differentiation, addition or multiplication of analytic functions. This simplifies the process of constructing new acyclicity constraints.\n- it is shown that the proposed family of acyclicity constraints does not suffer from vanishing gradients\n- the experimental section shows that analytic acyclicity constraints indeed can lead to better performances in DAG learning", "weaknesses": "- line 106-107: \"However, these methods are known to perform poorly when applied to normalized data since they rely on scale information across variables for complete DAG recovery.\" This statement is a bit confusing since the findings in [1] and [2] (which generalized [1]) show that it is not the acyclicity constraint that leads to bad performance but the loss itself (MSE in this case).\n- line 127-129: the claim that Eq. 4 converges if the spectral radius of $\\mathbf{A}$ is smaller than $r$ is correct. However, it would be good to provide some intuition (and maybe proof in the appendix) since most propositions are based on that claim/assumption.\n- line 154-156: while I agree that in terms of computational efficiency, it is desirable to have acyclicity constraints whose gradient is representable by itself, the evaluation of these constraint functions can still be costly if not carefully chosen. It would be good if the authors could elaborate a bit more on that and provide an intuition why analytic functions are computationally efficient beyond the gradient computation.\n- Eq. 10: Could the authors provide a more detailed derivation of the constraints? It's easier to see where the constraints in Eq. 10 come from if Prop. 1 and 2 are actually applied to the functions in Eq. 9.\n- Eq. 11: shouldn't it rather be $\\frac{\\partial f^s_{\\text{log}}(x)}{\\partial x} = -s (s -x)^{-1} = -s \\cdot f^s_{inv}$?\n- line 223: the authors claim that $(\\mathbf{I} - \\mathbf{B} / s)^{-1}$ can be cached. However, $\\mathbf{B}$ will change in each iteration and thus hast to be re-computed, right?\n- Prop. 5: Couldn't the fact that the gradient of the trace of $(\\mathbf{I} - \\mathbf{B} / s)^{-n}$ is smaller than that of $(\\mathbf{I} - \\mathbf{B} / s)^{-n-k}$ potentially lead to exploding, rather than vanishing gradients? A brief explanation on this would be appreciated.\n- Algorithms 1 and 2 are not properly discussed in the main text. Especially, it is not entirely comprehensible from the text what the parameters of the algorithms mean.\n- the text in Sec. 4 and 5 reads a bit \"rushed\" and is not polished.\n- I think the experiments miss the main message the authors try to convey: In Sec. 3 and 4 the authors provide a strong theoretical analysis of the acyclicity constraints. However, the experiments seem to be mostly based on the paragraph in line 356-367 where the authors briefly discuss the learning dynamics of differentiable DAG learners when the scale of the data is known/unknown which is not the main contribution of this paper. The authors only show that analytical DAG constraints yield better results on simulated data (except Tab 4 where results are not significant and standard deviation is missing). While this is an important result and demonstrates that the proposed family of acyclicity constraints positively affects learning on simulated data, it remains unclear from an empirical perspective if the initial goal of solving the vanishing gradient problem was achieved. To underline the paper's important theoretical contributions, the authors should revise their experiments and rather analyze the learning dynamics and gradient statistics over learning instead of focusing on evaluating the methods' accuracies with the new acyclicity constraints. Also, since the acyclicity constraint is not primarily affected by data scale, the focus of the experimental evaluation is questionable. \n\n# References\n[1] Reisach et al. Beware of the simulated dag! causal discovery benchmarks may be easy to game. NeurIPS 2021.\n\n[2] Seng et al. Learning Large DAGs is Harder Than You Think. ICLR 2024.", "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)