CoLiDE: Concomitant Linear DAG Estimation
Authors: Seyed Saman Saboksayr
Gonzalo Mateos
Mariano Tepper
Venue: ICLR 2024 poster (2024)
Keywords: directed acyclic graph
concomitant scale estimation
causal discovery
graph structure learning
Primary area: general machine learning (i.e., none of the above)
TL;DR: We bring to bear ideas from concomitant scale estimation in linear regression and propose a novel score function for regression-based inference of linear DAGs, demonstrating significant improvements relative to existing state-of-the-art methods.

Abstract: We deal with the combinatorial problem of learning directed acyclic graph (DAG) structure from observational data adhering to a linear structural equation model (SEM). Leveraging advances in differentiable, nonconvex characterizations of acyclicity, recent efforts have advocated a continuous constrained optimization paradigm to efficiently explore the space of DAGs. Most existing methods employ lasso-type score functions to guide this search, which (i) require expensive penalty parameter retuning when the $\textit{unknown}$ SEM noise variances change across problem instances; and (ii) implicitly rely on limiting homoscedasticity assumptions. In this work, we propose a new convex score function for sparsity-aware learning of linear DAGs, which incorporates concomitant estimation of scale and thus effectively decouples the sparsity parameter from noise levels. Regularization via a smooth, nonconvex acyclicity penalty term yields CoLiDE ($\textbf{Co}$ncomitant $\textbf{Li}$near $\textbf{D}$AG $\textbf{E}$stimation), a regression-based criterion amenable to efficient gradient computation and closed-form estimation of exogenous noise levels in heteroscedastic scenarios. Our algorithm outperforms state-of-the-art methods without incurring added complexity, especially when the DAGs are larger and the noise level profile is heterogeneous. We also find CoLiDE exhibits enhanced stability manifested via reduced standard deviations in several domain-specific metrics, underscoring the robustness of our novel linear DAG estimator.

Published as a conference paper at ICLR 2024
COLIDE: CONCOMITNT LINEAR DAG ESTIMATION
Seyed Saman Saboksayr, Gonzalo Mateos
Mariano Tepper
University of Rochester
Intel Labs mariano.. wwwpereinteee.commmm
{ssaboksa, gmateosb) www.eesseee..oorooere.e......
rwwhester. edu
ABSTRACT
We deal with the combinatorial problem of learning directed acyclic graph (DAG) structure from observational data adhering to a linear structural equation model (SEM). Leveraging advances in differentiable, nonconvex characterizations of acyclicity, recent efforts have advocated a continuous constrained optimization paradigm efficiently explore the space of DAGs. Most existing methods em- play lasso-type score functions to guide this search, which (i) require expensive penalty parameter retuning when the unknown SEM noise variances change across problem instances; and (ii) implicitly rely on limiting homoscedasticity assump- tions. In this work, we propose new convex score function for sparsity-awan learning of linear DAGs, which incorporates concomitant estimation of scale and thus effectively decouples the sparsity parameter from the exogenous noise levels. Regularization via smooth, nonconvex acyclicity penalty term yields CoLiDE (Concomitant Linear DAG Estimation), a regression-basse -based criterion amenable efficient gradient computation and closed-form estimation of noise variances methods
in heteroscedastic scenarios. Our algorithm outperforms state-of-the-ee without incurring added complexity, especially when the DAGs are larger and the noise level profile is heterogeneous. We also find CoLiDE exhibits enhanced stability manifested via reduced standard deviations several domain-specinf metrics, underscoring the robustness of our novel linear DAG estimator.
INTRODUCTION
Directed acyclic graphs (DAGs) have well-appreciated merits for encoding causal relationships within complex systems, as they employ directed edges to link causes and their immediate effects. This graphical modeling framework has gained prominence in various machine learning (ML) appli- cations spanning domains such as biology (Sachs al.) [0055] Lucas al. (004). genetics (Zhang al. 2013), finance (Sanford & Moosa 2012), and economics Pourret et al. (2008), to name few. However, since the causal structure underlying group of variables is often unknown, there a need to address the task of inferring DAGs from observational data. While additional interventional data are provably beneficial to the related problem of causal discovery (Lippe et al.) 2021 Squires et al.) 2020 Addanki al. 2020 Brouillard al. 2020), said interventions may be infeasible or ethically challenging to implement. Learning a DAG solely from observational data poses significant com- putational challenges. primarily due to the combinatorial acyclicity constraint which is notoriously difficult to enforce (Chickering 1996] Chickering et al. 2004). Moreover, distinguishing between DAGs that generate the same observational data distribution nontrivial. This identifiability chal- lenge may arise when data are limited (especially in siggh-dmensional settings). or, when candidate graphs in the search space exhibit so-termed Markov equivalence; see e.g., (Peters et 2017).
Recognizing that DAG learning from observational data is in general an NP-hard problem, recent ef- forts have advocated continuous relaxation approach which offers an efficient means of exploring the space of DAGs (Zheng et al., 2018 Ng et al. 2020) Bello et al.) (022). In addition to break- throughs in differentiable, nonconvex characterizations of acyelicity (see Section 3. the choice of an appropriate score function is paramount to guide continuous optimization techniques that search for faithful representations of the underlying DAG. Likelihood- based methods (Ng et al., 2020) enjoy desirable statistical properties, provided that the postulated probabilistic model aligns with the actual causal relationships (Hastie et al. 2009 Casella & Berger 202I). On the other hand, regression-based methods (Zheng et al. 2018 Bello al. 2022) exhibit computational efficiency,

---

Published as a conference paper at ICLR 2024
robustness, and even consistency (Loh & Bühhmann 2014), especially in high figh-dimensiona settings where both data scarcity and model uncertainty are prevalent (Ndiaye et al. 2017). Still, workhorse regression- based methods such as ordinary least squares (LS) on the assumption of homoscedas- ticity, meaning the variances of the exogenous noises are identical across variables. Deviations from this assumption can introduce biases in standard error estimates, hindering the accuracy of causal discovery (Long & Ervin 2000; Loh & Bühlmann 2014). Fundamentally, for heteroscedastic linear Gaussian models the DAG structure is non-identifiable from observational data (Peters al. 2017).
Score functions often include a sparsity-promoting regularization, for instance an fn-ormm penalty as in lasso regression [Tibshirani] 1996). This turn necessitates careful fine-tuning of the penalty pa- rameter that governs the trade-off between sparsity and data fidelity (Massias et al.) (2018). Theoreti- cal insights suggest an appropriately scaled regularization parameter proportional to the observation noise level (Bickel et al. 2009), but the latter is typically unknown. Accordingly, several linear re- gression studies (unrelated to DAG estimation) have proposed convex concomitant estimators based on scaled LS, which jointly estimate the noise level along with the regression coefficients (Owen) 2007] Sun & Zhang 2012 Belloni et al., 2011 Ndiaye et al. 2017). A noteworthy representative is the smoothed concomitant lasso (Ndiaye al. 2017), which not only addresses the aforementioned parameter fine-tuning predicament but also accommodates heteroscedastic linear regression scenar- ios where noises exhibit non-equal variances. While these challenges are not extraneous DAG learning, the impact of concomitant estimators is so far largely unexplored in this timely field.
Contributions. In this work, we bring to bear ideas from concomitant scale estimation in (sparse) linear regression and propose a novel convex score function for regression-band inference of linear DAGs, demonstrating significant improvements relative to existing state-of-the-art methods Our contributions can be summarized as follows:
We propose a new convex score function for sparsity- aawar learning of linear DAGs, which in- corporates concomitant estimation of scale parameters to enhance DAG topology inference using continuous first-order optimization. We augment this LS-based score function with a smooth, non- convex acyclicity penalty term to arrive at CoLiDE (Concomitant Linear DAG Estimation), simple regression- based criterion that facilitates efficient computation of gradients and estimation of exoge- nous noise levels via closed-form expressions. To the best of our knowledge, this is the first time that ideas from concomitant scale estimation permeate benefits to DAG learning.
In existing methods, the sparsity regularization parameter depends on the unknown exogenous noise levels, making the calibration challenging. CoLiDE effectively removes this coupling, leading to minimum (or no) recalibration effort across diverse problem instances. Our score function is fairly robust to deviations from Gaussianity, and relative to ordinary LS used in prior score-based DAG estimators, experiments show CoLiDE is more suitable for challenging heteroscedastic scenarios.
We demonstrate CoLiDE's effectiveness through comprehensive experiments with both simulated (linear SEM) and real-worle datasets. When benchmarked against tht DAG learning algorithms, CoLiDE consistently attains better recovery performance across graph ensembles and exogenous noise distributions, especially when the DAGS are larger and the noise level profile is heterogeneous. We also find CoLiDE exhibits enhanced stability manifested via reduced standard deviations in various domain-specifie metrics, underscoring the robustness of our novel estimator.
2 PRELIMINARIES AND PROBLEM STATEMENT
Consider a directed graph g EV.E, W), where V = f1.....d d} represents the set of vertices, and 3 V x V is the set of edges. The adjacency matrix W = [ wal E Raxd collects the edge weights. with Waj * 0 indicating a direct link from node i to node j. We henceforth assume that g belongs to the space ITD of DAGs, and rely on the graph to represent conditional independencies
among the variables in the random vector x = [x1...., xa]] T E RRd. Indeed, if the joint distribution p(x) satisfies a Markov property with respect to G e D, each random variable z, depends solely on
'hile preparing the final version of this manuscript, we became aware of independent work in the un- published preprint (Xue et al. 2023). which also advocates estimating (and accounting for) exogenous noise variances to improve DAG topology inference. However, dotears is two-step marginal estimator requiring additional interventional data, while our concomitant framework leads to new L.S -S-hased score function and, importantly, relies observational data only.
2

---

