DAG Learning on the Permutahedron
Authors: Valentina Zantedeschi
Luca Franceschi
Jean Kaddour
Matt Kusner
Vlad Niculae
Venue: ICLR 2023 poster (2023)
Keywords: directed acyclic graph
causal discovery
differentiable sorting

Abstract: We propose a continuous optimization framework for discovering a latent directed acyclic graph (DAG) from observational data. Our approach optimizes over the polytope of permutation vectors, the so-called Permutahedron, to learn a topological ordering. Edges can be optimized jointly, or learned conditional on the ordering via a non-differentiable subroutine. Compared to existing continuous optimization approaches our formulation has a number of advantages including: 1. validity: optimizes over exact DAGs as opposed to other relaxations optimizing approximate DAGs; 2. modularity: accommodates any edge-optimization procedure, edge structural parameterization, and optimization loss; 3. end-to-end: either alternately iterates between node-ordering and edge-optimization, or optimizes them jointly; We demonstrate, on real-world data problems in protein-signaling and transcriptional network discovery, that our approach lies on the Pareto frontier of two key metrics, the SID and SHD.

DAG LEARNING ON THE PERMUTAHEDRON
Valentina Zantedeschi
Luca Franceschi
Jean Kaddour
ServiceNow Research vzantedeschi @emmall....mm com
Amazon Web Services"
University College London,
f cranuluc@amazon.
Centre for AI
Joan. kaddour 200ucl ac.
Matt J. Kusner
Vlad Niculae
University College London, Centre for AI
Informatics Institute, University of Amsterdam v. ni cu culas@uva.n)) n1
m. w..kkusnrrrull... ac. uk
ABSTRACT
We propose continuous optimization framework for discovering a latent directed acyclic graph (DAG) from observational data. Our approach optimizes over the polytope of permutation vectors, the so-called Permutahedron, to learn topo- logical ordering. Edges can optimized jointly, or learned conditional on the ordering via a non-differentiable subroutine. Compared to existing continuous optimization approaches our formulation has a number of advantages including: 1. validity: optimizes over exact DAGs as opposed to other relaxations optimizing approximate DAGs; 2. modularity: accommodates any adge-optimization proce- dure, edge structural parameterization, and optimization loss; 3. end-to-end: either alternately iterates between node-ordering and edge-optimization, or optimizes them jointly. We demonstrate, real-world data problems in protein-signaling and transcriptional network discovery, that our approach lies on the Pareto frontier of two key metries, the SID and SHD.
INTRODUCTION
In many domains, including cell biology (Sachs al., 2005), finance (Sanford Moosa, 2012), and genetics (Zhang et 2013), the data generating process is thought to be represented by underlying directed acylic graph (DAG). Many models rely on DAG assumptions, c.g., causal modeling uses DAGs model distribution shifts, ensure predictor fairness among subpopulations, or learn agents more sample-efficiently (Kaddour et al., 2022). A key question, with implications ranging from better modeling to causal discovery, is how to recover this unknown DAG from observed data alone. While there are methods for identifying the underlying DAG if given additional interventional data (Eberhardt, 2007; Hauser & Bühlmann, 2014; Shanmugam et al., 2015; Kocaoglu al., 2017; Brouillard et al., 2020; Addanki et al., 2020; Squires et al., 2020; Lippe et al., 2022), it is not always practical or ethical to obtain such data (e.g.. one aims discover links between dietary choices and deadly discases).
Learning DAGs from observational data alone fundamentally difficult for two reasons. (i) Estima- tion: it is possible for different graphs to produce similar observed data, either because the graphs are Markov equivalent (they represent the same set of data distributions) or because not enough samples have been observed to distinguish possible graphs. This riddles the search space with local minima; (ii) Computation: DAG discovery is a costly combinatorial optimization problem over an exponentially large solution space and subject to global acyclicity constraints.
To address issue (ii), recent work has proposed continuous relaxations of the DAG learning problem. These allow one to well-studied continuous optimization procedures to search the space of DAGS given a score function (e.g., the likelihood). While these methods are more efficient than combinatorial methods, the current approaches have one more of the following downsides: 1. Invalidity: existing methods based on penalizing the exponential of the adjacency matrix (Zheng
Work done prior joining Amazon.

---

et al., 2018; Yu et al., 2019; Zheng et al., 2020: Ng et al., 2020; Lachapelle al., 2020; He et al., 2021) are not guaranteed to return a valid DAG in practice (see Ng et al. (2022) for a theoretical analysis), but require post-processing to correct the graph to DAG. How the learning method and the post-processing method interact with each other is not currently well-understood; 2. Non- modularity: continuously relaxing the DAG learning problem is often done leverage gradient-based optimization (Zheng et al., 2018; Ng et al., 2020; Cundy et al., 2021; Charpentier et al., 2022). This requires all training operations to be differentiable, preventing the use of certain well-studied black-box estimators for learning edge functions; 3. Error propagation: methods that break the DAG learning problem into two stages risk propagating errors from one stage to the next (Teyssier & Koller, 2005; Bühlmann et al., 2014; Gao et al., 2020; Reisach et al., 2021; Rolland et al., 2022).
Following the framework of Friedman & Koller (2003), we propose a new differentiable DAG learning procedure based on a decomposition of the problem into: (i) learning topological ordering (i.e., total ordering of the variables) and (ii) selecting the best scoring DAG consistent with this ordering. Whereas previous differentiable order-based works (Cundy et al., 2021; Charpentier et al., 2022) implemented step (i) through the usage of permutation matrices', we take a more straightforward approach by directly working in the space of vector orderings. Overall, make the following contributions to score-based methods for DAG learning:
We propose a novel vector parametrization that associates a single scalar value to each node. This parametrization is (i) intuitive: the higher the score the lower the node is in the order; (ii) stable, as small perturbations in the parameter space result in small perturbations in the DAG space.
With such parameterization in place, we show how to learn DAG structures end-to-end from observational data, with any choice of edge estimator (we do not require differentiabblitt)... To do so, we leverage recent advances in discrete optimization (Niculae et al., 2018; Correia et al., 2020) and derive a novel top-k oracle over permutations, which could be of independent interest.
We show that DAGs learned with our proposed framework lie on the Pareto front two key metrics (the SHD and SID) two real-world tasks and perform favorably on several synthetic tasks.
These contributions allow us to develop a framework that addresses the issues of prior work. Specif- ically, our approach: 1. Models sparse distributions of DAG topological orderings, ensuring all considered graphs are DAGs (also during training); 2. Separates the learning of topological orderings from the learning of edge functions, but 3. Optimizes them end-to-end, either jointly or alternately iterating between learning ordering and edges.
RELATED WORK
The work on DAG learning can be largely categorized into four families of approaches: (a) combina- torial methods, (b) continuous relaxation, (c) two-stage, (d) differentiable, based.
Combinatorial methods. These methods are either constraint-based, relying on conditional indepen- dence tests for selecting the sets of parents (Spirtes et al., 2000), or score-based, evaluating how well possible candidates fit the data (Geiger& & Heckerman, 1994) (see Kitson et al. (2021) for a survey). Constraint- based methods, while elegant, require conditional independence testing, which is known to be hard statistical problem Shah & Peters (2020). For this reason, we focus our attention in this paper on score-based methods. Of these, exact combinatorial algorithms exist only for small number of nodes d (Singh & Moore, 2005; Xiang & Kim, 2013; Cussens, 201 1), because the space of DAGs grows superexponentially in d and finding the optimal solution NP-hard to solve (Chickering, 1995). Approximate methods (Scanagatta et al., 2015; Aragam & Zhou, 2015; Ramsey et al., 2017) rely on global or local search heuristics in order to scale to problems with thousands of nodes.
Continuous relaxation. To address the complexity of the combinatorial search, more recent methods have proposed exact characterizations of DAGs that allow one to tackle the problem by continuous optimization (Zheng et al., 2018; Yu et al., 2019; Zheng al., 2020; Ng et al., 2020: Lachapelle et al., 2020; He et al., 2021). 1). To do so, the constraint on acyclicity is expressed as a smooth function (Zheng et al., 2018; Yu et al., 2019) and then used as penalization term to allow efficient optimization.
'Critically, the usage of permutation matrices allows maintain fully differentiable path from loss to parameters (of the permutation matrices) via Sinkhom iterations or other (inexact) relaxation methods.
2

---