Published as a conference paper at ICLR 2024
its parents PA, E V W # 0}. Here, we focus linear structural equation models (SEMs) to generate said probability distribution, whereby the relationship between each random variable and its parents is given by z; w x + 24. Vi E V. and 12,,....: zd] _ is vector of mutually independent, exogenous noises; see e.g. (Peters et al. 2017). Noise variables 2; can have different variances and need not be Gaussian distributed. For dataset X Rddx of i.i.d. samples drawn from P(x), the linear SEM equations can be written in compact matrix form as X W X+2 Z.
Problem statement. Given the data matrix X adhering a linear SEM, our goal is to learn the latent DAG G E DD by estimating its adjacency matrix W as the solution to the optimization problem
min
S(W) subject to g(W) E D, W
where S(W) is a data-dependent score function to measure the quality of the candidate DAG. Ir- respective of the criterion, the non-convexity comes from the combinatorial acyclicity constraint g(W) E D; see also Appendix for a survey of combinatorial search approaches relevant to solv- ing I but less related to the continuous relaxation methodology pursued here (see Section
proper score function typically encompasses a loss or data fidelity term ensuring alignment with the SEM as well as regularizers to promote desired structural properties on the sought DAG. For linear SEM, the ordinary LS loss 1 XX- W X is widely adopted, where || is the Frobenius 2n
norm. When the exogenous noises are Gaussian, the data log-likelihood can be effective alterna- tive (Ng et al. 2020). Since sparsity is a cardinal property of most real- world graphs it is prudent to augment the loss with an fi-norm regularizer to yield S(W) 1 X W X I|? + X [W]l Il2, where 2 0 is a tuning parameter that controls edge sparsity. The score function S(W) bears resemblance to the mulli-task variant of lasso regression (Tibshirani] (996). specifically when the response and design matrices coincide. Optimal rates for lasso hinge on selecting A x / d/n Bickel et al., 2009 Li et al. 2020). However, the exogenous noise variance o² is rarely known in practice. This challenge is compounded in heteroscedastic settings, where one should adopt a weighted LS score (see (Loh & Bühlmann] 2014) for an exception unlike most DAG learning meth- ods that stick to bias-inducing ordinary LS). Recognizing these limitations, in Section
2n
| we propose a novel LS-based score function to facilitate joint estimation of the DAG and the noise levels.
RELATED WORK
We briefly review differentiable optimization approaches that differ in how they handle the acyclicity constraint, namely by using an explicit DAG parameterization or continuous relaxation techniques.
Continuous relaxation. Noteworthy methods advocate an exact acyclicity characterization using nonconvex, smooth functions H : IR R + R of the adjacency matrix, whose zero level set is D. One can thus relax the combinatorial constraint G(W) E D by instead enforcing H(W) = and tackle the DAG learning problem using standard continuous optimization algorithms (Zheng et al.,
2018 Yu et al. 2019 Wei al. 2020 Bello al., (022). The pioneering NOTEARS formulation introduced Hexpm
d. where denotes the Hadamard (element-wise) product and Tr(-) is the matrix trace operator ZZhnng et al. 2018). Diagonal entries of powers of W o W encode information about cycles in G. hence the suitability of the chosen function. Follow-up work suggested more computationally efficient acyclicity function Hpply( (W) = Tr( WwW.WW/.. W where I is the identity matrix (Yu et al. 2019;; see also (Wei et al. 2020) that studies the general Recently, Bello et al. (2022) proposed the
Tr(e" WoW)
family H((W) Σd_ C₁ Tr((W> w)d). Cg log-determinant acyclicity characterization Hidel (1, s) dloog(s) - (oness)))oo((e(((ii s1 - Wo o WW), s RR outperforming prior relaxation methods in terms of (nonlinear) DAG recovery and efficiency. Although global optimality results are so far elusive, progress is being made (Deng et al. 2023b).
Order-based methods. Other recent approaches exploit the neat equivalence G G(W) E W IT UII, where II E 1) dxd a permutation matrix (essentially encoding the variables' causal ordering) and U Radd is an upper-triangular weight matrix. Consequently, one can search over exact DAGs by formulating an end-to-end differentiable optimization framework to minimize S(II UII) jointly over II and U, or in two steps. Even for nonlinear SEMs, the challenging pro- cess of determining the appropriate node ordering is typically accomplished through the Birkhoff polytope of permutation matrices, utilizing techniques like the Gumbel-Sinkhun approximation (Cundy et al. 202T), or, the SoftSort operator (Charpentier al. 2022). Limitations stemming from
3

---

Published as a conference paper at ICLR 2024
misaligned forward and backward passes that respectively rely hard and soft permutations
well documented (Zanteeesch et al. 2023). The DAGuerreotype approach in (Zantedeschi et al.)
2023 instead searches over the Permutahedron of vector orderings and allows for non-smooth edge weight estimators. TOPO (Deng et al. 2023a) performs a bi-level optimization, relying on topolog- ical order swaps the outer level. Still, optimization challenges towards accurately recovering the causal ordering remain, especially when data are limited and the noise level profile is heterogeneous.
While this work is framed within the continuous constrained relaxation paradigm to linear DAG learning, preliminary experiments with TOPO (Deng et al. 2023a) show that CoLiDE also benefits order-based methods (see Appendix E.8
We leave a full exploration as future work.
CONCOMITANT LINEAR DAG ESTIMATION
Going back our discussion lasso-type score functions in Section minimizing S(W) 2 +I
W W X + A W subject to smooth acyclicity constraint (W) = 0 (as in c.g., NoTears (hhngg al. 2:18)): (i) requires carefully retuning when the unknown SEM noise vari- ance changes across problem instances; and (ii) implicitly relies on limiting homoscedasticity as- sumptions due to the ordinary LS loss. To address issues (i)- here we propose new convex score function for linear DAG estimation that incorporates concomitant estimation of scale. This way, we obtain a procedure that is robust (both in terms of DAG estimation performance and parameter fine- tuning) to possibly heteroscedastic exogenous noise profiles We were inspired by the literature of concomitant scale estimation in sparse linear regression (Owen 2007,) Belloni et al. 2011 Sun & Zhang 2012: Ndiaye et al.] 2017); see Appendix A for concomitant lasso background that informs the approach pursued here. Our method is dubbed CoL.DDE (Concomitant Linear DAG Estimation).
Homoscedastic setting. We start our exposition with a simple scenario whereby all exogenous vari- ables 21, 2d in the linear SEM have identical variance o². Building on the smoothed concomi- tant lasso (Ndiaye et al., (007), we formulate the problem of jointly estimating the DAG adjacency matrix W and the exogenous noise scale a as
min 1 2nd
do
X
W
X
+ A// subject to H(W) O,
Wo0
2
S(W.r)
where H Rdxdd_ + R is a nonconvex, smooth function, whose zero level set is D as discussed in Section 3 Notably, the weighted, regularized LS score function S(W, σ)) is now also a function
of and it can be traced back to the robust linear regression work of Huber (198I). Due to the rescaled residuals, A in (2) decouples from a as minimax optimality now requires À x /log n (Li al. 2020 Belloni et al. 201T). A minor tweak to the argument in the proof of (Owen 2007 Theorem 1) suffices to establish that S(W, .r) is jointly convex in W and a. Of course, (2) is still a nonconvex optimization problem by virtue of the acyelicity constraint H(W) = Huber [1981)
included the term do 2 so that the resulting squared scale estimator is consistent under Gaussianity. The additional constraint a 2 To safeguards against potential ill -l-posed scenarios where the estimate ā approaches zero. Following the guidelines in (Ndiaye et al. 2017), we set
[ x 10-2.
Vdn
With regards to the choice of the acyclicity function, we select HiderWW,s
d llg(s)
log(det(s1 - W o W)) based on its more favorable gradient properties in addition to several other compelling reasons outlined in (Bello et al. 2022 Section 3.2). Moreover, while LS-based lin- ear DAG learning approaches are prone to introducing cycles in the estimated graph, it was noted that a degermeennnnnt term arising with the Gaussian likeeiiiddood objective tends to mitigate this undesirable effect (Ng et al. 2020). Interestingly, the same holds true when Haaer is chosen larizer, but this time without being tied to Gaussian assumptions. Before moving on to optimization issues, we emphasize that our approach is in principle flexible to accommodate other convex loss functions beyond LS (c.g., Huber's loss for robustness against heavy-tailed contamination), other acyclicity functions, and even nonlinear SEMs parameterized using e.g., neural networks. All these are interesting CoLiDE generalizations beyond the scope of this paper.
regu-
Optimization considerations. Motivated by our choice of the acyclicity function, to solve the constrained optimization problem (2) we follow the DAGMA methodology by Bello et al. (2022). Therein, it is suggested to solve sequence of unconstrained problems where Hute is dualized and
4

---

Published as a conference paper at ICLR 2024
viewed as regularizer. This technique has proven more effective in practice for specific problem as well, when compared to, say, an augmented Lagrangian method. Given a decreasing sequence of values Mk 0, at step & of the COLIDE-EV (equal variance) algorithm one solves
do + +|IWil + Hudar W 88),,
min W.02%
j4ke
X
W
X
(3)
2no
where the schedule of hyperparameters Mr 2 0 and s& 0 must be prescribed prior to implemen- tation; see Section [4.1] Decreasing the value of μu enhances the influence of the acyclicity function Hase(W (W, s) in the objective. Bello et al. (022) point out that the sequential procedure (3) is remi- niscent of the central path approach of barrier methods, and the limit of the central path 14 - is guaranteed yield DAG. In theory, this means additional post-processing (e.g., edge trimming) is needed. However, in practice we find some thresholding is required to reduce false positives.
Unlike Bello et al. (0022), CoLiDE-EV jointly estimates the noise level and the adjacency matrix W for each P4g- To this end, one could solve for a in closed form [cf. (4) and plug back the solution in (3), obtain DAG-regularized square-root lasso (Belloni et al.) 2011) type of problem in W. We did not follow this path since the resulting loss |X- W X I is not decomposable across samples, challenging mini- mini-batch based stochastic optimization if it were needed for scalability (in Appendix B we show that our score function is fully decomposable and present corresponding preliminary experiments). A similar issue arises with GOLEM (Ng et al., 2020), where the derived Gaussian profile likelihood yields non-separable log-sum loss. Alternatively, and similar to the smooth concomitant lasso algorithm Ndiaye et al. (2017), we rely on (inexact) block coordinate descent (BCD) iterations. This cyclic strategy involves fixing to its most up-to-date value and minimizing (3) inexactly w.r.t. W, subsequently updating a in closed form given the latest W via
Ô = max (VIr(L-W) W) cov(X) ()) /d,o0)
(4)
where cXv(X) == _XX is the precomputed sample covariance matrix. The
interplay between noise level and DAG estimation should be apparent. There are several ways to inexactly solve the W subproblem using first-order methods. Given the structure of (3). an elegant solution is to rely on the proximal linear approximation. This leads to an ISTA-type update for W, and overall a provably convergent BCD procedure in this nonsmooth, nonconvex setting (Yang eet al. 2020 Theorem 1). Because the required line search can be computationally taxing, we opt for a simpler heuristic which is
single step of the ADAM optimizer (Kingma & Ba 2015) to refine W. We observed that running multiple ADAM iterations yields marginal gains, since we are anyways continuously re-updating W in the BCD loop. This process is repeated until either convergence is attained, or, a prescribed maximum iteration count T3 is reached. Additional details on gradient computation of (3) w.r.t. W and the derivation of (4) can be found in Appendix B
Heteroscedastic setting. We also address the more challenging endeavor of learning DAGs in het- eroscedastic scenarios, where noise variables have non-equal variances (NV) of ....... o. Bringing
to bear ideas from the generalized concomitant multi-task lasso (Massias et al. 2018] and mimicking the optimization approach for the EV case discussed earlier, we propose the CoLiDE-NV estimator
1 1 Tr ((X - W W X w w + Tr(E) + +||WII + Hiae 888k). 2n.
min W W.EZE.
(4k
(5) Note that Σ = diag(o 01................... aa) is a diagonal matrix of exogenous noise standard deviations (hence not a covariance matrix). Once more, we set En = / Vdiag (cov (X)) x 10-2, where v(-) is meant to be taken element-wise. closed form solution for Σ given W is also readily obtained,
max (Vdiag ((I WTT cov(X)(I W)), Σo
(W)(
(6)
A summary of the overall computational procedure for both CoLiDE variants is tabulated under Algorithm see Appendix B for detailed gradient expressions and a computational complexity dis- cussion. CoLiDE's per iteration cost is O(d²), ), on par with of-te- thette-artttttt DAG learning methods.
The first summand in (5) resembles the consistent weighted LS estimator studied in Loh & Bühlmann (2014), but therein the assumption is that exogenous noise variances are known up to
5

---

Published as a conference paper at ICLR 2024
Algorithm 1: CoLiDE optimization
Function CoLiDE-EV update:
In: data X and hyperparameters and H = [(1K, 8n, Tx))
Update W with one iteration of
Out: DAG W and the noise estimate (EV) or Σ (NV). Compute lower-hounds de Eg.
first-order method for (
Compute using
Initialize = 0, a on x 102 or Σ Eg x 103.
Function CoLiDE-NV update:
foreach (1uR, Sk; T&) E H do
Update W with one iteration of
for ( I..... Ti do
first-order method for 5
LApply CoLiDE-EV or NV updates using puu and sk.
Compute Σ using
a constant factor. CoL removes this assumption by jointly estimating W and Σ, with marginal added complexity over finding the DAG structure alone. Like GOLEM (Ng et al. 2020) and for general (non-identifiable) linear Gaussian SEMs, as n - oc CoLiDE-NV probably yields DAG that quasi-equivalent to the ground truth graph (see Appendix for futher details).
4.1 FURTHER ALGORITHMIC DETAILS
Initialization. Following Bello et al. (022). we initialize W = 0 as it always lies in the feasibility region of Hidet- We also initialize σ = σg x 102 in CoLiDE-EV or Σ = Eo 10a for CoLiDE-NV.
Hyperparameter selection and termination rule. To facilitate the comparison with the state-of- the-art DAGMA and to better isolate the impact of our novel score function, we use the hyperparam- eters selected by Bello et al. (2022) for CoLiDE-EV and NV. Hence, our algorithm uses À 0.05 and employs K = 4 decreasing values of μa E [1.0, 0. 0.01, 0,001], and the maximum number of 10º, 7 x 10"]. Furthermore, early stop- ping is incorporated, activated when the relative error between consecutive values of the objective function falls below 10-6 The learning rate for ADAM is 3 x 10-4
BCD iterations is specified Ta [2 x 10 2 x 10
We adopt 8 [1, 0.9, 0.8, 0.7].
Post- post-processing Similar to several previous works (hhng et al. 2018 Ng et al. 2020 Bello et al.) 2022), conduct final thresholding step to reduce false discoveries. In this post-processing step, edges with absolute weights smaller than 0.3 are removed.
5 EXPERIMENTAL RESULTS
We now show that CoLiDE leads to statc-of-the-ar DAG estimation results by conducting a compre- hensive evaluation against other state-of of-te-of-thheart approaches: GES (Ramsey al.) 2017), GOLEM [gg et al. 2020), DAGMA (Bello et al. (002). SortNRegress (Reisach et al. 2021) ). and DAGuerreo- type (Zantedeschi et al. 2023). We omit conceptually important methods like No Tears (Zheng al., 2018), that have been outclassed in practice. To improve the visibility of the figures, we report in each experiment the results of the best performing methods, excluding those that perform particu- larly poorly. We use standard DAG recovery metrics: (normalized) Structural Hamming Distance (SHD), SHD over the Markov equivalence class (SHD-C), Structural Intervention Distance (SID), True Positive Rate (TPR), and False Discovery Rate (FDR). We also assess CoLiDE's ability to esti- mate the noise level(s) across different scenarios. Further details about the setup are in Appendix [D]
DAG generation. As standard (Zheng al. 2018), we generate the ground truth DAGS utilizing the Erdös-Rényi or the Scale-Free random models, respectively denoted as ERK or SFk, where i the average nodal degree. Edge weights for these DAGs drawn uniformly from a range of feasible edge weights. We present results for & = 4, the most challenging setting (Bello et al. 2022).
Data generation. Employing the linear SEM, simulate n 1000 samples using the homo- or heteroscedastic noise models and drawing from diverse noise distributions, i.c., Gaussian, Exponen- tial, and Laplace (see Appendix D). Unless indicated otherwise, report the aggregated results of ten independent runs by repeating all experiments 10 times, each with a distinct DAG.
Appendix B contains additional results (e.g., other graph types, different number of nodes, and k configurations). There, CoLiDE prevails with similar observations as the cases discussed next.
Homoscedastic setting. We begin by assuming equal noise variances across all nodes. We em- ploy ER4 and SF4 graphs with 200 nodes and edge weights drawn uniformly from the range EEl 2 0.5) U |0.5, (Zheng et al. 2018 Ng al. 2020 Bello al. 2022).

---

Published as a conference paper at ICLR 2024
Classiaa
Expanesiial
Luplace
Nommatzed
COLIDEEV
GOLEM-EV
DAGMA
SunNNRegress
Notnz
Figure I: Mean DAG recovery performance, plus/minus one standard deviation, is evaluated for ER4 (top row) and SF4 (bottom row) graphs, each with 200 nodes, assuming equal noise variances. Each column corresponds to different noise distribution.
Table 1: DAG recovery results for 200-node ER4 graphs under homoscedastic Gaussian noise.
variance
variance
GOLEM
DAGMA
CoLiDE-NV
CoLiDE-EV
GOLEM
DAGMA
CoLiDE-NV
CoLiDE-EV
SHD SID
468.6 144.0
100.14118
111.9499
87.3:33.7
336.6:233..0
194,4=662
1574:44.2
105.6+131.55
22260:3951
4389=1204
5333#872
4010=11169
14472:49203
6582-1227
6067=1088
4444=1586
SHD-C FDR
473.6:2144.8
101.2:41.0
113.6:29.2
88.1=33.8
341,0=234,9
199993361
161.02143.5
107.1:151.
0.28=0.10
0.07=0033
0.08=0.02
0.06=00022
0.21=0.13
0.1540.02
0.12-0.03
0.08=0.04
IPR
0.66:20.09
0.94=0001
0.9320.01
0.95:10.01
0.76:20.18
0.92:0.01
0.93=0.01
0.95=0.01
In Figure [ we investigate the impact of noise levels varying from 0.5 to 10 for different noise distributions. CoLiDE-EV clearly outperforms its competitors, consistently reaching lower SHD. Here, it is important to highlight that GOLEM, being based on the profle log-likeidhood for the Gaussian case, intrinsically addresses noise estimation for that particular scenario. However, the more general CoLiDE formulation still exhibits superior performance even in GOLEM's specific scenario (left column). We posit that the logarithmic nonlinearity in GOLEM's data fidelity term hinders its ability to fit the data. CoLiDE's noise estimation provides more precise correction, equivalent to a square-root nonlincarity [see ( in Appendix [A.I], giving more weight to the data accurately.
fidelity term and consequently allowing fit the data
In Table I we deep-dive in two particular scenarios: (1) when the noise variance equals 1, as in prior studies (Zheng et al.) 2018 Ng al.) [2020] Bello et al. (022), and (2) when the noise level is increased to Note that SortNRegress does not produce state-of-the-art results (SHD of 652.5 and 615 in each case, respectively), which relates to the non-triviality of the problem instance. These cases show that CoLiDE's advantage over its competitors is not restricted to SHD alone, but equally extends all other relevant metrics. This behavior is accentuated when the noise variance is set to 5, as CoLiDE naturally adapts to different noise regimes without any manual tuning of its hyperparameters. Additionally, CoLiDE-EV consistently yields lower standard deviations than the alternatives across all metrics, underscoring its robustness.
CoLiDE-NV, although overparametrized for homoscedastic problems, performs remarkably well, either being on par with CoL.DE-EEV or the second-best alternative in Figure and and Table This is of particular importance as the characteristies of the noise are usually unknown in practice, favoring more general and versatile formulations.
Heteroscedastic setting. The heteroscedastic scenario, where nodes do not share the same noise variance, presents further challenges. Notably, this setting is known to be non-identifiable (Ng al.) 2020) from observational data. This issue is exacerbated as the number of nodes grows as, whether estimating the variances explicitly implicitly, the problem contains d additional unknowns, which renders its optimization harder from a practical perspective. We select the edge weights by uniformly drawing from [-1, -0.25] U [0.25, 1] and the noise variance of each node from (0.5, 10).

---

Published as a conference paper at ICLR 2024
Gaussian
Exponential
Laplace
325
380
ER4
225
Normalized SHD
180
COLIDEEV
130
COLIDE-N
GOLEM-EV
GOLEM-NV
DAGMA
SorNRs
S5
Number of nodes
Mean DAG recovery performance, plus/minus one standard deviation, under heteroscedas- tic noise for both ER4 (top row) and SF4 (bottom row) graphs with varying numbers nodes. Each column corresponds different noise distribution.
Figure 2:
Table DAG recovery results for 2o0-node ER4 graphs under heteroscedastic Gaussian noise.
GOLEM-EV
GOLEM-NV
DAGMA
Regress
CoLiDE-EV
CoLiDE-NV
SHD
642.9111.5
481.3145.8
470.21.00.6
397.8-
426.5149.6
390.7135.6
19628.11088
2569912214
1456
22560 1749
23326111200
2273411767
SHID-C
630.1153.5
509.5150066
490.2.146.8
402.1 26.3
449.4144.9
407.9 407.9137.5
FDR
0.35110.05
0.30.10.33
0.31
0.20.10.01
0.2910.04
0.25
TPR
0.33.10.09
0.60.10.07
0.64
0.62.10.02
0.68
0.68110.01
This compressed interval, compared to the one used in the previous section, has a reduced signal-to- noise ratio (SNR) (Zheng et al. 2018 Reisach et al. 2021). which obfuscates the optimization.
Figure 2 presents experiments varying noise distributions, graph types, and node numbers. CoLiDE- NV is the clear winner, outperforming the alternatives virtually all variations. Remarkably, CoLiDE-EV performs very well and often is the second-best solution, outmatching GOLEM-NV, considering that an EV formulation is clearly underspecifying the problem.
These trends are confirmed in Table 2 where both CoLiDE formulations show strong efficacy on ad- ditional metrics (we point out that SHD is often considered the more accurate among all of them) for the ER4 instance with Gaussian noise. In this particular instance, SertNRegress is competitive with CoLiDE-NV in SID, SHD-C, and FDR. Note that, Col CoLiDENNV consistently maintains lower devi- ations than DAGMA and GOLEM, underscoring its robustness. All in all, CoLiDE-NV is leading the pack in terms of SHD and performs very strongly across all other metries.
In Appendix B.S we include additional experiments where we draw the edge weights uniformly from f-om|-2, 0.5|
|0.5, 2). Here, CoLiDE-NV and CoLiDE-EV are evenly matched as the simplified problem complexity does not warrant the adoption of a more intricate formulation like CoLiDE-NV.
Beyond continuous optimization. We also tested the CoLiDE objective combination with TOPO Deng et al. (2023a), recently proposed bi-level optimization framework. Results in Ap- pendix E.8 that CoLiDE yields improvements when using this new optimization method.
Noise estimation. An accurate noise estimation is the crux of our work (Section D. leading to new formulation and algorithm for DAG learning. method's 's ability estimate noise variance reflects its proficiency in recovering accurate edge weights. Metrics such SHD and TPR prioritize detection of correct edges, irrespective whether the edge weight closely matches the actual value. For methods that do not explicitly estimate the noise, we can evaluate them a posteriori by estimating the DAG first and subsequently computing the noise variances using the residual variance formula 2² 1
dn X W T X |3 in the EV case, or, o,² || we x|3 in the NV case.
We examine whether CoLiDE surpasses other state-of-thh-art ahe approaches in noise estima- tion. To this end, we generated 10 distinct ER4 DAGs with 200 nodes, employing the homo and heteroscedastic settings described earlier in this section. Assuming Gaussian noise distribution,
8

---

Published as a conference paper at ICLR 2024
Homoscedastic noise
Heteroscedastic noise
0.3
Relative crror
CoLiDE-EV
Relative
0.2
0.45
CoLIDE-NV
DAGMA
0.40
1001000
2000
5000
5001000
2000
5000
Number samples
Namber of samples
Figure 3: Mean relative noise estimation errors, plus/minus one standard deviation, as function of the number of samples, aggregated from ten separate ER4 graphs, each comprising 200 nodes.
Table 3: DAG recovery performance on the Sachs dataset (Sachs et al. 2005).
GOLEM-EV
GOLEM
DAGMA
SonNRegress
DAGuementype
CoLiiDE
CoLiDENN
SHD
22
15
16
13
14
13
13
12
SID
49
47
56
47
SHD-C
19
11
15
13
12
11
13
14
FDR
0.83
0.66
0.5
0.61
0.57
0.5
0.54
0.53
TPR
0.11
0.11
0.05
0.29
0.17
0.23
0.29
0.35
generated different numbers of samples to assess CoLiDE's performance under limited data scenar- ios. We exclude GOLEM from the comparison due its subpar performance compared DAGMA. In the equal noise variance scenario, depicted in Figure 3 (left), CoLiDE-EV consistently outper- forms DAGMA across varying sample sizes. CoLiDE-EV is also slightly better than CoLiDE-NV, due to its better modeling of the homoscedastic case. In Figure (right), we examine the non non-equal noise variance scenario, where CoLiDE-NV demonstrates superior performance over both CoLiDE- EV and DAGMA across different sample sizes. Of particular interest is the fact that CoLiDE-NV provides a lower error even when using half as many samples as DAGMA.
Real data. To conclude our analysis, we extend our investigation to real-world example, the Sachs dataset (Sachs et al. 2005), which is used extensively throughout the probabilistic graphical models literature. The Sachs dataset encompasses cytometric measurements of protein and phospholipid constituents in human immune system cells (Sachs et al. 2005). This dataset comprises 11 nodes and 853 observation samples. The associated DAG is established through experimental methods as outlined by |aahss et al.] (2005). and it enjoys validation from the biological research community. Notably, the ground truth DAG in this dataset consists of 17 edges. The outcomes of this experiment are consolidated in Table 3 Evidently, the results illustrate that CoLiDE-NV outperforms all other state-of-the- methods, achieving SHD of 12. To our knowledge, this represents the lowest achieved SHD among continuous optimization-ben based techniques applied to the Sachs problem.
CONCLUDING SUMMARY, LIMITATIONS, AND FUTURE WORK
In this paper we introduced CoLiDE, a framework for learning linear DAGs wherein we simulta- neously estimate both the DAG structure and the exogenous noise levels. We present variants of CoLiDE estimate homoscedastic and heteroscedastic noise across nodes. Additionally, estimat- ing the noise eliminates the necessity for fine-tuning the model hyperparameters (e.g., the weight on the sparsity penalty) based on the unknown noise levels. Extensive experimental results have vali- dated CoLiDE's superior performance when compared to other state-of-the tfte-of-the-rtt art methods in diverse synthetic and real-worls settings, including the recovery of the DAG edges as well as their weights.
The scope of our DAG estimation framework is limited to observational data adhering to linear SEM. In future work, we will extend it to encompass nonlinear and interventional settings. Therein, CoLiDE's formulation, amenable to first-order optimization, will facilitate symbiosis with neural networks to parameterize SEM nonlinearities. Although CoLiDE's decomposability is demonstra- ble property, further experimental results needed to fully assert the practical value of stochastic optimization. Finally, we plan to introduce new optimization techniques to realize CoLiDE's s full potential both in batch and online settings, with envisioned impact also order-based methods.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This research was partially conducted while the first author was an intern at Intel Labs. The authors would like thanks the anonymous reviewers for their thoughtful feedback and valuable suggestions on the original version of this manuscript, which led to markedly improved revised paper.
REPRODUCIBILITY STATEMENT
We have made our code publicly available at htttps:///iith/e wtt.:h...oo..aaattt////o/ooooott ide CoLiDE-related algorithmic details including hyperparameter selection are given in Section 4.1 Additional implementation details are in Appendix
REFERENCES
Raghavendra Addanki, Shiva Kasiviswanathan, Andrew McGregor, and Cameron Musco. Efficient intervention design for causal discovery with latents. In Proc. Int. Conf. Mach. Learn., pp. 63-73. PMLR, 2020.
Albert-Láásıóó Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286 (5439):509-512, 1999.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: Learning DAGs via M-matrices and log-determinant acyclicity characterization, In Proc. Adv. Neural. Inf. Process. Syst., volume 35, pp. 8226-8239, 2022.
Alexandre Belloni, Victor Chemozhukov, and Lie Wang. Square-root lasso: pivotal recovery of sparse signals via conic programming. Biometrika, 98(4))791-806, 2011.
Peter J Bickel, Ya' acov Ritov, and Alexandre B Tsybakov. Simultaneous analysis of Lasso and Dantzig selector. Ann. Stat., 37:1705-1732, 2009.
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. In Proc. Adv. Neural. Inf. Process. Syst., volume 33, pp. 21 21865-28887,, 2020.
Peter Bühlmann, Jonas Peters, and Jan Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. Ann. Stut., 2266):252--2556 2014.
George Casella and Roger L Berger. Statistical Inference. Cengage Learning, 2021.
Bertrand Charpentier, Simon Kibler, and Stephan Günnemann. Differentiable DAG sampling. In Proc. Int. Conf. Learn. Representations, 2022.
Eunice Yuh-Jie Chen, Yujia Shen, Arthur Choi, and Adnan Darwiche. Learning bayesian networks with ancestral constraints. Proc. Adv. Neural. Inf. Process. Syst., volume 29, 2016.
David Maxwell Chickering. Learning Bayesian networks is NP-complete. Learning from Data: Artificial Intelligence and Statistics pp. 121-130, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Mach. Learn. Res., 3(Nov):507-554, 2002.
David Maxwell Chickering, David Heckerman, and Chris Meek. Large-sample learning of Bayesian networks is NP-hard. J. Mach. Learn. Res., 1:13003330. 2004.
Chris Cundy, Aditya Grover, and Stefano Ermon. BCD Nets: Scalable variational approaches for Bayesian causal discovery. In Proc. Adv. Neural. Inf. Process. Syst. volume 34, pp. 7095-7110, 2021.
Chang Deng, Kevin Bello, Bryon Aragam, and Pradeep Kumar Ravikumar. Optimizing NOTEARS objectives via topological swaps. In Proc. Int. Conf. Mach. Learn., pp. 7563-7595 PMLR, 2023a.
10

---

Published as a conference paper at ICLR 2024
Chang Deng, Kevin Bello, Pradeep Kumar Ravikumar, and Bryon Aragam. Global optimality in bivariate gradient- baasee DAG learning. In ICML 2023 Workshop: Sampling and Optimization in Discrete Space, 2023b.
AmirEmad Ghassami, Alan Yang, Negar Kiyavash, and Kun Zhang. Characterizing distribution equivalence and structure learning for cyclic and acyclic directed graphs. In Proc. Int. Conf. Mach. Lear.. pp. 3494-3504. PMLR, 2020.
Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, volume Springer, 2009.
Peter Huber. Robust Statistics. John Wiley & Sons Inc., New York, 1981.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. Proc. Int. Conf.
Learn. Representations, 2015.
Neville K Kitson, Anthony C Constantinou, Zhigao Guo, Yang Liu, and Kiattikun Chobtham. A survey of Bayesian network structure learning. Artif. Intell. Rev., 56:8721-8814 2023.
Xinguo Li, Haoming Jiang, Jarvis Haupt, Raman Arora, Han Liu, Mingyi Hong, and Tuo Zhao. On fast convergence of proximal algorithms for sqrt-lasso optimization: Don't 't worry about its nonsmooth loss function. In Uncertainty in Artificial Intelligence, pp. 49-59, PMLR, 2020.
Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclic- ity constraints. In Proc. Int. Conf. Learn. Representations, 2021.
Po-Ling Loh and Peter Bühlmann. High-dimensional learning of linear causal networks via inverse covariance estimation. J. Mach. Learn. Res., 15(():3065-3105, 2014.
J Scott Long and Laurie H Ervin. Using heteroscedasticity consistent standard errors in the linear regression model. Am. Stat., 54(3):217-222 2000.
Peter JF Lucas, Linda Van der Gaag, and Ameen Abu-Hanna. Bayesian networks in biomedicine and health-care. Artif. Intell. Med., (003))2001214 2004.
Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online learning for matrix factor- ization and sparse coding. J. Mach. Learn. Res., 11(1), 2010.
Mathurin Massias, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon. Generalized concomi- tant multi-task lasso for sparse multimodal regression. Proc. Int. Conf. Artif. Intell. Statist., pp- 998-1007. PMLR, 2018.
Eugene Ndiaye, Olivier Fercoq, Alexandre Gramfort, Vincent Leclère, and Joseph Salmon, Efficient smoothed concomitant lasso estimation for high dimensional regression. In Journal of Physics: Conference Series, volume 904, pp. 012006, 2017.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. In Proc. Adv. Neural. Inf. Process. Syst., volume 33, pp. 17943-17954, 2020.
Siqi Nie, Denis D Mauá, Cassio De Campos, and Qiang Ji. Advances in learning bayesian net- works of bounded treewidth. In Proc. Adv. Neural. Inf. Process. Syst., volume 27, 2014.
Art Owen. robust hybrid of lasso and ridge regression. Contemp. Math., 443(7):59-72, 2007.
Young Woong Park and Diego Klabjan. Bayesian network learning via topological order. J. Mach. Learn. Res., 18(1):355133482, 2017.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference: Founda- tions and Learning Algorithms. The MIT Press, 2017.
Olivier Pourret, Patrick Na, Bruce Marcot, et al. Bayesian Networks: A Oractical Guide to Appli- cations. John Wiley & Sons, 2008.
11

---

Published as a conference paper at ICLR 2024
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romen and Clark Glymour. A million vari- ables and more: The fast greedy equivalence search algorithm for learning sigh-dimensional graphical causal models, with application to functional magnetic resonance images. Int. J. Data Sci. Anal., 3:121-129, 2017.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated DAG! Causal discovery benchmarks may be easy game. In Proc. Adv. Neural. Inf. Process. Syst, volume 34, pp. 27772-27784, 2021.
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.
Andrew D Sanford and Imad A Moosa. A Bayesian network structure for operational risk modelling in structured finance operations. J. Oper. Res. Soc., 63:431-444, 2012.
Chandler Squires, Yuhao Wang, and Caroline Uhler. Permutation-based causal structure learning with unknown intervention targets. In Conf. Uncertainty Artif. Intell., pp. 1039-1048. PMLR, 2020.
Nicolas Städler, Peter Bühlmann, and Sara Van De Geer. l-ppenalization for mixture regression models. Test, 19:209-256, 2010.
Tingni Sun and Cun-Hui Zhang. Scaled sparse linear regression. Biometrika, 9944)8779-888, 2012.
Robert Tibshirani. Regression shrinkage and selection via the lasso. J. R. Stat. Soc., B: Stat. Methodol,, 58(1):267-22) -288, 1996.
Jussi Viinikaa, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. Towards scalable Bayesian learn- ing of causal DAGs. In Proc. Adv. Neural. Inf. Process. Syst., volume 33, pp. 6584-6594, 2020.
Matthew J. Vowels, Necati Cihan Camgoz, and Richard Bowden. D' like DAGs? A survey structure learning and causal discovery. ACM Computing Surveys, -3-4) 2022.
Dennis Wei, Tian Gao, and Yue Yu. DAGs with no fears: A closer look at continuous optimization for learning Bayesian networks. In Proc. Adv. Neural. Inf. Process. Syst., volume 33, pp. 3895- 3906, 2020.
Albert Xue, Jingyou Rao, Sriram Sankararaman, and Harold Pimentel. dotears: Scalable, consistent DAG estimation using observational and interventional data. arXiv preprint: arXiv:2305.1922 19215 [stat.ML), pp. 1-37, 2023.
Yang Yang. Marius Pesavento, Zhi hhi-Quan Luo, and Bjöm Ottersten. Inexact block coordinate descent algorithms for nonsmooth nonconvex optimization. JEEE Trans. Signal Process., 68:947-961, 2020.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In Proc. Int. Conf. Mach. Learn., pp. 7154-7163. PMLR, 2019.
Valentina Zantedeschi, Luca Franceschi, Jean Kaddour, Matt Kusner, and Vlad Niculae. DAG learn- ing on the permutahedron. In Proc. Int. Conf. Learn. Representations, 2023.
Bin Zhang, Chris Gaiteri, Liviu-Gabriel Bodea, Zhi Wang. Joshua McElwee, Alexei Podtelezh- nikov, Chunsheng Zhang, Tao Xie, Linh Tran, Radu Dobrin, al. Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer's disease. Cell, 153(3):707-720, 2013.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric Xing. DAGs with no tears: Con- tinuous optimization for structure learning. In Proc. Ady. Neural. Inf. Process. Syst., volume 31, 2018.
12

---

Published as a conference paper at ICLR 2024
A ADDITIONAL RELATED WORK
While the focus in the paper has been on continuous relaxation algorithms to tackle the linear DAG learning problem, numerous alternative combinatorial search approaches have been explored as well; see (Kitson et al. 2023) and (Vowels et al. 2022) for up-to-date tutorial expositions that also survey host of approaches for nonlinear SEMs. For completeness, here we augment our Section review of continuous relaxation and order-based methods to also account for discrete optimization alternatives in the literature.
Discrete optimization methods. A broad swath of approaches falls under the category of score- based methods, where (e.g., BID, BIC, BDe, MDL) scoring functions are used to guide the search for DAGs in D. e.g., (Peters et al. 2017] Chapter 7.2.2). A subset of studies within this domain introduces modifications to the original problem by incorporating additional assumptions regarding the DAG or the number of parent nodes associated with each variable (Nie et al., 2014 Chen et al., 2016 Vinikka et al. 2020). Another category of discrete optimization methods is rooted in greedy search strategies or discrete optimization techniques applied to the determination of topological
ders (Chickering] 2002 Park & Klabjan 2017). GES (Ramsey et al.] 2017) is a scalable greedy algorithm for discovering DAGs that chose one of our baselines. Constraint-based meth- ods represent another broad category within discrete optimization approaches (Bühlmann et al.) 2014). These approaches navigate D by conducting independence tests among observed variables; see e.g., (Peters et al., 2017] Chapter 7 7.2.1). Overall, many of these combinatorial search methods exhibit scalability issues, particularly when confronted with high- -igh-dimessional settings arising with large DAGs. This challenge arises because the space of possible DAGs grows at a superexponential rate with the number of nodes, e.g.. Chickering (2002.
A.1 BACKGROUND ON SMOOTHED CONCOMITANT LASSO ESTIMATORS
Consider linear regression setting y HO++ é in which we have access a response vector y Rd and a design matrix H E d comprising explanatory variables or features. To obtain a sparse vector of regression coefficients 0 Rº, one can utilize the convex lasso estimator (Tibshirani] 1996)
1 Ily H 2 + 1l
Ô argmin 0
2d
which facilitates continuous estimation and variable selection. Statistical guarantees for lasso hinge on selecting the scalar parameter > to be proportional to the noise level. In particular, under the assumption that E ~ (O, La)) solving 7 with X x ov
log
yields the minimax optimal d
solution for parameter estimation in high dimensions Bickel et al. 2009 Li et al. 2020). However, most cases having knowledge of the noise variance a luxury we may not possess.
To address the aforementioned challenge, a promising solution involves the simultaneous estima- tion of both sparse regression coefficients and the noise level. Over the years, various formulations have been proposed to tackle this problem, ranging from penalized maximum likelihood approaches
(Städler et al. 2010) to frameworks inspired by robust theory (Huber 1981 Owen 2007). Several of these methods are closely related; in particular, the concomitant lasso approach in (Owen] 2007) has been shown to be equivalent to the so-termed square-rov root lasso (Belloni et al., 2011). and rediscov-
ered the scaled lasso estimator (Sun & Zhang 2012). Notably, the smoothed concomitant lasso
(Ndiaye et al.] (017) stands out as the most recent and efficient method suitable for high-dimensional settings. The approach in (Ndiaye et al.) 2017
jointly estimate the regression coefficients 0 and the noise level a by solving the jointly conver problem
min 1 Ily 0.0
HO
2
+
llll (1
(8)
2do
2
where dg is predetermined lower bound based on either prior information or proportional to the initial noise standard deviation, e.g., do vd
xylla x 10-2 Inclusion of the constraint or 2 ag is
motivated in (Ndiaye et al.) 2017) to prevent ill-conditioning as the solution ô approaches zero To solve (8) BCD algorithm is adopted wherein one iteratively (and eyclically) solves (8) for
with fixed o,, and subsequently updates the value of a using a closed-form solution given the most up-to-date value of 0; see (Ndiaye al. 2017) for further details.
13

---

Published as a conference paper at ICLR 2024
Interestingly, disregarding the constraint 2 du and plugging the noise estimator ò Iyy H0 /2/ in (8) yields the squared-root lasso problem (Belloni al. 2011 i.e.,
min 0
ly HO + X||0||1-
(9)
This problem is convex but with the added difficulty of being non-smooth when y HO. It has been proven that the solutions to problems (8) and (9) are equivalent (Ndiaye et al. 2017). Moreover, one
can show that if ~ V(0,021, and X
log the solution of (9) is minimax optimal, and d'
notably, it is independent of the noise scale
Remark. Li et al. (0200 noticed that the non-d non-differentiabiliyy of the squared-root lasso is not an issue, the sense that subgradient can be used safely, one guaranteed avoid the singularity. For DAG estimation, due to the exogenous noise in the linear SEM, we are exactly in this situa- tion. However, we point out that this alternative makes the objective function not separable across samples, precluding stochastic optimization that could be desirable for scalability. Nonetheless, leave the exploration this alternative as future work.
A limitation of the smoothed concomitant lasso in (8) is its inability to handle multi-task settings where response data Y E dgdxn is collected from d diverse sources with varying noise levels. Statistical and algorithmic issues pertaining to this heteroscedastic scenario have been successfully addressed in (Massias et al., 20188) by generalizing the smoothed concomitant lasso formulation to perform joint estimation of the coefficient matrix @ R2"" and the square- rooot of the noise covari- ance matrix Σ = diag(nn ( .a) Rdxd, Accordingly, the so-termed generalized concomitant multi-task lasso given by
tr(E)
min e 2nd
Y - HOl 2-1
12,1 +IE En}
(10)
2d
Similar to the smoothed concomitant lasso, the matrix Eg either assumed to be known a priori, or, it is estimated from the data, e.g., via Eg
Yllr 10-22 as suggested in (Massias et al. 2018
Vnd
B ALGORITHMIC DERIVATIONS
In this section, we compute the gradients of the proposed score functions with respect to W, derive closed-form solutions for updating and Σ, and discuss the associated computational complexity.
B.1 EQUAL NOISE VARIANCE
Gradients.
We introduced the score function S(W. .)
2nd - do + A|| in we calculate the gradient of f(W. o) w.t..
1
Defining f(W,a) 1 X W W, for fixed 0 ô.
X de
The smooth terms in the score function can be rewritten
X + dā
f(W.0) Ø)
X - W
(11)
2no
+ dō
Tr
((X
W X) (X - W
x))
(12)
2nd
2
dā
Tr (XT
(I-W) WI W ) +
(13)
2nd
2
dō
= 26
Tr ((I W) cov (X) W)) + 2
(14)
where (X) 1xx is the sample covariance matrix. Accordingly, the gradient of f(W, is
= cov ccov()][ W]. Vw/(W,0) (XXT)
(XXT)L W
(15)
nă nō
(16)
14

---

Published as a conference paper at ICLR 2024
Closed-form solution for 5 in A. We update by fixing W to W and computing the minimizer in closed form. The first-order optimality condition yields
Of(W.o)
-1 2no²
d
do
X
-
0. W
X
2
2
(17)
Hence,
à² 1
X W X12
(18)
nd
1
Tr (X W (X wTX))
(19)
nd
1 nd
II X (I (I- W) - )))))
(20)
1 nd
Tr
(1I (I- (I- (I-W))
(21)
d Tr (l- cov( w))
1
(22)
Because of the constraint
2 co. the minimizer is given by
δ = max (VT-(T- (TT-W W) cov(X)(I W)) W))/d/o0)
(23)
Complexity. The gradient with respect to W involves subtracting two matrices of size d d, result- ing in a computational complexity of (dd)) Additionally, three matrix multiplications contribute to complexity of Ondd dd). The subsequent element-wise division adds (dd)) Therefore, the main computational complexity of gradient computation is O (dd), on par with sfate-of-the-art con- tinuous optimization methods for DAG learning. Similarly, the closed-form solution for updating o has computational complexity of O(d*). involving two matrix subtractions (O(2d))) four matrix multiplications (O(nd + 2d²)), and matrix trace operations (d))..
Online updates. The closed-form solution for Ô in (4) can be used in batch setting when all samples are jointly available. The proposed score function allows to design a solution for this scenario, as the updates naturally decompose in time (i.e., across samples). In the mini-batch setting, at each iteration t, we have access to randomly selected subset of data X, E dx with n, < n as the size of the mini-batch. Consequently, we could utilize mini-batch stochastic gradient descent as the optimization method within our framework. Given an initial solution Wa, it is possible to conceptualize online algorithm with the following iterations for = 1, 2.
1. Compute the sample covariance matrix cov(X,) (t 1) cov cov(X,,1) + 1 no
xXXX)
2. Compute W, using (X,) and first-order update in Algorithm 3. Compute the noise estimate G; using cov(X.) and (4).
Alternatively, we could keep sufficient statistics for the noise estimate. In this case, the update would proceed follows. We first compute the residual 6, n ||, - W and update the
sufficient statistic e = ct-1 + 4. Finally, we compute the noise estimate ā: mmx (Vie,.o))
The use of this type of sufficient statisties in online algorithms has a long-standing tradition in the adaptive signal processing and online learning literature (e.g.. Mairal et al., 2010 Although the above iterations form the conceptual sketch for an online algorithm, many important details need to be addressed to achieve full online solution and leave this work for the future. Our goal here was to illustrate that resulting DAG structure and scale estimators are decomposable across samples.
To empirically explore the mentioned sufficient statistics, conducted an experiment using mini- batches of data. The performance of these optimization methods is tested on an ER4 graph of size 1000 i.i.d samples using linear SEM. We assume the noise
d = 50, where we generate distribution is Gaussian, and the noise variances are equal. We explore two distinct batch sizes, namely 100 and 500, corresponding to using 10% and 50% of the data at each iteration, respectively. To monitor how well these algorithms track the output of the CoLiDE-EV algorithm, denoted as W
and o², we caleulate the relative error |[W- WE at each iteration t. Similarly, the relative error [
15

---

Published as a conference paper at ICLR 2024
10"
10"
Relative arror
Relarive crta
10"
10.
10a
10000
15000
15000
20000
Number of Nlrations
Number of illaaaiooss
Batch size 100
Batch size 500
Figure 4: Tracking performance of mini-batch stochastic gradient descent in relation to the output of the original CoLiDE-EV algorithm. The left plot illustrates the tracking of the output graph W while the right plot represents the tracking of the noise level o*.
for the noise scale is computed |=,-o") As outlined in Section address CoLiDE-EV for sequences of decreasing Pk- To enhance visual clarity, we focus on W, and og for the smallest μik-
la*l
The experimental results, shown in Figure 4 demonstrate that mini-batch stochastic gradient de- scent with varying batch sizes adeptly follows the output of CoLiDE-EV for both the DAG adja- cency matrix (left) and the noise level (right). While these findings showcase promising results, we acknowledge the need for further evaluation to conclusively assert the decomposability benefits of our method. This experiment serves as an initial exploration, providing valuable insights for subsequent analysis of the online updates.
B.2 NON-EQUAL NOISE VARIANCE
Gradients. We introduced the score function S(W, Σ) in 5 Upon defining
g W. Σ 2n
Tr ((X W w
Σ - w xTx) + TT))),
(24)
derive the gradient of g (W, Σ) w.r.. W, while maintaining Σ Σ fixed. To this end, note that g( (WW Σ) can be rewritten as
g( (W, Σ Tr
XX - W
X (x x W())
+ Tr(E)
(25)
2n
2
Tr
Σ
(X X - wxx) Tr(E)
(26)
2n
T (1-w)) + + Tr($)
Tr
(I - W
(27)
2
S [1 w) cov(X)(I w)) + TT(E).
Tr
(28)
Accordingly, the gradient of g( E)
Vwg(W, = 2n
T XX +XX
+X TW +
(29)
-XX
XX W
(30)
=-cov(X) [ W Σ
(31)
16

---

Published as a conference paper at ICLR 2024
Closed-form solution for Σ in (6. We update Σ by fixing W to W and computing the minimizer in closed form. The first-order optimality condition yields
- X
Vxg( W. Σ)
Wx)XX
W
Σー I = 0.
(32)
2n
Hence,
X W X) WTX)
=²=
(33)
(I- W) cov coov(X((1) W).
(34)
Note that we can ascertain Σ is diagonal matrix, based on SEM assumptions of exogenous noises being mutually independent. Because of the constraint Σ 2 Eo. the minimizer is given by
Σ = =max (Vdlag ((I W)T cov(()(( (
(35)
where v(-) and max(-) indicates an element- wise operation, while the operator diag(-) extracts the diagonal elements of a matrix.
Complexity. The gradient with respect to W entails the subtraction of two d x d matrices, re- sulting in (dd)) complexity. The inverse of diagonal matrix, involved in the calculation, also incurs a complexity of Odd). Additionally, four matrix multiplications contribute of Ond + )d The subsequent element-wise division introduces additional (df)). Conse- quently, the principal computational complexity of gradient computation Odd)), aligning with the computational complexity of state-of-the-and art methods. The closed-fome solution for updating E has a computational complexity of Odd)), involving two matrix subtractions (Od2d)), four matrix multiplications (O(nd + 2d))). and element-wise operations the diagonal elements of x d matrix (O(2d)).
complexity
Online updates. As in the homoscedastic case, CoLiDE-NY. The main difference is that we now need a collection of d sufficient statistics, one for each node in the graph. For j 1,..., and all 1,2, now have
devise online algorithm based
E,e ((x)) (W. (WTxx))))
(36)
Cj Cja-1 + EgAt
(37)
3 max (Vi (Veessoo)
(38)
where [(x:),...... (x4)) and e 0 for all j. Using Xt, the noise estimates (ō σ }
and W,-1, we can compute W., e.g., using a first-order update as in Algorithm Although these elements provide the foundation for an online algorithm, we leave the completion of this new technique as future work. any case, this shows the updates decompose across samples.
B.3 GRADIENT OF THE
ACYCLICITY FUNCTION
We adopt Husel (, s) d log(s) det W W)) as the acyelicity function in CoLiDE's formulation. As reported in Bello et al. 2022), the gradient of Hike (W, s) is given by
VHule W) 2 [sI - W W W.
(39)
The computational complexity incurred in each gradient evaluation is (dd)) This includes matrix subtraction (O)"),, four element-wise operations (O(4d)) and matrix inversion (O(d)),
Complexity summary. All in all, both CoLiDE variants, CoLiDE-EV and CoLiDE-NY, incur a per iteration computational complexity of Odd"). While CoLiDE concurrently estimates the unknown noise levels along with the DAG topology, this comes with no significant computational complexity overhead relative state-of-the-art continuous relaxation methods for DAG learning.
17

---

Published as a conference paper at ICLR 2024
C GUARANTEES FOR (HETEROSCEDASTIC) LINEAR GAUSSIAN SEMS
Consider a general linear Gaussian SEM, whose exogenous noises can have non-equal variances (i.e., the heteroscedastic case). This a non-identifiable model (Peters et al. 2017). meaning that the ground-truth DAG cannot be uniquely recovered from observational data alone. Still, we argue that the interesting theoretical analysis framework put forth in (Ghassami al. 2020 and (Ng al, 2020) as well as its conclusions carry over to CoLiDE. The upshot is that just like GOLEM, the solution of the CoLiDE-NV problem with μ4 - oo asymptotically (in n) will be a DAG equivalent to the ground truth. The precise notion of (quasi)-eqiivalencc among directed graphs was introduced by Ghassami et al. 2020), and we will not reproduce all definitions and technical details here. The interested reader referred to (Ng et al. 2020 Section 3.1).
Specifically, it follows that under the same assumptions in (Ng et al. 2020 Section 3.1), Theorems I and 2 therein hold for CoLiDE-NV when μ4 + xC Moreover, Ng et al., 2020, Corollary 1) holds as well. This corollary motivates augmenting the score function in (5) with the DAG penalty Hused [, sk), obtain DAG solution quasi-equivalent to the ground truth DAG in lieu of the termed "Triangle assumption' The proofs of these results are essentially identical to the ones in (Ng al. 2020 Appendix B), so in the sequel we only highlight the minor differences.
Let 9* and @ be the ground truth DAG and the precision matrix of the random vector x E R so that the generated distribution x ~ (0, - Let W and Σ be the weighted adjacency matrix and the diagonal matrix containing exogenous noise variances of. respectively. As - sc, the law of large numbers asserts that the sample covariance matrix cov(X) + e almost surely. Then, we drop the penalty term, 14 - oc, and considering X such that both trace terms in the score function in (5) dominate asymptotically, CoLiDE-NV's optimality optmatlly condition implies [ef. (34)
Σ (I- O-- W).
(40)
Note that when μ4 oc, CoL is convex and we are thus ensured to attain global optimality. This means that we will find a pair (W, Σ) of estimates such that (I W [I - (I-W)- e, and denote the directed graph corresponding to W by ((W). One can readily show that under the mild assumptions in (Ng et al. 2020), ((W)& is quasi equivalent to G*. The remaining steps of the
proofs of both Theorems are identical to those laid out by Ng et al. (2020 Appendix B).
D IMPLEMENTATION DETAILS
In this section, we provide a comprchensive description of the implementation details for the exper- iments conducted to evaluate and benchmark the proposed DAG learning algorithms.
D.1 COMPUTING INFRASTRUCTURE
All experiments were executed on 2-core Intel Xeon processor E5-2695v2 with a clock speed of 2.40 GHz and 32GB RAM. For models like GOLEM and DAGuerreotype that necessitate GPU processing, we utilized either the NVIDIA A100, Tesla V100, or Tesla T4 GPUs.
D.2 GRAPH MODELS
In our experiments, each simulation involves sampling a graph from two prominent random graph models:
Erdós-Rényi (ER): These random graphs have independently added edges with equal probability. The chosen probability is determined by the desired nodal degree. Since ER graphs are undirected, we randomly generate permutation vector for node ordering and orient the edges accordingly.
Scale-Free (SF): These random graphs are generated using the preferential attachment process (Barabási & Albert 1999). The number of edges preferentially attached is based on the desired nodal degree. The edges are oriented each time a new node attached, resulting in sampled directed acyelic graph (DAG).
18

---

Published as a conference paper at ICLR 2024
D.3 METRICS
We employ the following standard metrics commonly used in the context of DAG learning:
Structural Hamming distance (SHD): Quantifies the total count edge additions, dele- tions, and reversals required to transform the estimated graph into the true graph. Normal- ized SHD obtained by dividing this count by the number of nodes.
SHD-C:
Initially, we map both the estimated graph and the ground truth DAG to their corresponding Completed Partially Directed Acyclic Graphs (CPDAGs). A CPDAG repre- sents a Markov equivalence class encompassing all DAGs that encode identical conditional independencies. Subsequently, we calculate the SHD between the two CPDAGs. The nor- malized version of SHD-C is obtained by dividing the original metric by the number of nodes.
Structural Intervention Distance (SID); Counts the number of causal paths that are dis- rupted in the predicted DAG, When divided by the number of nodes, we obtain the normal- ized SID.
True Positive Rate (TPR): Measures the proportion of correctly identified edges relative to the total number of edges in the ground-truth DAG.
False Discovery Rate (FDR): Represents the ratio of incorrectly identified edges to the total number of detected edges.
For all metrics except TPR, lower value indicates better performance.
D.4 NOISE DISTRIBUTIONS
We generate data by sampling from set of linear SEMs, considering three distinct noise distribu-
tions:
Gaussian: z; ~~,00,, of), i I,...,d, where the noise variance of each node of²
Exponential: 2; Exp(A:) =1,..-. 1,...,d, where the noise variance of each node is X,
Laplace: 2 ~ Laplace (0, 6.), i L..... d, where the noise variance each node is 26
D.5
BASELINE METHODS
To assess the performance of our proposed approach, we benchmark it against several state-of-the- art methods commonly recognized as baselines this domain. All the methods are implemented using Python.
GOLEM: This likelihood- based method was introduced by Ng et al. (2020). The code for GOLEM is publicly available on GitHub at htttps:/ggitubb com. ignavier/
golem Utilizing Tensorflow, GOLEM requires GPU support. We adopt their recom- mended hyperparameters for GOLEM. For GOLEM-EV, we set \ 2 10-2 and X, 5, and for GOLEM-NY, A₁ = 2 x 10-3 and A2 = 5; refer to Appendix F of Ng et al. (020)) for details. Closely following the guidelines in Ng al. 2020), initialize GOLEM-NV with the output of GOLEM-EV.
DAGMA: We employ linear DAGMA introduced in Bello al. (2022
The code is We adhere to their rec- ommended hyperparameters: 7 4, 19. 7). 0.05, = 0.1, and (1))) 1; consult Appendix C Bello et al. (2022) for specifics.
available at /t/ttps://ggtttt thub. www.c////////// nsbe 1 / da
GES: This method utilizes fast greedy approach for learning DAGs, outlined in Ram- sey et al. (2017). The implementation leverages the py-causal Python package and can be accessed at /ttps::/github. com. cttp///diccc///yycccuu... the hyperparameters as follows: scoreId =
We configure data dataType cont ionuun and www--ccc-ssseeoo maxDegree 5,
False.
19

---

Published as a conference paper at ICLR 2024
ER4
SF4
ER4
SF4
Chusiim
Noematzed 110
Expunamta
FDR
100
Exponentta
Lapiaco
Lapbbco
Noise variance
Noise variance
EV
GOLEM-EV
DAGMA
Figure 5: DAG recovery performance assessed for ER4 and SF4 graphs with 200 nodes, assuming equal noise variances. Each row represents distinct noise distribution, and shaded area depicts the standard deviation. The first two columns display SID, while the remaining columns focus FDR.
This method adopts two-step framework involving node ordering based increasing variance and parent selection using the Least Angle Regressor (Reisach et al. 2021). The code is publicly available at Varsortability
RogrNRegress:
cotp: Scr Sttp///rrrddd/////
DAGuerreotype: This recent approach employs two-step framework to learn node or- derings through permutation matrices and edge representations, either jointly or in a bi- level fashion (Zantedeschi et al. 2023). The implementation accessible at hetps: //github. wtt..cc//vvaateddeeccc htwptesscci//AAureeetttooo DAGuerreotype and can utilize GPU process- ing. We employ the linear version of DAGuerreotype with sparseMN. as the operator for node ordering learning. Remaining parameters
defaults per their paper. For real datasets, we use bi-level optimization with the following hyperparameters: K 100, pruning.reg=0. 01, 1r.theta-0, and standards ze-True, For synthetic sim- ulations, we use joint optimization with K = 10.
E ADDITIONAL EXPERIMENTS
Here, we present supplementary experimental results that were excluded from the main body of the paper due to page limitations.
E. ADDITIONAL METRICS FOR THE HOMOSCEDASTIC EXPERIMENTS
For the homoscedastic experiments in Section for 200-node graphs at various noise levels. To offer a comprehensive analysis, we now introduce additional metrics- SID and FDR-in Figure We have also incorporated normalized SHD-C, as illustrated in Figure 6 The results in Figures 5 and 6 demonstrate the consistent superiority of CoLiDE-EV over other methods both metrics, across all noise variances.
exclusively presented the normalized SHD
E.2 SMALLER GRAPHS WITH HOMOSCEDASTIC ASSUMPTION
For the homoscedastic experiments in Section 5 we exclusively examined graphs comprising 200 nodes. Here, for more comprehensive analysis, extend investigation to encompass smaller
20

---

Published as a conference paper at ICLR 2024
Gaession
Exponenial
Laplace
Nomalize SHD
COLIDE-E
COLIDE-N
GOLEM-EY
DAGMA
SortNRegrest
Figure 6: DAG recovery performance assessed for ER4 and SF4 graphs with 200 nodes, assuming equal noise variances. Each represents distinct noise distribution, and the shaded area depicts the standard deviation.
graphs, specifically those with 100 and 50 nodes, each with nodal degree of 4. Utilizing the Linear SEM model, as detailed in Section 5 and assuming equal noise variances across nodes, we gener- ated the respective data. The resulting normalized SHD and FDR are depicted in Figure 7 where the first three rows present results for 100-node graphs, and the subsequent rows showcase findings for 50-node graphs. notable that the number of nodes decreases, GOLEM's performance im- proves, benefitting from a greater number of data samples per variable. However, even in scenarios with fewer nodes, CoLiDE-EV consistently outperforms other state-of-the-art methods, including GOLEM.
E.3 ADDITIONAL METRICS FOR THE HETEROSCEDASTIR EXPERIMENTS
For the heteroscedastic experiments in Section 5 exclusively displayed normalized SHD for various node counts and noise distributions. To ensure comprehensive view, we now include normalized SID and FDR in Figure 8 and normalized SHD-C in Figure 9 representing the same graphs analyzed in Section 5 Figures g and ( demonstrate the consistent superiority of CoLiDE- multiple metrics.
NV over other state-of-the-art methods
E.4 DATA STANDARDIZATION HETEROSCEDASTIE SETTINGS
In simulated SEMs with additive noise, a phenomenon highlighted in (Reisach al. 2021 indicates that the data generation model might inadvertently compromise the underlying structure, thereby un- expectedly influencing structure learning algorithms. Reisach al. (202T) argue that, for generically sampled additive noise models, the marginal variance tends to increase along the causal order. In response to this observation, they propose a straightforward DAG learning algorithm named Sort- NRegress (see Appendix D.5), which we have integrated into our experiments to demonstrate the challenging nature the task hand. To mitigate the information leakage caused by such inadver- tent model effects, Reisach et al. (221 suggest standardizing the data of each node have zero mean and unit variance.
We note that such standardization constitutes non-linear transformation that can markedly im- pact performance, particularly for those methods developed for linear SEMs. For our experiments here, utilized the same graphs presented in the heteroscedastic settings section of the paper but standardized the data. The recovery performance of the algorithms tests is evaluated using the nor- malized SHD and SHD-C metrics, illustrated in Figure Despite the overall degradation in per-
21

---

Published as a conference paper at ICLR 2024
ER4
SF4
ER4
SP4
Chussria
Chusiia
0.4
0.8
Expunamtt
Exponenttal
FDR
Laplace
Lapbbco
Noise variancs
Noise variance
ER4
SF4
ER4
SF4
Nomalired
FDRR
Laplace
Laplace
Noise variance
Noise variance
COLIDE-E
GOLEM-BV
DAGMA
Figure 7: DAG recovery performance for graphs with 100 nodes (the first three rows) and 50 nodes. Each graph designed with a nodal degree 4, and equal noise variances are assumed across nodes. The initial two columns display normalized SHD, while the subsequent columns present FDR. The shaded areas indicate standard deviations.
formance resulting from such standardization, CoLiDE-NV continues to outperform state-of-the-arr algorithms this challenging task.
We also conducted similar standardization on 10 small and sparse ER1 graphs, each comprising 20 nodes, with a substantial dataset of n 100, 000 i.i.d. observations. This time, we consid- ered the greedy algorithm GES (see Appendix D.5) for comparisons. Prior to standardization, the performance of the algorithms in terms of SHD was as follows: CoLiDE-EV: 12.8; CoLiDE-NV; 14.9; and GES: 14.9. Post-standardization. the SHD values changed CoLiDE-EV: 18.8; CoLiDE-
NV: 19.5; and GES: 14.9. Notably, the performance GES remained unaffected by standardization. While CoLiDE exhibited superior performance compared to others based continuous relaxation (see Figure [0). displayed some sensitivity to standardization. indicating that it is not as robust in this aspect greedy algorithms. We duly acknowledge this limitation of our approach.
22

---

Published as a conference paper at ICLR 2024
ER4
SF4
ER4
SF4
Noematzzen 122
Expununtas
125
Eapunentta
FDR
0.55
Lapplco
Lapbcco
20
Number ofnodes
Number of nodes
CoLiDE-EV
COLIDENY
GOLEM-EV
GOLEM-N
DAGMA
Figure 8: DAG recovery performance evaluated for ER4 and SF4 graphs with varying node counts under different noise distributions (each row). The figures display normalized SID (first two columns) and FDR (last two columns). We consider varying noise variances across nodes, and the shaded areas represent standard deviations.
Ganssian
Exponetial
Laplace
3.25
300
235
230
235
250
280
235
75
280
50
Noraalzed
3.0
2.99
233
235
235
20-
2.0
3.0
1.75
135
135
Number nodes
Son oss
Figure 9: DAG recovery performance evaluated for ER4 and SF4 graphs with varying node counts under different noise distributions (each row). The figures display normalized SHD-C. We consider varying noise variances nodes, and the shaded represent standard deviations.
E.5 HIGH SNR SETTING WITH HETERROSCEDASTICIT ASSUMPTION
For the heteroscedastic experiments in Section 5 we discussed how altering edge weights can impact the SNR, consequently influencing the problem's 's complexity. challenging scenario by setting edge weights to be drawn from -,,50.5] U [0.5, 5]. creating an easier counterpart to the heteroscedastic experiments in Section 5 We maintain the assumption of varying noise variances across nodes. The normalized SHD and FDR for these experiments are presented in Figure [[] As illustrated Figure [T] CoLiDE variants consistently outperform other state-of-the-art methods in this experiment. Interestingly, CoLiDE-ERY performs on par with
this section, we explore less
23

---

Published as a conference paper at ICLR 2024
ER4
SF4
BR4
Chustaa
Nortallinn liredd SED-C
Expunants
Eaponentaa
Normaltzos
Lappaco
Lamboce
Number ofnodes
Number of nodes
CoLiDE-EV
GOLEM-EV
SonNRogu
Figure 10: DAG recovery performance for ER4 and SF4 graphs when data standardization is ap- plied, considering varying noise variances for each node. These experiments were conducted across different numbers of nodes and noise distributions (each row). The first columns display nor- malized SHD, while the rest present normalized SHD-C. Shaded areas represent standard deviations.
ER4
SP4
ER4
SP4
Exponemtiall
Vormaltzed
FDR
Lapaase
83
Number of'ndds
Number of nodes
CoLiDE-EV
COLIDENV
GOLEM-EV
GOLEM-NV
DAGMA
Figure 11: DAG recovery performance in high SNR setting for ER4 and SF4 graphs, considering varying noise variances for each node. These experiments were conducted across different numbers of nodes and noise distributions (each row). The first two columns display normalized SHD, while rest present FDR. Shaded areas represent standard deviations.
CoLiDE-NV in certain cases and even outperforms CoLiDE-NV in others, despite the mismatch in noise variance assumption. We attribute this to the problem's 's complexity not warranting the need for intricate framework like CoLiDE-NV in instances.
24

---

Published as a conference paper at ICLR 2024
BR2
BRZ
Chussta
Expunamtas
Esponentta
FDR
Laplace
Lapicco
8:51
Noise variancs
Noise variance
CoLiDE-EV
GOLEM-EV
DAGMA
DAG recovery performance assessed sparse graphs with a nodal degree of 2 and 500 nodes, assuming homoscedastic noise. The first two columns depict normalized SHD, and the following columns depict FDR. Shaded regions represent standard deviations.
Figure 12:
E.6 LARGER AND SPARSER GRAPHS
is a standard practice to evaluate the recovery performance of proposed DAG learning algorithms
on large sparse graphs, demonstrated in recent studies (Ng al. 2020 Bello et al. 2022). There- fore, we conducted a similar analysis, considering both Homoscedastic and Heteroscedastic settings, with underlying graph nodal degree of 2.
Homoscedastic setting. We generated 10 distinct ER and SF DAGs, each consisting of 500 nodes and 1000 edges. The data generation process aligns with the setting described in Section [5] suming equal noise variance for each node. We replicated this analysis across various noise levels and distributions. The performance results, measured in terms of Normalized SHD and FDR, presented in Figure T2 CoLiDE-EV consistently demonstrates superior performance in most sce- narios. Notably, GOLEM-EV performs better on sparser graphs and, in ER2 with an Exponential noise distribution and high variance, it marginally outperforms CoLiDE-EV. However, CoLiDE- maintains consistent level of performance across different noise distributions and variances in homoscedastic settings.
Heteroscedastic setting. Similarly, we generated sparse DAGs ranging from 200 nodes to 500 nodes, assuming nodal degree of 2. However, this time, generated data under the heteroscedas- ticity assumption and followed the settings described in Section 5 The DAG recovery performance summarized in Figure [3] where report Normalized SHD and FDR across different noise distri- butions, graph models, and numbers of nodes. In most cases, CoLiDE-NY outperforms other state- of-the-art algorithms. In few cases where CoLiDE-NV the second-best, CoLiDE-EV emerges as the hest-performing method, showcasing the superiority of CoLiDE's variants different noise scenarios.
E.7 DENSE GRAPHS
To conclude our analysis, we extended our consideration denser DAGs with a nodal degree of 6.
Homoseedastic setting. We generated 10 distinct ER and SF DAGs, each comprising 100 and 200 nodes, assuming a nodal degree of 6. The data generation process was line with the approach detailed in Section 5 assuming equal noise variance for each node. This analysis replicated across varying noise levels and distributions. The performance results, assessed in terms of normal- ized SHD and FDR, presented in Figure 14 CoLiDE-EV consistently demonstrates superior
25

---

Published as a conference paper at ICLR 2024
BR2
SF2
BR2
22
14
022
Expunumtts
Esponentaa
Normllizaa
FDR
62
Lappaco
94
Number ofnodes
Number of nodes
CoLiDE-EV
COLIDENY
GOLEM-EV
GOLEM-NV
DAGMA
Figure 13: DAG recovery performance in large, sparse graphs with nodal degree of 2, considering varying noise variances for each node. The experiments were conducted across different numbers of nodes and noise distributions (each row). The first two columns display normalized SHD, while the rest present FDR. Shaded areas represent standard deviations.
performance across scenarios for DAGs with 100 nodes (Figure 14 bottom three rows), and most cases for DAGs with 200 nodes (Figure [4 - top three rows). Interestingly, in the case of 200-node DAGs, we observed that GOLEM-EV outperforms CoLiDE-EV in terms of normalized SHD in certain scenarios. However, closer inspection through FDR revealed that GOLEM-EV fails to detect many of the true edges, resulting in higher number of false positives. Consequently, SHD metric tends to converge to the performance of CoLiDDE-VV due to the fewer detected edges overall. This suggests that CoLiDE--EV excels in edge detection. It appears that CoL.iDE-EV encounters challenges in handling high- high-noise scenarios in denser graphs compared to sparser DAGs like ER2 or ER4, Nonetheless, it still outperforms other state-of-the-at methods
Heteroscedastic setting. this analysis, we also incorporate the heteroscedasticity assumption while generating data outlined in Section 5 Dense DAGs, ranging from 50 nodes to 200 nodes and assuming a nodal degree of 6, were created. The DAG recovery performance summarized in Figure 13 presenting normalized SHD (first two columns) and FDR across varying noise distri- butions, graph models, and node numbers. CoLiDE-N consistently exhibits superior performance compared other state-of-the-art algorithms in the majority of cases. In few instances where CoLiDE-NV is the second-best GOLEM-EV emerges as the best-performing method by small margin. These results reinforce the robustness of CoLiDE-NV across spectrum of scenarios, span- ning sparse to dense graphs.
E.8 INTEGRATING COLIDE WITH OTHER OPTIMIZATION METHODS
CoLiDE introduces novel convex score functions for the homoscedastic and heteroscedastic settings in (2) and (5). respectively. Although we propose an inexaxt BCD optimization algorithm based on continuous relaxation the NP-hard problem, our convex score functions can readily optimized using other techniques. For instance, TOPO (Deng et al. (223a) advocates a bi-level optimization algorithm. In the outer level, the algorithm performs a discrete optimization over topological orders by iteratively swapping pairs of nodes within the topological order of a DAG. At the inner level, given a topological order, the algorithm optimizes score function. When paired with convex score function, TOPO is guaranteed find a local minimum and yields solutions with lower scores. The integration of TOPO and CoLiDE combines the former's appealing convergence properties with the latter's attractive features discussed at length in the paper.
26

---

Published as a conference paper at ICLR 2024
ER6
SP6
ER6
SP6
Expunantis
Espomentta
FDR
Laplaco
Lapboco
Noise varianco
Noise variance
ER6
SF6
ER6
SF6
M--
FDR
Laplace
Lappace
Noise variance
Noise varianco
COLIDEEV
GOLEM-BV
DAGMA
Figure 14: DAG recovery performance for graphs with 200 nodes (first three rows) and 100 nodes (bottom three rows). Each graph is constructed with a nodal degree of 6, assuming equal noise vari- ances across nodes. The first two columns display normalized SHD, while the subsequent columns present FDR. Shaded areas represent standard deviations.
Using the available code for TOPO in wtlhttts://gittai com/ wttp//uuutaanoTOOOO here we pro- vide initial assessment of this novel combination, denoted as CoLiDE-TOPO. We generate 10 different ER graphs with nodes from linear Gaussian SEM, where the noise variance of each node was randomly selected from the range [0.5, 10]. For the parameters of TOPO, adhered to the values (sp 10, Ssmall = 100, and Starge 1000) recommended by Deng et al. (2023a). The results in Table show that utilizing TOPO
optimizer enables CoLiDE-NV attain better performance, with up to 30% improvements in terms of SHD and for this particular setting.
Granted, the selection of the solver depends the problem at hand. While discrete optimizers excel over smaller graphs, continuous optimization tends to be preferred for larger DAGs (with hundreds of nodes). Our experiments show that Col iiEE offers advantages in both of these settings.
27

---

Published as a conference paper at ICLR 2024
E86
SP6
ER6
SP6
Chissia
Eaponentta
FDR
Nomall
0.50
425
Lapiaco
28
16
20g
36
18
Number of nodes
Numbor of nodes
CoLiDE-EV
CO
GOLEM-EV
GOLEM-NV
DAGMA
Figure 15: DAG recovery performance in dense graphs featuring a nodal degree of 6, considering differing noise variances each node. The experiments were conducted across varying numbers of nodes and noise distributions (each row). The initial two columns illustrate normalized SHD, while the subsequent columns present FDR. Shaded areas depict standard deviations.
Table 4: DAG recovery results for 50-node ER4 graphs under heteroscedastic Gaussian noise.
CoLiDE
CoLiDE-TONe TOPO
TOPO
EV
EV
NV
SHD
100.7+29.6
100.9+23.5
128.4+38.0
88.8121.4
66.9:11544
SHD-C
102.4+13..0
102.8+23.8
131.1#37.3
91.31:0.7
68.8
FDR
0.26+t0.10
0.26+10.07
0.36=#0.88
0.23+0006
0.16#10.04
TPR
0.69+10.06
0.68+00.06
0.72=#0.06
0.7410.04
0.74:10.04
E.9 EXAMINATION
ROLE OF SCORE FUNCTIONS
Here we conduct empirical investigation the role of score functions in recovering the underly- ing DAG, by decoupling their effect from the need of a (soft or hard) DAG constraint. Given a DAG G(W), the corresponding weighted adjacency matrix can be expressed as W = II' UII, where II E {0, 1) ppdx is permutation matrix (effectively encoding the causal ordering of variables), and U Iddx is an upper-triangular weight matrix. Assuming we know the ground truth and param- eterizing the DAG as W = II UI so that an acyclicity constraint is no longer needed, the DAG learning problem boils down minimizing the score function S(W) solely with respect U
min S II UTI
(41)
Notice that 41 is convex optimization problem if the score function is convex. Subsequently, reconstruct the DAG structure using the true permutation matrix and the estimated upper- trianulaa weight matrix W ÚI.
In this setting, evaluate the CoLiDE-EV and CoLiDE-NV score functions in and In comparison, we assess the profile log- likelihod score function used by GOLEM Ng et al. 2020).
As reminder, the GOLEM-NV score function tailored for heteroscedastic settings, is given by
d
log (||x: w; xx1x1 + log (I det([I - W))) + W|]
(42)
2 i=1
28

---

Published as a conference paper at ICLR 2024
Gonssion
Exponential
Laplace
Noemalzzd canaltzed SID
Noemalized semalized
Number modees
COLiDBEV
CoLiDE-NY
GOLEM-NV
Figure 16: DAG recovery performance in ER4 DAGs under the assumption of known permutation matrices. The experiments were conducted for varying numbers of nodes and noise distributions (each column). The first row illustrates normalized SHD, while the subsequent row presents SHD- C. Shaded areas depict standard deviations.
where R" represents node i's data, and W = [wi...., ..]] Raxd The GOLEM-EV score function, simplified version of (42). given by
d log [|X- W x (.) log det W) W
(43)
2
In all cases select À 0.002 by following the guidelines in Ng al. 2020).
is important to note that the Gaussian log-likelihood score function used by GOLEM incorporates
a deieedeepn term - Σd = log of, as discussed in Ng et al. 2020 Appendix C), which is bears some similarities with our CoL-
profiled out during optimization. This noise-dependent iDE formulation. Indeed, after dropping the log (1 det W)) term in 42 that vanishes when W corresponds to a DAG, the respective score functions are given by:
CoLiDE-NV: S(W, Σ)
Tr ((X
W
X
X
W
T
X))
o X
W
2n
2 i=1
GOLEM-NV: S(W, (W
Tr (X
W X E
X W X))
Σ 1=1
log ? A/W
2n
2
While the similarities can be striking, notice that in CoLiDE the squared linear SEM residuals are scaled by the noise standard deviations o,. On the other hand, the negative Gaussian log-likelihood uses the variances of for scaling as in standard weighted LS squares. This observation notwithstand- ing, the key distinction lies in the fact that the negative Gaussian log-likelihood lacks convexity in the making it non-jointly convex. In contrast, CoLiDE is jointly convex in W and Σ due to Huber's 's concomitant estimation techniques (see Appendix A.T), effectively convexifying the Gaus- sian log-likeling ligliikel;hood; refer to Oowen (0077] Section 5) for more detailed discussion in the context of linear regression.
Going back to the experiment, we generate 10 different ER4 DAGs along with their permutation matrices. These graphs have number of nodes ranging from 10 to 50. Data generated using linear SEM under the heteroscedasticity assumption, akin to the experiments conducted in the main body of the paper. The normalized SHD and SHD-C the reconstructed DAGs are presented in Figure 16 Empirical results suggest that the profile log- ligeliieldhood score function outperforms CoLiDE in smaller graphs with 10 or 20 nodes. However, as the number of nodes increases, CoLiDE (especially the CoL.DDE LDEEEVE score function), outperforms the likelihood-hanl score functions.
29


Reviews:
[{"review_id": "73qQmYNNx7", "reviewer": "Reviewer_xkAR", "summary": "The paper studies the problem of DAG structure learning from a score-based viewpoint for linear models. \nThe authors propose a new score function that also estimates the noise levels and experimentally show that it can lead to better accuracies by leveraging recent advances in continuous non-convex characterizations of DAGs.", "questions": "* With respect to my point in the weaknesses section, in my opinion, it would be more enlightening to show that the proposed score function leads to identify the true underlying DAG. The current contribution feels like just \"another score function\" with no guarantees of identifiability. The non-equal noise variances was also studied in Loh and Buhlmann (2014) where they proposed a weighted LS that would lead to identifiability of the true DAG, this weighted LS depends on the noise levels as well, I wonder if jointly optimizing such objective would also lead to accuracy improvements.\n\n* I wonder if the authors experimented with non-linear models as well?  Given that I would consider this work to be \"empirical\", it would be good to use these ideas into nonlinear models as well. \n\n* I will also note, a recent method called TOPO by Deng et al. (2023) \"Optimizing NOTEARS objectives via topological swaps\" shows improvements in structure estimation for score-based methods. Their theory suggests that given a convex score (as in this paper) their optimization algorithm would guarantee a local optimum. It would be interesting to see if using the proposed convex score + TOPO can  obtain even more accurate DAGs, specially for non-equal variances. Finally, the same authors have provided initial insights into global optimality of continuous DAG learning methods which can also help to motivate this line of work in the continuous-constrained framework, see Deng et al (2023) \"Global Optimality in Bivariate Gradient-based DAG Learning\".", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 4, "contribution": 2, "strengths": "* The paper is clearly written and the contributions are easy to digest.\n* The proposed score leads to structure improvements w.r.t. sota methods.", "weaknesses": "* Significance: The paper considers only linear models, hindering the significance of the proposed loss function.\n* Novelty: The authors borrow ideas from concomitant lasso, and straightforwardly apply it to the score function for DAG learning. While it is totally okay with borrowing ideas from prior work, it feels that this is indeed the only technical contribution of the paper. The optimization part feels identical to prior work expect for the extra noise terms.", "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": "fhxIb4GEri", "reviewer": "Reviewer_wJd3", "summary": "The paper introduces a new continuous optimization problem for DAG learning. It leverages results from the concomitant scale estimation literature to learn a weighted adjacency matrix while estimating the scale of the exogenous noise variables. The experiments clearly show that optimizing the new objective (by inexact block coordinate descent) instead of the original l1-regularized objective of DAGMA results in better estimation of the graph across various settings.", "questions": "(minor) There is a sign typo in the second-last equation of page 14.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper tackles an important problem of interest to the general ICLR community. \n\n2. The proposed regularization is general enough that it can be plugged in many of the continuous optimization problems recently proposed for learning DAGs. The work's impact is hence potentially high as it could improve performance of many state-of-the-art methods.\n\n3. The paper is generally well presented. Its claims are well supported by an extensive empirical analysis that illustrates the DAG recovery capabilities of the method on several settings and for noise estimation.", "weaknesses": "1. Although the adjacency matrix $W$ can be efficiently updated with stochastic gradient steps, the closed-form for the noise scale is evaluated on the full data because it is not decomposable. This makes the method scale poorly to big data. This limitation should be highlighted in the text or an efficient approximation could be discussed and empirically evaluated.\n\n2. It is not clear how Problem 2 is obtained, i.e., under which assumptions the noise-dependent terms appear in the objective. It would be useful to report such derivation in the appendix. This would in particular allow for verifying if the sparsity inducing term $||W||_1$ can be replaced by the score-equivalent term $||W||_0$, used e.g. in [Brouillard et al. 2020, Zantedeschi et al. 2023], without loss in noise estimation performance.\n\n3. The paper does not describe how $\\lambda$ was tuned. In Section 4.1 it is only mentioned that it was \"empirically determined\".\n\n4. Sortnregress should be reported also in Figure 1. Currently the text reports sortnregress results for two values of the noise scale and for a single graph type, but it wouldn't hurt the readability of Figure 1 to add all the results for that baseline for all the settings. Plotting such results would allow to clearly see which settings are trivial and which are of interest.\n\nI would be inclined to increase my rating if these points are addressed.", "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": "190cR9FVC3", "reviewer": "Reviewer_VEfR", "summary": "The work proposes a new differentiable structure learning method for learning linear acyclic model that eliminates the assumption of equal error variances needed by several existing differentiable methods based on least squares. Building upon existing idea on smoothed concomitant lasso, the proposed method develops a regression-based score function that includes concomitant estimation of scale and decouples the sparsity parameter from the exogenous noise levels. Experiments with simulated and real-world datasets are provided.", "questions": "- Does the method work after data standardization (see the study by Reisach et al. (2021)? Since the method is specifically for heteroscedastic setting, this experiment should be included to support the claim.\n- For heteroscedastic Gaussian noise, the paper should compare the recovery results of Markov equivalence classes instead of DAGs, since the true DAG cannot be identified in theory.\n- Regarding performance of DAGMA and GOLEM:\n    - For DAGMA, did the authors try using the log-likelihood in the heteroscedastic setting? The authors of DAGMA paper consider such log-likelihood for nonlinear setting, but could be straightforwardly done for linear setting.\n    -  For GOLEM, did the authors use the EV version to initialize the NV version, as suggested by their paper? Also, Section 5.1 says that GOLEM is based on profile-log-likelihood--I think a more straightforward comparison with Eq. (5) is their version without profiling.\n    - Did the paper try to tune the hyperparameters for these two methods, since the settings considered here are quite different from their papers?\n- What does \"decouples the sparsity parameter from the exogenous noise levels\" mean? I did not manage to find any elaboration or explanation of it.", "limitations": "", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The problem considered is highly relevant because it is important to relax the assumption of equal error variances to handle heteroscedastic noises.", "weaknesses": "The formulation (5) in the heteroscedastic setting lacks identification guarantee. It is unclear which specific settings it is theoretically correct for. For the linear Gaussian setting, one should use Gaussian likelihood, e.g., in GOLEM, while for linear non-Gaussian setting, one should use non-Gaussian likelihood, e.g., in NOTEARS-ICA.\n\nThere are some possible issues with the experiments, elaborated in the next section.", "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)