However, this procedure no longer guarantees the absence of cycles any stage of training, and solutions often require post-proesssing. Concurrently to this work. Bello et al. (2022) introduce a dog-determinant characterization and optimization procedure that is guaranteed to return DAG at convergence. In practice, this relies on thresholding for reducing false positíves in edge prediction.
Two-stage. The third prominent line of works learns DAGs in two-stages: (i) finding an ordering of the variables, and (ii) selecting the best scoring graph among (or marginalizing over) the structures that are consistent with the found ordering (Teyssier & Koller, 2005; Bühlmann al., 2014; Gao et al., 2020; Reisach et al., 2021; Rolland et al., 2022). As such, these approaches work over exact DAGs, instead of relaxations. Additionally, they work on the space of orderings which is smaller and more regular than the space of DAGs (Friedman & Koller, 2003; Teyssier & Koller, 2005), while guaranteeing acyclicity. The downside of these approaches that errors in the first stage can propagate to the second stage that is unaware of them.
Differentiable, obder- The final line of works uses the two-stage decomposition above, but addresses the issue of error propagation by formulating an end-to-end differentiable optimization approach (Friedman Cundy et al. (2021) optimizes node ordering using the polytope of permutation matrices (the Birkhoff polytope) via the Gumbel-Sinkhorn approximation (Mena et al., 2018). This method requires (d time and (d²) memory complexities. To lower the time complexity, Charpentier suggest to leverage another operator (SoftSort, Prillo & Fisenschlos, 2020) that drops a constraint on the permutation matrix (allowing row-stochastic matrices). Both methods introduce mismatch between forward (based on the hard permutation) and backward (based on soft permutations) passes. Further, they require all downstream operations be differentiable, including the edge estimator.
Koller, 2003; Cundy et al., 2021; Charpentier et al., 2022). In particular,
al. (2022)
3 SETUP
3.1 THE PROBLEM
Let X Rudd be a matrix of observed data inputs generated by an unknown Structural Equation Model (SEM) (Pearl, 2000). An SEM describes the functional relationships between d features as edges between nodes in DAG g E D D[d] (where D[d] the space of all DAGs with d nodes). Each feature Xj E R" is generated by some (unknown) function ʃJ of its (unknown) parents pa(j) as: X,- = (x(((())) ). We keep track of whether an edge exists in the graph G using an adjacency matrix A E {0, ] (i.e., A,y= =1 if and only if there is a (directed) edge > j). For example, a special case is when the structural equations are linear with Gaussian noise,
y, / (X, A,) = X (w,, o A₂) + E; 8 8~~ N(0,0)
where w, E Rd, A, is the jthh column of A, and is the noise variance. This is just add intuition; our framework is compatible with non-linear, non-Gaussian structural equations.
3.2 OBJECTIVE
Given X, our goal to recover the unknown DAG that generated the observations. To do so we must learn (a) the connectivity parameters of the graph, represented by the adjacency matrix A, and (b) the functional parameters = { that define the edge functions ( _ Score-based methods (Kitson et al., 2021) learn these parameters via a constrained non-linear mixed-integer optimization problem
Ex(x,, (X Az)) + XQ(D),
min AEDIA
(2)
where e: IR" x RD - IR is loss that describes how well each feature x, is predicted by f As in eq. (1), the adjacency matrix defines the parents of each feature x, as follows pa(j) = {i S
[d] \.j A.j = 1). Only these features will be used by f to predict xj. via X o A,. The constraint A E D[d] enforces that the connectivity parameters A describes a valid DAG. Finally, ?( is
regularization term encouraging sparseness. So long as this regularizer takes the same value for all DAGs within Markov equivalence class, the consistency results of Brouillard et al. (2020, Theorem 1) prove that the solution to Problem (2) is Markov equivalent to the true DAG, given standard assumptions.
3

---

equation parameters
(0) min E(
predicted
real data
structure parameters
graph
(unlnown)
0.6
0,1
50
min
0.6 &
Figure 1: DAGuerreotype: Our end-to-end approach to DAG learning works by (a) learning a sparse distribution over node orderings via structure parameters 0, and (b) learning sparse predictor w" to estimate the data, Any black-box predictor can be used to learn w"; differentiability is not necessary.
3.3 SPARSE RELAXATION METHODS
In developing our method, we will leverage recent works in sparse relaxation. In particular, we will make use of the SparseMAP (Niculae et al., 2018) and Top-k SparseMAX (Correia et al., 2020) operators, which we briefly describe below. At high level, the goal of both these approaches is to relax structured problems the form a* := arg a s a so that da* /Ös is well-defined
(where A = {a R a 0, Σ Og 1} is the D-dimensional simplex). This will allow to be learned by gradient- based methods. Note that both approaches require querying an oracle that finds the best scoring structures. We are unaware of such an oracle for DAG learning, i.e. for ID [dd being the vertices of A. However, we will show that by decomposing the DAG learning problem, we can find an oracle for the decomposed subproblem. We will derive this oracle Section 4 and prove its correctness.
Top-k sparsemax (Correia al., 2020). This approach works by (i) regularizing a, and (ii) constrain- ing the number of non-zero entries of a to be at most as follows k: arg max Jo lo s ||ax To solve this optimization problem, top-k sparsemax requires an oracle that returns the structures with the highest scores s ar.
MAP (Niculae et al., 2018). Assume has low-dimensional parametrization where B D and < D. SparseMAP relaxes ax = arg m Ba by regularizing the lower-dimensional 'marginal space' arg max o Ba |Ba l The relaxed problem can be solved using the active set algorithm (Nocedal & Wright, 1999), which iteratively queries an oracle
B r,
for finding the best scoring structure (r B TTa( (1) TB iteration +
4 DAG LEARNING VIA SPARSE RELAXATIONS
A key difficulty when learning DAG structures is that the characterization of the set of all valid DAGs: as soon some edges are added, other edges are prohibited. However, note the following key observation: any DAG can be decomposed as follows (i) Assign to each of the d nodes rank and reorder nodes according to this rank (this is called a lopological ordering); (ii) Only allow edges from lower nodes in the order to higher nodes, i.e., from node i node j if x; Xj. This approach lies at the core of our method for DAG learning, which we dub DAGuerreotype, shown in Figure 1. In this section the structural and the edge equations parameters leveraging the sparse relaxation methods introduced above and study the computational complexity of our method.
derive the framework, present global sensitivity result, show how to learn both
4.1 LEARNING ON THE PERMUTAHEDRON
Given d nodes, let Eg be the set of all permutations of node indices {1...., d). Given vector Idd, let v" := [ ..... (d) be the vector of reordered according the permutation a E Σ.. Similarly, for a matrix M R, let M be the matrix obtained by permuting the rows and columns of M by
4

---

Let Dc[d] be the set of complete DAGs (i.e., DAGs with all possible edges). Let R E {0, 1]axd be the binary strictly upper triangular matrix where the upper triangle is all equal to 1. Then Deld] can be fully enumerated given R and Σg, follows:
Beld] = {R* E Ed, R E {0. ddddd = R <t.
(3)
Therefore it is sufficient to learn a in step (i), and then learn which edges to drop in step (ii).
The vector parameterization. Imagine now that 0 R defines score for each node, and these scores induce an ordering o()): the smaller the score, the earlier the node should be the ordering. Formally, the following optimization problem finds such ordering:
E arg max 0 p" where p = [1, 2, ...
(4)
oEEa
Note that simple oracle solves this optimization problem: sort 0 into increasing order (as given by The Rearrangement Inequality (Hardy et al., 1952, Thms. 368-369)). We emphasize that we write 'E' in eq. (4) since the r.h.s can be set: in fact, this happens exactly when some components of 0 are equal. Beside being intuitive and efficient, the parameterization of [o given by 0 + R² allows us to upper bound the structural Hamming distance (SHD) between any two complete DAGs (I2 ) by the number of hhper-ppance of "equal coordinates" (H4, {x E Rd that are traversed by the segment connecting and l' (see Figure 4 in the Appendix for a schematic). More formally, state the following theorem.
and R²((")
Theorem 4.1 (Global sensitivity). For any O E IR and O' R
SHD (R R=(0")
ΣΣ
Sn, (0 + 1(0' - -0) dt
(5)
JEE(0,11
where (x) is the (generalized) Dirac delta that evaluates to infinity if x E A and otherwise.
In particular, Theorem 5 shows that we can expect that small changes in 0 (e.g. due to gradient- based iterative optimization) lead to small changes in the complete DAG space, offering result that is reminiscent to Lipschiz--smoothnee commentary (also compared to the parameterization based on permutation matrices) Appendix C.
for smooth optimization. We defer proof and further
Learning 0 with gradients. Notice that we cannot take gradients through eq. (4) because (a) (0) is even function (due to it possibly being set), and (b) even if we restrict the parameter space to have ties, the mapping is piece-wise constant and uninformative for gradient-based learning. To circumvent these issues, Blondel et al. (2020) propose to relax problem (4) by optimizing over the convex hull of all permutations of p. that is the the order-d Permutahedron P[d] == conv [p" a E E,), and adding convex regularizer. These alterations yield a class of differentiable mappings (soft permutations), indexed by E R+,
7 2 11
(0) arg max 0 μ PEd
which, in absence of ties, are exact for T + 0. This technique is, however, unsuitable for our case, as the (00)ss do not deseribe valid permutations except when taking values on vertices of P|d]. Instead, we show next how to obtain meaningful gradients whilst maintaining validity adapting the sparseMAP or the top-k sparsemax operators to our setting.
Leveraging sparse relaxation methods. Let D = d! be the total number permutations of d elements, and A D be the D-dimensional simplex. We can (non-uniquely) decompose μ = Eoer, an for some a E A. Plugging this into eq. (6) leads to
a SparseMPP (0) arg max 0 Eowalpo] 7 ||Eaaa 1p DEAD
2
We can recognize in (7) an instance of the SparseMAP operator introduced in Section 3.3. Among all possible decomposition, will favor sparse ones. This is achieved by employing the active set algorithm (Nocedal & Wright, 1999) which only requires access to an oracle solving eq. (4), ie.. sorting 0.
5

---

Alternatvely, because the only term in the regularization that matters for optimization is regularize it alone, and directly restrict the number non-zero entries of a to some kt > 2 as follows
jop-k sparsemax
(0)
arg max
θ E [po]
||a||
aFA Rel lkk
where we assume ties are resolved arbitrarily.
Algorithm 1: Top-A permutations.
This is formulation of top-k sparsemax introduced in Sec- tion 3.3 that we can efficiently employ to learn DAGs pro- vided that we have access to fast algorithm that returns set of k permutations with highest value of go(a) = 0 pa. Algorithm 1 describes such an oracle, which, to our knowl- edge, has never been derived before and may be of inde- pendent interest. The algorithm restricts the search for the best solutions to the set of permutations that are one adjacent transposition away from the best solutions found so far. We refer the reader Appendix A for notation and proof correctness of Algorithm 1.
Data: & [d!] 0 R Result: top-k permutations Tx(O) P(0) for En
arg Ea (o(o)); while (+) S do a En
arg urggaxxee PIE 90(o;
P(0) + ( JE (d [d-1]}: (0) (e)
end
Remark: Top-k sparsemax optimizes over the highest-scoring permutations, while sparseMAP draws any set of permutations the marginal solution can be decomposed into. As an example, for 0 = 0, sparseMAP returns two permutations, and its inverse o On the other hand, because at 0 all the permutations have the same probability, top-k sparsemax returns an arbitrary subset of k: permutations which, when using the oracle presented above, will lie the same face of the permutahedron. In Appendix D we provide an empirical comparison of these two operators when applied to the DAG learning problem.
4.2 DAG LEARNING
In order to select which edges to drop from R, we regularize the set of edge functions ( }/ to be sparse via S2 (( in eq. (2), i.e. || or ||] Incorporating the sparse decompositions of eq. (7) or eq. (8) into the original problem in eq. (2) yields
min Ea (0) 0.4
(x,, fo, (X (R)))) + (
(9)
where (Re) is the jth column of R and a² the top-k sparsemax or the SparseMAP distribution. Notice that for both sparse operators, in the limit т + 0_ the distribution a* puts all probability mass on one permutation: (0) (the sorting of 0), and thus eq. (9) is a generalization of eq. (2).
We can solve the above optimization problem for the optimal 0, 4 jointly via gradient-based opti- mization. The downside of this, however, is that training may move towards poor permutations purely because & is far from optimal. Specifically, at each iteration, the distribution over permutations a (0) is updated based on functional parameters & that are, on average, good for all selected permutations. In early iterations, this approach can be highly suboptimal it cannot escape from high-error local minima. To address this, we may push the optimization over & inside the objective:
min E Ba (0)
x fx por (Xo()))))
(10)
S.L. 4" (a) = arg min
Σ (x) (xo (R + x0
This is a bi-level optimization problem (Franceschi et al., 2018; Dempe & Zemkoho, 2020) where the inner problem fits one set of structural equations { 1 per (0). Note that, opposed to many other settings (e.g. in meta-learning,, the outer objective depends on only through the distribution a. (0), and not through the inner problem +. In practice, this means that the outer

---

optimization does not require gradients (or differentiability) of the inner solutions at all, saving computation and allowing for greater flexibility in picking à solver for fitting (a))).² For example, for ( Φ11 1 we may invoke any Lasso solver, and for () = ||ol we may use the algorithm of Louizos et al. (2017), detailed in Appendix B. The downside of this bi-level optimization is that it is only tractable when the support of a (0) has few permutations. Optimizing for & and o jointly is more efficient.
4.3 COMPUTATIONAL ANALYSIS
The overall complexity of our framework depends on the choice of sparse operator for learning the topological order and on the choice of the estimator for learning the edge functions. We analyze the complexity of learning topological orderings specific to our framework, and refer the reader to previous works for the analyses of particular estimators (e.g., Efron et al. (2004)). Note that, independently from the choice of sparse operator, the space complexity of DAGuereeotyee is least of the order of the edge masking matrix R. hence (dd)). This in line with most methods based on continuous optimization and can be improved by imposing additional constraints on the in-degree and out-degree of a node.
SparseMAP.
Each SparseMAP iteration involves length-d argsort and Cholesky update of s-by-s matrix, where s is the size of the active set (number of selected permutations), and by Carathéodory's convex hull theorem (Reay, 1965) can be bounded by d+ 1. Given a fixed number of iterations K (as in our implementation), this leads to a time complexity of (Kd²) and space complexity O(s² + sd) for SparseMAP. Furthermore, we warm-staan the sorting algorithm with the last selected permutation. Both in theory and in practice, this is better than the (dd) complexity of maximization over the Birkhoff polytope.
Top-k sparsemax. Complexity for top-k sparsemax is dominated by the complexity of the top-k oracle. our implementation, the top-k oracle has an overall time complexity O (K²dd) and space complexity (Kd)) (as detailed in Appendix A) when searching for the best K permutations. When K fixed, as in our implementation, this leads to an overall complexity of the top-k sparsemax operator (d2). In practice, K has to be of an order smaller than Vd for our framework to more efficient than existing end-to-end approaches.
4.4 RELATIONSHIP TO PREVIOUS DIFFERENTIABLE ORDER-BASHR METHODS
The advantages of method over Cundy et al. (2021); Charpentier et al. (2022) are: (a) parametrization. based on sorting, improves efficiency in practice (Appendix D) and has theoretically stabler learning dynamics as measured by our bound on SHD (Theorem (C.1)); (b) our method allows for any downstream edge estimator, including non-differentiable ones, critically allowing for off-the-shell estimators; (c) empirically our method vastly improves over both approaches in terms of SID and especially SHD on both real-world and synthetic data (Section 5).
EXPERIMENTS
5.1 EXPERIMENTAL SETUP
Datasets Reisach et al. (2021) recently demonstrated that commonly studied synthetic benchmarks have a key flaw. Specifically, for linear additive synthetic DAGs, the marginal variance of each node increases the 'deeper' the node is in the DAG (i.e., all child nodes generally have marginal variance larger than their parents). They empirically show that a simple baseline that sorts nodes by increasing marginal variance and then applies sparse linear regression matches or outperforms state-of-the-ae DAG learning methods. Given the triviality of simulated DAGs, here we evaluate all methods on two real-world tasks: Sachs (Sachs et al., 2005). dataset of cytometric measurements of phosphorylated protein and phospholipid components in human immune system cells. The problem consists of d = 11
This computational advantage shared with the score- fore-ffnniion estimator (SFE, Rubinstein, 1986; Williams, 1992; Paisley et al., 2012: Mohamed et al., 2020). We are not aware of any applications of SFE to permutation learning, likely due to the #P-completennss of marginal inference over the Birkhoff polytope (Valiant, 1979; Taskar, 2004).
7

---

sachs
syntren
180
SID
140
15
10
40
60
80
100
170
140
160
SHD
SHD
ours-spmax (linear)
outs-spmax-jointt (linear)
NPVAR (GAM)
Golem (NEV)
spraxspmax (non-linear)
surs-spmax-joont (non-linear)
VI-DP-DAG (softsort)
ours-spmax (LARS)
sortnregress
NoTears (linear)
BCONets
Figure 2: SHD vs SID on real datasets Sachs and SynTReN. Results are averaged over 10 seeds and the solutions lying the Pareto front circled.
nodes, 853 observations and of the graph reconstructed by Sachs et al. (2005) as ground-truit DAG, which contains 17 edges; Syn'ReN (den Bulcke et al., 2006), set of 10 pseudo-real transcriptional networks generated by the SynTRen simulator, each consisting of 500 simulated gene expression observations, and DAG of d = 20 nodes and of e edges with e E {20, ... 25]. We use the networks made publicly available by Lachapelle et al. (2020). In Appendix D we, however, compare different configurations of our method also on synthetic datasets, as they constitute an ideal test-bed for assessing the quality of the ordering learning step independently from the choice of structural equation estimator.
framework against state-of-the-an methods: NoTears (both its lin- ear (Zheng et al., 2018) and nonlinear (Zheng et al., 2020) models), the first continuous optimization method, which optimizes the Frobenius reconstruction loss and where the DAG constraint is enforced via the Augmented Lagrangian approach; Golem (Ng et al., 2020), another continuous optimiza- tion method that optimizes the likelihood under Gaussian non-equal variance error assumptions regularized by NoTears's DAG penalty: CAM (Bühlmann et al., 2014), two-stage approach that estimates the variable order by maximum likelihood estimation based on an additive structural equation model with Gaussian noise; NPVAR (Gao et al., 2020), an iterative algorithm that learns topological generations and then prunes edges based on node residual variances (with the Generalized Additive Models (Hastie & Tibshirani, 2017) regressor backend to estimate conditional variance); sortnregress (Reisach et al., 2021), 1) two-steps strategy that orders nodes by increasing variance and selects the parents of node among all its predecessors using the Least Angle Regressor (Efron et al., 2004); BCDNets (Cundy et al., 2021) and VI-DP-DAG (Charpentier et al., 2022), the two differentiable, probabilistic methods described in Section 2. Before evaluation, post- pest-process the graphs found by NoTears and Golem by first removing all edges with absolute weights smaller than 0.3 and then iteratively removing edges ordered by increasing weight until obtaining DAG, the learned graphs often contain cycles.
Baselines We benchmark
Metries We compare the methods by two metrics, assessing the quality of the estimated graphs: the Structural Hamming Distance (SHD) between true and estimated graphs, which counts the number of edges that need be added removed or reversed to obtain the true DAG from the predicted one: and the Structural Intervention Distance (SID, Peters & Bühlmann, 2015), which counts the number of causal paths that are broken in the predicted DAG. is standard in the literature to compute both metrics, given their complementarity: SHD evaluates the correctness of individual edges, while SID evaluates the preservation of causal orderings. We further remark here that these metrics privilege opposite trivial solutions. Because true DAGs usually sparse (i.e., their number of edges is much smaller than the number possible ones), SHD favors sparse solutions such as the empty graph. On the contrary, given topological ordering, SID favors dense solutions (complete DAGs in the limit) as they are less likely to break causal paths. For this reason, we report both metrics and highlight the solutions the Pareto front in Figure 2.
and training details We set the hyper- parameters of all methods to their default values. For method tuned them by Bayesian Optimization based the performance, in terms
8

---

Legend:
True Edge
Correct prediction
Missing Edge
Wrong prediction
True DAG
Daguerretyppe (ours)
sortnnegreess
VI-DP-DAG
Figure 3: Sachs. True DAG and DAGs predicted by the best-performing methods. We plot on the left of the bar correct and missing edges and on the right of the bar wrong edges found by each method. DAGuerreotype strikes a good balance between SID and SHD; other methods focus overly on one over the other by either predicting too few (sortnregress)
many edges (VI-DP-DAG).
of SHD and SID, averaged over several synthetic problems from different SEMs. For our method, we optimize the data likelihood (under Gaussian equal variance error assumptions, derived in Ng et al. (eq. 2, 2020)) and we instantiate / to masked linear function (linear) or as a masked MLP as for -nTeinernonn.narr. Because our approach allows for modular solving the functional parameters &, we also experiment with using Least Angle Regression (LARS) (Efron al., 2004). We additionally apply 12 regularizations on {0, Φ to stabilize training, and we standardize all datasets to ensure that all variables have comparable scales. More details provided in Appendix D. The code for running the experiments is available at https: /ttpt//gg//tt/////// com wwtp:///waateeeeccc///eeeeeeeeee
5.2 RESULTS
We report the main results in Figure 2, where we omit some baselines (e.g., DAGuerreotype with sparseMAP) for sake of clarity. We present the complete comparison in Appendix D. together with additional metrics and studies on the sensitivity DAGuerreotype depending on the choice of key hyper-parameters. To provide better idea of the learned graphs, we also plot in Figure the graphs Icarned by the best- pest-performing methods on Sachs. We observe that the solutions found by NPVAR and the matrix-exponente in terms of SHD, but have the worst SIDs. This can be explained by the high sparsity of their predicted graphs (see the number of edges in Tables and 2). On the contrary, the high density of the solutions of VI-DP-DAG makes them among the best in terms of SID and the worst in terms of SHD, DAGuerreotype provides solutions with good trade-off between these two metries and which lie on the Pareto front. Inevitably its performance strongly depends on the choice of edge estimator for the problem hand. For instance, the linear estimator better suited for Sachs than for SynTReN. In Appendix D, we assess our method's performance independently from the quality of the estimator with experiments on synthetic data where the underlying SEM is known, and the estimator can be chosen accordingly.
regularized methods (NoTears and Golem) are the best
CONCLUSION, LIMITATIONS AND FUTURE WORK
In this work, we presented DAGuerreotype, a permutation-basee method for end-to-end learning of directed acyclic graphs. While our approach shows promising results in identifying DAGs, the optimization procedure can still be improved. Alternative choices of estimators, such as Generalized Additive Models (Hastie & Tibshirani, 2017) as done in CAM and NPVAR, could be considered to improve identification. Another venue for improvement would be to include interventional datasets at training. would be interesting to study this context whether our framework is more sample- efficient, i.e., allows to learn the DAG with fewer interventions observations.

---

ACKNOWLEDGEMENTS
We are grateful to Mathieu Blondel, Caio Corro, Alexandre Drouin and Séhastien Paquet for discus- sions. Part of this work carried out when VZ was affiliated with INRIA-London and University College London. Experiments presented in this paper were partly performed using the 5000 testbed, supported by scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see ht https: / www. grid5000. fI). VN acknowledges support from the Dutch Research Council (NWO) project VI. Veni.212.2222
5000
REFERENCES
Raghavendra Addanki, Shiva Kasiviswanathan, Andrew McGregor, and Cameron Musco. Efficient intervention design for causal discovery with latents. In International Conference on Machine Learning, pp. 63-73. PMIR, 2020.
Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery data mining, pp. 2623-2631, 2019.
Bryon Aragam and Qing Zhou. Concave penalized estimation of sparse gaussian bayesian networks. The Journal of Machine Learning Research, 16, 2015.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: learning dags via m- m-matrices and log-determinant acyclicity characterization, Advances in Neural Information Processing Systems, 2022.
Mathieu Blondel, Olivier Tehoul, Quentin Berthet, and Josip Djolonga. Fast differentiable sorting and ranking. In International Conference on Machine Learning, pp. 950-959, PMLR, 2020.
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. In NeurlPS, 2020.
Peter Bühlmann, Jonas Peters, and Jan Ernest. Cam: Causal additive models, high-dimensinn order search and penalized regression. The Annals of Statistics, 42(6))2526 25526, 2014.
Bertrand Charpentier, Simon Kibler, and Stephan Günnemann. Differentiable DAG sampling. In International Conference on Learning Representations, 2022.
David Maxwell Chickering. Learning bayesian networks np-complete. In Doug Fisher and Hans- Joachim Lenz (eds.), Learning from Data Fifth International Workshop Artificial Intelligence and Statistics, AISTATS 1995, Key West, Florida, USA, January, 1995. Proceedings. Springer, 1995.
Gonçalo Correia, Vlad Niculae, Wilker Aziz, and André Martins. Efficient marginalization of discrete and structured latent variables via sparsity. Advances in Neural Information Processing Systems, 33, 2020.
Chris Cundy, Aditya Grover, and Stefano Ermon Bed nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34, 2021.
James Cussens. Bayesian network learning with cutting planes. In Uncertainty in Artificial Intelli- gence, 201
Stephan Dempe and Alain Zemkoho. Bilevel optimization: Advances and next challenges. Springer, 2020.
Tim Van den Bulcke, Koen Van Leemput, Bart Naudts, Piet van Remortel, Hongwu Ma, Alain Verschoren, Bart De Moor, and Kathleen Marchal. Syntren: a generator synthetic gene expression data for design and analysis of structure learning algorithms. BMC Bioinform., 7:43, 2006.
Frederick Eberhardt. Causation and intervention. Unpublished doctoral dissertation, Camegie Mellon University, pp. 93, 2007.
10

---

Bradley Efron, Trevor Hastie, lain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of statistics, 32(2):407-499, 2004.
Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning, pp. 1568-1577, PMLR, 2018.
Nir Friedman and Daphne Koller. Being bayesian about network structure. A bayesian approach to structure discovery in bayesian networks. Machine learning, 50, 2003.
Ming Gao, Yi Ding, and Bryon Aragam. polynomial-time algorithm for learning nonparametric
causal graphs. In Advances in Neural Information Processing Systems, 2020.
Dan Geiger and David Heckerman. Learning gaussian networks. Uncertainty Proceedings 1994, pp. 235-243, Elsevier, 1994.
Godfrey Harold Hardy, John Edensor Littlewood, and György Pólya. Inequalities. Cambridge University Press, 1952.
Trevor Hastie and Robert Tibshirani. Generalized additive models. Routledge, 2017.
Alain Hauser and Peter Bühlmann. Two optimal strategies for active learning of causal models from interventional data. International Journal of Approximate Reasoning, 55(4):926-939, 2014.
Yue He, Peng Cui, Zheyan Shen, Renzhe Xu, Furui Liu, and Yong Jiang. DARING: differentiable causal discovery with residual independence. In KDD, 2021.
Jean Kaddour, Aengus Lynch, Qi Liu, Matt J. Kusner, and Ricardo Silva. Causal machine learning: A survey and open problems. arXiv preprint arXiv:2206. 150.: 2022.
Diviyan Kalainathan and Olivier Goudet. Causal discovery toolbox: Uncover causal relationships in python. arXiv preprint arXiv: (rXXiv:900.02222 920.00278.. 2019.
Neville Kenneth Kitson, Anthony C. Constantinou, Zhigao Guo, Yang Liu, and Kiattikun Chobtham. A survey of bayesian network structure learning. CoRR, abs/2109. 1415, 2021.
Murat Kocaoglu, Alex Dimakis, and Sriram Vishwanath. Cost-optimal learning causal graphs. In International Conference Machine Learning, pp. 1875-1884. PMLR, 2017.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural DAG learning. In ICLR, 2020.
Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclicity constraints. In International Conference on Learning Representations, 2022.
Christos Louizos, Max Welling, and Diederik P Kingma. Learning sparse neural networks through 1_0 regularization. arXiv preprint arXiv: 171201171.. 2017.
Gonzalo Mena, David Belanger, Scott Linderman, and Jasper Snoek. Learning latent permutations with gumbel-sinkhom networks. In International Conference on Learning Representations, 2018.
Elad Michael, Tony A. Wood, Chris Manzie, and Iman Shames. Global sensitivity analysis for the linear assignment problem. In ACC, pp. 3387-3392. IEEE, 2020.
Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih. Monte carlo gradient estimation in machine learning. J. Mach. Learn. Res., (1(132) 2020.
Ignavier Ng. AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear dags. In NeurlPS, 2020.
Ignavier Ng, Sébastien Lachapelle, Nan Rosemary Ke, Simon Aacoste-Julien. and Kun Zhang On the convergence of continuous constrained optimization for structure learning. In AISTATS, Proceedings of Machine Learning Research, 2022.

---

Vlad Niculae, Andre Martins, Mathieu Blondel, and Claire Cardie. Sparsemap: Differentiable sparse structured inference. In International Conference on Machine Learning, pp. 3799-3808. PMLR, 2018.
Mathias Niepert, Pasquale Minervini, and Luca Franceschi. Implicit mle: Backpropagating through discrete exponential family distributions. Advances in Neural Information Processing Systems, 34, 2021.
Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 1999.
John Paisley, David M. Blei, and Michael I. Jordan. Variational bayesian inference with stochastic search. In Proc. ICML, 2012.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché-Buu E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 8024-8035 Curran Associates, Inc., 2019.
Max Paulus, Dami Choi, Daniel Tarlow, Andreas Krause, and Chris Maddison. Gradient estimation with stochastic softmax tricks. Advances in Neural Information Processing Systems, 33:5691 15770, 2020.
Judea Pearl. Causality: Models, reasoning and inference. 2000.
Jonas Peters and Peter Bühlmann. Neural computation, 27(3):771-799, 2015.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Sebastian Prillo and Julian Eisenschlos. Softsort: A continuous relaxation for the argsort operator. In
ICML, 2020.
Joseph D. Ramsey, Madelyn Glymour, Ruben Sanchez- Romero, and Clark Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. Interna- tional Journal of Data Science and Analytics, 3, 2017.
John Reay. Generalizations of theorem of Carathéodory. Number 54. American Mathematical Soc., 1965.
Alexander G Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated dag! varsortability in additive noise models. Advances Neural Information Processing Systems, 34, 2021.
Paul Rolland, Volkan Cevher. Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In International Conference on Machine Learning, pp. 18741-18753. PMLR, 2022.
Reuven Y Rubinstein. The score function approach for sensitivity analysis of computer simulation models. Mathematics and Computers in Simulation, 28(5):351-37,, 1986.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721): 523-529, 2005. doi: 10.1126/science. 105809.
Andrew D Sanford and Imad A Moosa. A bayesian network structure for operational risk modelling in structured finance operations. Journal of the Operational Research Society, ((4)):33-444, 2012.
Mauro Scanagatta, Cassio P de Campos, Giorgio Corani, and Marco Zaffalon. Learning bayesian networks with thousands of variables. Advances Neural Information Processing Systems, volume 28, 2015.
12

---

Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. The Annals of Statistics, 48(3):1514-1533 2020.
Karthikeyan Shanmugam, Murat Kocaoglu, Alexandros G Dimakis, and Sriram Vishwanath. Learning causal graphs with small interventions. Advances in Neural Information Processing Systems, 28, 2015.
Ajit P Singh and Andrew W Moore. Finding optimal Bayesian networks by dynamic programming. Citeseer, 2005.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Chandler Squires, Yuhao Wang, and Caroline Uhler. Permutation-basee causal structure learning with unknown intervention targets. In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), 2020.
Ben Taskar. Learning structured prediction models: A large margin approach. PhD thesis, Stanford University, 2004.
Mare Teyssier and Daphne Koller. Ordering Ordeeing-based search: A simple and effective algorithm for learning bayesian networks. In UAI, pp. 548-549. AUAI Press, 2005.
Leslie G Valiant. The complexity of computing the permanent. Theor. Comput. Sci., 8(2): 189-211, 1979.
Ronald Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine Learning. 8(3-4):229-255 1992.
Jing Xiang and Seyoung Kim. A lasso for learning sparse bayesian network structure for continuous variables. In Advances in Neural Information Processing Systems, volume 26, 2013.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In ICML, Proceedings of Machine Learning Research, 2019.
Bin Zhang, Chris Gaiteri, Liviu- Gabriel Bodea, Zhi Wang, Joshua McElwee, Alexei A Podtelezhnikov, Chunsheng Zhang, Tao Xie, Linh Tran, Radu Dobrin, et al. Integrated systems approach identifies genetic nodes and networks late-onset alzheimer's disease. Cell, 153(3):007-720, 2013.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Dags with NO TEARS: continuous optimization for structure learning. In Advances in Neural Information Processing Systems, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, pp. 3414-3425, PMLR, 2020.
13

---

A TOP-K ORACLE
In this section, we propose an efficient algorithm for finding the top-k scoring rankings of a vector and prove correctness. We leverage this result in order apply the sparsemax operator (Correia et al., 2020) to our DAG learning problem. We are not aware of existing works on this topic and believe this result is of independent interest.
A.1 NOTATION AND DEFINITIONS
Let us denote the objective function go(o) (0. evaluating the quality a permutation a, and denote one of its maximizers as a! arg max Io( corresponding to the argsort of 9. Notice that we can equivalently write this objective as (0) = (0", .p) by applying the permutation to 0. Definition A.1 (Top-k permutations). We denote by (+) C Ea sequence of K highest-scoring permutations according to ge i.e., (0) = {o', ,0" ,.....o oh) and for any of & Tx(0):
g 2 go .... go(o*) 2 200(d')
In this section, we will make of the following definition and lemma. Definition A.2 (Adjacent transposition). oj := (j j+1) denotes 2-cycle of the components of the vector o, which transposes (flips) the two consecutive elements σg and Ø Lemma A.1. Let a E End be a permutation. Exactly one of the following holds.
1. go(a j) - go(o")
2. There exists an adjacent transposition (+]+ such that aj satisfies go(o) > go(o"))
Proof. Denote by 0" the permutation of 0 by a and let J = E [d-1] 03 > 14 If If then 0" is in increasing order, by the Rearrangement Incquality (Hardy et al., 1952, Thms. 368-369) a I-best permutation. Otherwise, applying the adjacent transposition (/ j++1) increases the score:
go j) - (o = 07 0f+1 00
A.2 BEST-FIRST SEARCH ALGORITHM
Algorithm (2) finds the set of k-best scoring permutations given 0. Starting from an optimum of ge. the algorithm grows the set of candidate permutations P(0) by adding all those that are one adjacent transposition away from the ith-best permutation at iteration i. It then selects the best scoring permutation P(0) be the (ip-(i+1) solution. Theorem A.2 proves that an i+ 1th-best solution must lie in this set P(O), hence that Algorithm (2) is correct.
Theorem A.2 (Correctness of Algorithm (2)), Given a sequence of k-l-hest permutations Th-1 _0',0? .o.....a"" - there exists k-th best gt i.e., one satisfying (o(a*=1) > gola*)) 2 (o) any of & (0), with the property that o² = o'j for some i kk....,k-1] k-1] and j fl,...,,d11 d-I..
Proof. Let h be a k-th best permutation.
Case 1. go(o * (o")). Invoking Lemma A.1 we have go( > (o( Since the inequality strict, must have otj E (0).
Case 2. go (a* = (²*) In this case, we have at least & permutations tied for first place. Any two permutations with equal score can only differ in indices that correspond to ties in 0. Therefore, any two tied permutations are connected by a trajectory of transpositions of equal score. This is a face of the permutahedron, of size > K, containing (0). This face must contain at least one permutation that is one adjacent transposition away from one of T (0), and we may as well take this one the k-th best instead of *
14

---

Algorithm 2: Top-k permutations.
Data: k E (1,...,d)) d 0 E Rd
Result: top-k permutations (0)
P(O) + {or ER arg ma, (0(0)) /* initialize set of candidates
with an optimum
while [T) < do
ER arg () (0) / retrieve best scoring solution
among the candidates that has not been retrieved yet P(O) + (aj
/+ add its one adjacent
transposition away permutations to the candidated
(a( < (x) U {o} /* update set of best permutations
*/
end
A.3 COMPUTATIONAL ANALYSIS
Finding 7 requires sorting vector of size d (hence Oo d) complexity). Then, at each iteration Algorithm (2), the maximum among the best candidates P(0) needs to found, which requires (ddk) complexity P(0)) < (d 2)k, and (dd) adjacent flip operations are applied to it.
The most expensive operation is checking that the selected best candidate not already in Th- (0). At worst, this requires going through the whole P(0) and comparing them in order to all top-k solutions, so no more than x P(0) comparisons of cost O(d) each.
This leads overall time complexity of dK). The space complexity is dominated by the size of P(0), which contains at most Kdd d vectors of size d, leading to O(Kd).
B LO REGULARIZATION
In the linear and non-linear variants of DAGuerreotype, we implement the regularization term Ωę of the inner problem of Eq. (10) with approximate LO regularizer. The exact L0 norm
|wll
1 wyyk0
(11)
counts the number of non-zero entries of w IR and, when used regularizer, it favors sparse solutions without injecting other priors. However, its combinatorial nature and its non-differentiability reparameterize (11) by
make its optimization intractable. Following Louizos et al. (2017), introducing a set of binary variables z E (0, 1) 1) and letting w = z, sothat ||w/|o Σ 24. Next, reformulate
let p(z; π) Bernoulli(zr) where E [0, 1j For linear SEMs, the Inner Problem (10) as follows:
min Σ (aj-plss(m))
(x, (M (x) (x
1,)) +
Tji
(12)
where now the decision (inner) variables are intended matrices and { 2 0 is hyperparameter. In the non-linear (MLP) case, we achieve sparsity at the graph level by group-regularizing the parameters corresponding to each input variable. In the experiments, we optimize (12) using the sae-s Monte Carlo straight-thoough estimator and set the final functional parameters
(T- (m=
W W. MAP : w. H π
(13)
where H is the Heaviside function. We leave the implementation of more sophisticated strategies, such as the relaxation with the hard-concrete distribution presented in (Louizos et al., 2017) or other estimators (e.g. Paulus al., 2020; Niepert et al., 2021), to future work.
15

---

1.5
1.0
0.5
0.0
-0.5
1.5
1.0
-1.0
0.5
-1.5
0.0
0.5
1.5
1.0
0.5
11.
0.0
1.0
0.5
1.0
-1.5
1.5
Figure 4: Representation of the degeneracy hyper-planes in R33 and of two-parameter points (in red) whose connecting segment intersects two hyyer-planes.
C CHARACTERIZATION AND SENSITIVITY OF THE VECTOR PARAMETRIZATION
In this section, we provide some intuition into the behavior of our score vector parameterization in the space of complete DAGs. More precisely, look the Maximum A Posterioni (MAP) complete DAG
M
Dold] where o² E arg max (0,ps)
(14)
and how it varies as a function of the score vector E Rd. The permutation 0" is MAP state (or mode) of both the SparseMAP and the sparsemax distributions (as well the standard categori- cal/softmax distribution). For brevity, shall rename M ( in the following.
We choose to work on the space of complete DAGs because there a one-to-one correspondence between topological orderings and complete DAGs. This is not generally true when analyzing the space of DAGs, as permutation does not uniquely identify a DAG and vice versa.
C.1 PRELIMINARY
For the results of this section, we will use the following relationship between transpositions (adjacent not) and SHD.
Proposition C.0.1 (SHD difference after flip). Consider 0 Rd and O' which is obtained by applying a flip (i j) with the convention that < j. All the edges from nodes between i and directed towards i or need to be reversed. Thus, the SHD between complete DAGs SHD SHD(M(0), M ((0)) 2(j i) - 1. (because no new undirected edges are added, no undirected edges are removed, and 2 edges are reversed).
From Proposition (C.0.1) deduce that the SHD difference after applying adjacent flip is (H(
C.2 ANALYSIS
Recall that o* sorts the elements of in increasing order. Then, the points E Rd where arg (0, Po)) is non-singleton (degeneracy) exactly those that have least tie among their
16

---

entries (i.e., 3 i, such that 0, =,)). Following this simple observation, we can populate Rd with hyper-planes (of dimensionality d - 1), and call them Hi,; with i < J. such that A0 E H. 0 = 0 Intersections of such hyper-planes are lower dimensional subspaces where more than 2 entries of 0 equal. For d = 3, Figure 4 depicts the (2) 3 hyper-planes in orange, green and blue, and their intersection (a line) in gray.
The hyper-planes [HH.shis 's delimit exactly d! open cones of the form C = {0 E R4 < 0g < < 07 a E Ed}. i.., each cone is the space of all points that give the same MAP permutation and do not contain any ties.
This partition of the space allows us to reason about the sensitivity of our MAP estimates to changes to the parameters. For instance, the points of a cone do not have any ties, changes to the score vector do not affect SHD: 40 E C, O' E C (HD(M(O). M(0)) = 0.
We first analyze the sensitivity to single-entry changes, hence assessing how much entry of 0 can be individually perturbed without changing its MAP. Proposition C.0.2 provides the entry-wise ranges in which the optimal solution set does not change.
Proposition C.O.2 (Entry-wise intervals). For any Rd and Et E. arg max (0" .p) arg max (0 + e;0s, P) if and only if:
E 12, I2,....dd 11. 8, E (0 - 02", 0, 02¹),
= L, (-00, 0,, - -7"):
ii=d, (0 07',00) 00)
where denotes the i-th standard unit vector and a* E arg max (0, p²).
We can relate this result to changes in terms of SHD, by making the following observation: if we gradually increase one coordinate 0; initially ranked of, its rank changes in the optimal ordering soon it becomes greater than the coordinate right after it in the ordering, entailing adjacent transposition between the two coordinates? Greater perturbations entail a longer sequence of adjacent transpositions. By Proposition C.O.1, this implies an increase in SHD) of for each transposition. Visually, the SHD increases by every time the perturbed vector crosses a hyper-plane along its perturbation direction, this corresponds to swapping two components that are adjacent when optimally sorted.
For comparison, we can apply the same reasoning the matrix parametrization of the linear assign- ment problem, deployed by Cundy et al. (2021) for learning permutation matrices. In this context, we can leverage the edge sensitivity analysis reviewed in Michael al. (2020, Equations 6-7). As a general remark, it is not as intuitive to determine the entry-wise intervals for this problem, not only because we have d2 variables instead of d but especially because it requires solving a different assignment problem per entry. Furthermore, the minimal perturbation that changes the MAP does not necessarily result in flipping two adjacent nodes, entailing an SHD relative to the optimal permutation of at least 1. Consider, for example, the following 3 x 3 matrix parameter:
16
16
15
5
16
10
16
9
10
Its optimal permutation matrix corresponds to the rank (3, 1, with score of 29. The minimal perturbations on 11,3 or 23,2 that change the MAP result in the solution (2, 1,3) which has an SHD of 3 w.r.t. the optimal (and the second best score of 31). For problems of scale larger than this example, the resulting SHD can take values up to 2d - 3 for non-adjacent transpositions.
We generalize and formalize this relationship between perturbations and SHD changes for general vector perturbations (global sensitivity). Theorem (C.1) upper bounds the SHD between complete DAGS of any pair of score vectors by the number of hyper-planes crossed by the segment connecting them.
AA similar reasoning also applies when decreasing the value of one coordinate instead.
17

---

Theorem C.1 (Global sensitivity). For any O E R2 and O' E Rd
(M)( - xE
SH -0) (0''00--0) dt
(15)
where (x(2) is the (generalized) Dirac delta that evaluates to infinity if x E A and otherwise.
Proof. Let first consider the case where 0 C and O' C", i.e., the two vectors do not contain ties. Let us denote (respectively o') the permutation that sorts the components of point of C' (C") by increasing order. The minimal-length sequence of adjacent flips that need to be applied to a to obtain of has length equal to the number of times the segment connecting their parameters crosses a hyper-plane. Then the SHD between their complete DAGs equals the minimal number of adjacent flips, which proves the result. When either A or l' lies on hyper-plane, the minimal required number of adjacent flips might be smaller, hence the upper bound in Theorem C.I.
In Figure the segment connecting the two red dots intersects the green and blue hyperplanes and hence the resulting complete DAGs will have an SHD of at most 2 (in fact, exactly 2 in this case).
This intuitive characterization links the SHD distance in the complete DAG space a partition of R resulting from the "sorting" operator (i.e. the MAP, or maximizer of the linear program) of the score vector O. This implies that, during optimization, if 0₂ are the interior of any cone (which happens almost surely) then it is very likely that updating the parameters results in small changes to the SHD unless the parameters have all similar values.
A deeper analysis of the vector parameterization is the object of future work, as we hope it can fuel further improvements in the optimization algorithm, such as better initialization strategies or reparameterization. For instance, optimizing in S (over polar coordinates) would expose the role of the radius r as similar to the temperature parameter for the distributions we consider (the smaller the radius, the higher the "temperature"). Typically, the temperature is left constant during training annealed, suggesting this might be advantageous also in our scenario. We leave the exploration of this strategy to future work.
D ADDITIONAL EXPERIMENTS
We provide detailed description of the experimental setup and report additional results. The method is implemented in (PyTorch, Paszke et al., 2019), and the code used for carrying out the experiments included in the supplementary material. All experiments were run on machine with 16 cores, 32Gb of RAM and an NVIDIA
10-SXM4-80GG GPU.
DAGuerreotype 's optimization and evaluation For our method, we optimize the data likelihood (under Gaussian equal variance eтror assumptions, as derived in Ng et al. (Equation 2, 2020)) We optimize the bilevel Problem (10) when not specified otherwise. We report results for the following three variants of DAGuerreotype, where the edge estimator of the graph is instantiated
(linear) with LO regularization, 1 (X, Ay) = (eb o A2) and wy E IRd;
(non-linear) with LO regularization, 170h (X, Ay) - hy 9 (X, with h2 locally connected
MLP with one hidden layer, 50 hidden units and sigmoid activation function, 33 linear layer with d x 50 hidden units (50 per parent i) and masking out all non-parents of j (according Aj), and Φbf+ |1
2 ll2 as for NoTears (non-linear);
(LARS) with Least Angle Regression (LARS) (Efron al., 2004), 15 (X, X (Φ)oA,) and Φ; E IRd.
We additionally apply a lz regularization on {0, Φ) to stabilize training, and we standardize all datasets to ensure that all variables have comparable scales.
With any variant, the outer problem is optimized for 5, 000 maximum iterations by gradient descent and early-stopped when approximate convergence is reached. When using the (linear) and (non-linear)
18

---

sachs
spmax (linear)
52
ours-spmax (non-linear)
ours-spmax
50-
ours spmax joint (linear)
ours jamax-join::: (non-linear)
40
ours-SpMAP (linear)
oUrs-SpMAP (nor-lineat)
(LARS)
Ours-spMAP-joint (linear)
ours-spMAAP-jin (non- linear)
sortnregress
42
NPVAR (GAM)
CAM
10-
No Tears (linear)
NoTears (non neear)
Galem (NEV)
VI-P-DAA (softsort)
VI-DP DAG (sinkhorn)
15
SHD
40
BCDN
Figure 5: SHD vs SID on Sachs. The solutions lying on the Pareto front are colored in orange and the others in blue.
back-ends, the graph of each permutation is optimized also by gradient descent, for 1, 000 epochs and also with an carly-stopping mechanism based on approximate convergence. After training, the graph for the mode permutation further fine-tuned, and the final evaluation is carried with this model.
We also experiment with the approximated, but faster, joint optimization of 0 and Φ in Problem (2) instead of the bilevel formulation, when we add the suffix joint the method name. In this case, as we need a differentiable graph estimator for jointly updating the ordering and the graph, we instantiate the method only with the (linear) and (non-linear) back-ends These models are optimized by gradient descent for 5, 000 epochs and with an early-stopping mechanism based on approximate convergence.
The default hyper-parameters of our methods were chosen as follows. We set the sparse operators' temperature
1 and K = 100, the strength of the la regularizations to 0.0005, and tuned the learning rates for the outer and inner optimization [10~6, 1, The tuning was carried out by Bayesian Optimization using (Optuna, Akiba et al., 2019) for 50 trials on synthetic problems, consisting of data generated from different types of random graphs (Scale-Free, Erdös Réniii BiPartite) and of noise models (e.g., Gaussian, Gumbel, Uniform) with 20 nodes and 20 or 40 expected edges. For each setting, three datasets are generated by drawing a DAG, its edge weights uniformly in [-2, 0..5 U [0.5, 2). and 1, 000 data points. A of hyper-paramettrs is then evaluated by averaging its performance on all the generated datasets, and the tuning is carried minimize SHD and SID jointly. The default value of tlper-parameter was then set be the average value among those lying on the Pareto front and rounded up to have single significant digit.
[10- 10-1] and pruning strength À E
Baseline optimization and evaluation
All baseline methods are optimized using the codes released by their authors, apart from CAM that is included in the Kalainathan & Goudet (Causal Discovery Toolbox, 2019). Before evaluation, we post- process the graphs found by NoTears and Golem by first removing all edges with absolute weights smaller than 0.3 and then iteratively removing edges ordered by increasing weight until obtaining DAG, as the learned graphs often contain cycles. For the probabilistic baselines (VI-DP-DAC and BCDNets), we make use of the mode model (in particular, the mode permutation matrix) for evaluation.
We set the typer-parametery of all methods to their default values, released together with the source code, apart from the parameters of the Least Angle Regressor module of sortnregress that uses the Bayesian Information Criterion for model selection.
Additional results on real-world tasks In Figures 5, 6. and Tables and 2 extend the evalu- ation on real-world tasks provided in the main text. More precisely, we report the results also for DAGuerreotype with sparseMAP, and provide additional metrics for comparison: the FI score using the existence of an edge as the positive class and the number of predicted edges to assess the density of the solutions.
19

---

syntren
sura-spmax (linear)
ours: sps-xppax (non linear)
ours-spmax (LARS)
180
ours-spmax-jint (linea))
purs-spmaxxjint (non-lineat)
ours: SDMP (linear)
ours-spMoc (non-lincar)
ours-spMAP (LARS)
ours-spMMP-joint. (linear)
curs curs-spNAP-joint (non-linear)
sortnregress
NPVAR (GAM)
CAM
(linea))
NoTbaars (non linear)
Golom INEV)
VI-O--DAG (sofsort)
V--P-DDG (sinkhorn)
46
60
6b
SHD
100
125
140
160
BCDNets
Figure 6: SHD SID on SynTReN. The solutions lying the Pareto front are colored in orange and the others in blue.
Sparsemax vs sparseMAP comparison In Figures 8 and we report an analysis of the effect of DAGuerreotype 's hyper-parameter K and of the sample size n on the quality of the learned DAG. Recall that K corresponds to the maximal number of selected permutations for sparsemax and to the maximal number of iterations of the active set algorithm for sparseMAP, and that is the principal parameter that controls the computational cost of the ordering learning step in our framework. This analysis is carried out data generated by linear SEM from scale-free graph with equal variance Gaussian noise. We choose this simple setting for two reasons: the true DAG can be identified from (enough) observational data only (Proposition 7.5, Peters al., 2017); provided with the true topological ordering, the LARS estimator can identify the true edges. In this setting, we can then assess the quality of sparsemax and sparseMAP independently from the quality of the estimator. For reference, in Figures 8 and 7 we also report the performance of optimizing LARS with a random ordering or with a true one.
We observe that sparsemax and sparseMAP provide MAP orderings that are significantly better than random ones for K 2 and for any sample size. For K big enough, these orderings give DAGS that are almost good as when knowing the variable ordering. Furthermore, apart when K 2, increasing the sample size generally results an improvement (although moderate) in the performance of sparsemax and sparseMAP. However, the gap from true's performance does not reduces when increasing n or K, which can be explained by the non-convexity of the search space and DAGuerreotype getting stuck in local minima. We further observe that sparsemax generally provides better solutions than sparseMAP's at comparable training time. The only settings where this is not the case are for d = 30 and K >> 35. Notice that sparseMAP's performance peaks at K = 35 and degrades for higher K in this setting. This phenomenon can be due to the inclusion of unnecessary orderings
MPP'seMPP's set, which ends up hurting training.
We further study in Figures 9 and 10 the effect of the pruning strength, (controlled by X) on the perfor- mance of both operators on the real datasets. For this experiment, we instantiate DAGuerreotype with the linear estimator and train it by joint optimization. As a general remark, when strongly penalizing dense graphs (higher X) SHD generally improves and SID degrades. On Sachs, the two operators do not provide significantly different results, while on SynTReN we find that sparsemax provides better SHD for comparable SID.
Additional results on synthetic data We report in Figures 11 and an additional comparison of DAGuerreotype with several sfate-of-the-and baselines on synthetic problems of varying number of nodes d and = 1, 000 samples generated from scale-free DAGs with 2d expected number of edges and different noise models: (Gaussian) linear SEM with equal variance Gaussian noise; (Gumbel) linear SEM with equal variance Gumbel noise; (MLP) 2-layer neural network SEM with sigmoid activations and equal variance Gaussian noise. To limit the varsortabilityy of the generated problems, the parameters of the SEMs uniformly drawn from [-0.5, -0.1] U 0.5). The resulting problems are still varsortable on average, as demostrated by the great performance of sortnregress, and by the fact that by initializing DAGuerreotype parameters 0 with the marginal variances of
20

---

Table 1: Sachs. We report Structural Hamming Distance (SHD, the lower the better), Structural Interventional Distance (SID, the lower the better), (F1, the higher the better), and the number of predicted edges for all methods.
Method
SHD 4
SID 4
F1 f
edges
NoTears (linear)
16.0
52.0
0.095
NoTears (non-linear)
14.0
50.0
0.320
Golem (NEV)
15.0
51.0
0.190
VI-DP-DAG (softsort)
43.0
37.0
0.265
51
VI-DP-DAG (sinkhorn)
42.0
38.0
0.269
50
BCDNets
39.0
46.0
0.196
sortnregress
13.1
50.8
0.377
9.0
CAM
30.0
47.0
0.28
33.0
NPVAR (GAM)
13.2
51.2
0.344
9.6
DAGuerreotype -ppaaaa (linear)
13.8
51.6
0.323
7.7
SSMMP (linear)
13.8
49.1
0.309
7.6
DAGuerreottpe--ppaa (LARS)
13.8
51.7
0.316
9.5
DAGuerreooype-spMAPP (LARS)
14.4
49.6
0.275
9.2
DAGuerreotype-pppax (non-linear)
13.5
50.9
0.348
7.1
DAGuereeotype-spMAA (non-linear
14.1
48.1
0.32
7.4
DAGuerreatype -spmaxx-oooo (linear)
16.9
44.9
0.36
16.1
(linear)
14.1
51.5
0.244
6.8
DAGuerreotype -pmaxxx-ooos (non-linear)
14.7
50.7
0.265
8.6
DAGuerreotype-spMMPPPooo (non-linear)
14.4
51.1
0.236
5.0
Table 2: Syn7ReN. We report Structural Hamming Distance (SHD, the lower the better). Structural Interventional Distance (SID, the lower the better), Topological Ordering Pearson Correlation (TOPC, the higher the better), (F1, the higher the better) number of predicted edges all averaged over the 10 networks.
Method
SHD
SID 4
F1
edges
NoTears (linear)
32.4
184.0
0.157
17.7
NoTears (non-linear)
40.0
187.9
0.165
28.1
Golem (NEV)
30.5
191.7
0.152
15.4
VI-DP-DAG (softsort)
165.4
114.8
0.109
175.7
VI-DP-DAG (sinkhorn)
164.0
121.4
0.104
173.6
BCDNets
117.8
148.6
0.113
119.0
sortnregress
86.4
156.0
0.151
89.8
CAM
82.6
116.0
0.192
86.6
NPVAR (GAM)
36.1
191.0
0.184
26.0
DAGuerreotype-spe -spmax (linear)
51.4
171.0
0.182
44.5
DACherreeoyppe (linear)
51.1
161.0
0.211
47.2
DAGuerreotype -spmax (LARS)
77.9
103.0
0.238
86.0
DAGuereeotpee-pp... (LARS)
78.7
116.0
0.212
83.8
DAGuencotyyee-spnag (non-linear)
62.8
121.0
0.245
65.1
DAGuerreotype-ssp.... (non-linear)
61.9
134.0
0.255
65.5
DAGuemeettype-sppaxx (linear)
152.0
138.0
0.104
161.0
A.Guerreottpe-spMAPP... (linear)
122.0
131.0
0.139
131.0
Satterreeotppe-spajxjjjoo (non-linear)
154.0
90.9
0.149
169.0
DAGuereeotyee-pMAPPP (non-linear)
138.0
118.0
0.143
150.0
21

---

120
1350
187
1000
SHD
27
150
750-
5e0-
250
random
48
192
187
200
100
12
0.0
04-
4.2
3000
15000
time (3)
592
2000
100000
802
1000
permulattinn aatons
100
60-
80
20-
250
1400-
e0gess
600
660-
50
K
Figure 7: Comparison of different strategies for learning topological orderings, on data (1, 000 samples) generated from linear SEM with Scale-Free graph and Gaussian noise, and varying number of nodes d. In order from top to bottom, we plot SHD, SID, F1, training time in seconds, number of permutations, and number of predicted edges (all the end of training) as function of the sparse operators' parameter K, which corresponds to the maximal number of sampled permutations for sparsemax and the maximal number of iterations the active set algorithm for SparseMAP. We also include two simple variants of DAGuerreotype, where the ordering is fixed to one true ordering (true) or to a random one (random). Results are averaged over 10 seeds.
the nodes consistently improves its performance, compared to initializing them with the vector of all zeros. For these experiments we set バ 10, use the linear edge estimator and jointly optimize all DAGuerreotype 's parameters.
Compared other differentiable order-based methods, DAGuerreotype consistently provides signif- icantly better trade-off between SHD and SID, confirming our findings on the real-world data. Indeed, these baselines generally discover DAGs with high false positive rates. DAGuerreotype equipped with the sparseMAP operator also improves upon the linear continuous methods based the exponential matrix regularization, but when equipped with the top-k sparsemax operator its results on these settings depend good initialization of 0 (the marginal variances in this case) and worsen with the number of nodes. A higher value of l would be required to improve DAGuerreotype -sparsemax performance in these settings, as shown in Figure 8.
22

---

K+2
K.110
KF 35
K-o 100
SND
300
sample:522
sample 5122
sample S22
samplesits
Figure 8: Comparison of different strategies for learning topological orderings, samples of varying size (n [100, 5'000] on the z-axes) generated from a linear SEM with Scale-Free graph and Gaussian noise, and number of nodes d = 20. In order from top to bottom, we plot SHD, SID, F1, training time in seconds, number of permutations, and number of predicted edges (all at the end of training) for 4 values of the sparse operators' parameter K, which corresponds the maximal number of sampled permutations for sparsemax and to the maximal number of iterations of the active algorithm for SparseMAP. We also include simple variants of DAGuerreotype. where the ordering is fixed to one true ordering (true) or to random one (random). Results are averaged over 10 seeds.
SHD
SID
edges
50
32-
sparsemax sparsel66P
20
30
pruning
pruning
pruning A
pruning
Figure 9: Sachs. Effect of I.0 pruning intensity (controlled by on SHD, SID, and number of learned edges for DAGuerreotype jointly optimizing linear estimator and the topological ordering distribution either with sparsemax (sparsemax) or sparseMAP (sparseMAP).
In terms of running times, DAGuerreotype is aligned with NoTears and is generally faster than CAM, Golem and VI-DP-DAG. Of course DAGuerreoyype 's running times strongly depend on the value of K, the choice of edge estimator and the optimization of either the joint or bi-level problems.
23

---

edges
240
297
125
sparsomax
192
160
sparseMAP
50
102
0.00
18.9
20-9
pruning
pruning
pruning
pruning
Figure 10: SynTReN. Efffect of LO pruning intensity (controlled by A) on SHD, SID, Fl and number of learned edges for DAGuerceotype jointly optimizing a linear estimator and the topological ordering distribution either with sparsemax (sparsemax) sparseMAP (sparseMAP).
Gaussian d-10
Gaus sian d=20
Gaussian d=30
Gaussian d=100
236-
250
250-
350
230
13c-
1000
50
50
300
200
2000
SHD
SHD
SHD
SHD
Gumbel d-10
Gumbel d0-20
Gumbell d=30
Gumbel d=100
5000
50
4000
SI0
40
335-
2000
310
2000
30-
132
10
20
10
60
103
15c
160
201
130
420
1800
2000
3003
4200
SHD
SHD
SHD
SHD
d=10
MLP dm
MLP dm30
MLP d=100
2000
355
2500
335-
2000
25
230
2500
20
130
1000
20
10
25
30
106
250
302
230
300
1000
2000
2000
4000
5000
SHD
SHD
SHD
SHD
curs sparsemax zeros ours wpwrsepax-varraae.e.eme inceeees
curs sparseMtP varlances WH-DPDDG (sottsort)
Golem (EV)
X
porthrogress
No tears (liear))
CAM
(pinkhom))
Figure 11: SHD SID synthetic datasets generated from scale-free DAGs with Gaussian (top), Gumbel (middle) and MLP (bottom) SEMs. DAGuerreotype 's (ours) 0 is either initialized with zero vector (zeros) with the marginal variances (variances).
24

---

(S)
0.30
68-
3:00
2000
2300
3300
230
0-
1.-
68-
Gumbbol
cur SgorseMAP-varaanee
4300
VI DAS (softsort)
VINDNAA Sankhumi
(incar)
sottrragiiss
CAM
kn
SHD
SID
time (s]
Chussonn
4000
2200
68-
3000
2000
84
Cuumbel
8.00
3300
3000
2000
2000
200
ours-paasemma vartreeeevvvrtrrre
our--paa Mallaaaallllreeencccc
4300
66-
Vi (softsort)
2.00
1.00
20
20
SHD
SID
time (5)
2000
Coussian
2000
1100
1000
1000
Gumbel
2000
1000
outs-parreemaa-zers
Golem (EV)
2000
Nolkers (ineer)
MLP
2000
Figure 12: Comparison DAGuerreotype (ours) with related methods terms of SHD, SID, F1, training time (in seconds) synthetic datasets generated from scale-free DAGs with Gaussian, Gumbel and MLP SEMs. DAGuerreotype (ours) is either initialized with zero vector (zeros) with the marginal variances (variances). For of reading, split the full comparison (top) into two, to focus on the comparison with differentiable order-based methods (middle) and with differentiable methods based on the matrix exponential constraint (bottom).
25


Reviews:
[{"review_id": "sauj5H9suA", "reviewer": "Reviewer_bSZG", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper studies the problem of learning DAGs from observational data. Different from prior work under the continuous framework, the authors propose an approach to optimize over the polytope of permutations, where edges can be jointly or conditionally optimized given a particular topological ordering. The key to the approach is to leverage recent progress on sparse relaxation methods such as top-k sparsemax and sparseMAP. Some experiments on real data are provided.", "summary_of_the_review": "My main concerns are on about the experiments, it is not clear if the proposed approach works well for linear/nonlinear models, for different dimensions, and for different types of graphs. ", "strength_and_weaknesses": "Strengths: I generally find the approach interesting and a nice mix of ideas from the sparse relaxation literature to the DAG learning problem. The paper is also generally well-written modulo some unspecified notation. \n\nWeaknesses: The main weakness is the experimental section. Given that this is not a theoretical work, one should at least expect a more comprehensive set of experiments. The work of Reisach et al. (2021) is cited for skipping comprehensive synthetic experiments basically, and both real-world networks experimented on have a small number of nodes. The scaling problem when using OLS has been known for a while for linear regression, and hence the same is expected to happen for linear SEMs, this is why methods such as PC are robust to this problem since PC works by performing CI tests---Reisach et al. (2021). To alleviate the marginal variance problem, one can just use a bit smaller weights for the simulations. It would be great to see how the proposed method behaves against other baselines for ER, SF graphs for different numbers of nodes. Moreover, it is a must to compare the runtimes as well, given that other methods such as NOTEARS do not have a clear computational complexity, it is important to see what the actual performance is in terms of runtime. \n\nFrom the experiments, I would like to conclude: Is the proposed method better? Is it faster? Both? I cannot conclude that from the current set of experiments.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "sMka2zpVwrW", "reviewer": "Reviewer_kPDZ", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "The paper develops a continuous optimization method for learning DAGs by optimizing over the polytope of permutation vectors. Two versions of the methods are provided: (1) alternately iterates between node-ordering and edge-optimization, and (2) optimizes them jointly. The resulting method has several advantages over existing methods, such as (1) returning exact DAGs and (2) accommodating non-differentiable black-box estimators. Experiment results on two real-world datasets are provided.", "summary_of_the_review": "- The proposed method may be rather incremental compared to existing ones.\n- Lack of experiment results on synthetic benchmarks.", "strength_and_weaknesses": "Strength: \n- A new continuous approach to traverse the DAG space is provided. Solving continuous DAG learning task is a topic with significance.\n- The approach is flexible and accommodates non-differentiable black-box estimators.\n\nWeakness:\n- As the paper mentioned, the possibility of using differentiable approach to learn DAGs via ordering has been proposed by (Cundy et al., 2021; Charpentier et al., 2022). Therefore, the proposed method may seem incremental, i.e., optimizing over permutation vectors instead of permutation matrices as in existing works.\n    - I would suggest the authors to provide a detailed discussion of the advantages of their method over (Cundy et al., 2021; Charpentier et al., 2022). Also, can those existing methods be modified to accommodate non-differentiable black-box estimators?\n- The paper argued that existing works (Cundy et al., 2021; Charpentier et al., 2022) required expensive optimization over the permutation matrices. However, the paper currently lacks empirical studies to support it--either comparison of running time or scalability to very large graphs would be more compelling to demonstrate the efficiency/scalability of the proposed method.\n-  Although the paper mentioned a recent study by Reisach et al. (2021) regarding the flaw of synthetic dataset in DAG learning benchmark, I would still suggest the authors to compare their method to the baselines on those synthetic benchmarks to have an idea about how the method performs. This is because for real world datasets, the model assumption (e.g.,. linear or nonlinear additive noise models) may not hold.\n- In 'over-relaxation' part of Section 1, several methods based on penalizing the 'DAG-ness' of the adjacency matrix are actually *guaranteed to return a valid DAG* under some conditions, as shown by the following papers. Therefore the statement is not entirely true and should be corrected.\n    - Remark 2 of Ng et al. (2022). On the Convergence of Continuous Constrained Optimization for Structure Learning.\n    - Lemma 6 of Bello et al. (2022). DAGMA: Learning DAGs via M-matrices and a Log-Determinant Acyclicity Characterization.\n\n----------------\nI appreciate the detailed response from the authors, which has addressed my concern. I have increased my score accordingly.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "bCU3pCUPkyQ", "reviewer": "Reviewer_EZMX", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "In this paper, the authors look at the problem of recovering the DAG structure from samples from a DAG-structured distribution. Despite having seen a lot of important works spanning the last several decades, the problem has remained open and is an active area of research. The proposed solution in this paper consists of two steps i) finding a topological ordering of the DAG, ii) discovering edges according to the ordering. The authors manage to do both steps which is end-to-end on the one hand and accomodates any black-box edge discovery algorithm for step ii).\n\nExperiments conducted by the authors reveal that their algorithm are better in terms of both the metrics SHD and SID. In fact, their solutions learn at the Pareto frontier of these two metrics.", "summary_of_the_review": "Given the superior experimental performance of the proposed approach, I tend to accept the paper.", "strength_and_weaknesses": "Strength: the proposed solution is flexible and more efficient and more accurate than the existing methods. Comparison with the related work has been covered in great detail. The experimental evidence is extensive, consisting of the best algorithms known before and includes both benchmark and synthetic datasets.\n\nWeakness: Ultimately, the paper reads like a conglomeration of existing ideas. I could not clearly find a single new technique that the authors have found.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "WdrQlQMX-N", "reviewer": "Reviewer_SU9Y", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "In this paper, the authors propose a framework for DAG learning from data. The method consists of a two-step procedure, where the first one is based on an observation about optimizing the order of nodes over the permutahedron. Then the edges to keep are selected by using a regularization over the set of functions.\nThe method is tested on two real-world datasets with promising results.\n", "summary_of_the_review": "The authors of the paper propose an interesting two-step framework for DAG learning, where different methods can be used for the sub-tasks.\nIt is well written, and the results are promising.", "strength_and_weaknesses": "The strengths of the paper are the clarity, the observation of the optimization over the permutahedron, and the fact that it is actually a framework where one can choose different methods for the sub-tasks.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster