Markov Equivalence and Consistency in Differentiable Structure Learning
Authors: Chang Deng
Kevin Bello
Pradeep Kumar Ravikumar
Bryon Aragam
Venue: NeurIPS 2024 poster (2024)
Keywords: structure learning
graphical models
causality
non-convex optimization
regularization
Primary area: causal_inference
TL;DR: By replacing commonly used score functions with a regularized log-likelihood, continuous-constrained structure learning consistently recovers the Markov equivalence class.

Abstract: Existing approaches to differentiable structure learning of directed acyclic graphs (DAGs) rely on strong identifiability assumptions in order to guarantee that global minimizers of the acyclicity-constrained optimization problem identifies the true DAG. Moreover, it has been observed empirically that the optimizer may exploit undesirable artifacts in the loss function. We explain and remedy these issues by studying the behavior of differentiable acyclicity-constrained programs under general likelihoods with multiple global minimizers. By carefully regularizing the likelihood, it is possible to identify the sparsest model in the Markov equivalence class, even in the absence of an identifiable parametrization. We first study the Gaussian case in detail, showing how proper regularization of the likelihood defines a score that identifies the sparsest model. Assuming faithfulness, it also recovers the Markov equivalence class. These results are then generalized to general models and likelihoods, where the same claims hold. These theoretical results are validated empirically, showing how this can be done using standard gradient-based optimizers (without resorting to approximations such as Gumbel-Softmax), thus paving the way for differentiable structure learning under general models and losses. Open-source code is available at \url{https://github.com/duntrain/dagrad}.

Markov Equivalence and Consistency in Differentiable Structure Learning
Chang Dengi* Kevin Bellof. Pradeep Ravikumar! Bryon Aragami
*Booth School of Business, University of Chicago, Chicago, IL 60637 Machine Learning Department, Carnegie Mellon University, Pittsburgh, PA 15213
Abstract
Existing approaches to differentiable structure learning of directed acyelic graphs (DAGss) rely on strong identifiability assumptions in order to guarantee that global minimizers of the acyclicity- constrained optimization problem identifies the true DAG. Moreover, it has been observed empirically that the optimizer may exploit undesirable artifacts in the loss function. We explain and remedy these issues programs under
by studying the behavior of differentiable eeyclicity-constrainee general likelihoods with multiple global minimizers. By carefully regularizing the likelihood, it is possible to identify the sparsest model in the Markov equivalence class, even the absence an identifiable paramctrization. We first study the Gaussian case in detail, showing how proper regularization of the likelihood defines a score that identifies the sparsest model. Assuming faithfulness, it also recovers the Markov equivalence class. These results are then generalized to general models and likelihoods, where the same claims hold. These theoretical results are validated empirically, showing how this can be done using standard gradient-based optimizers (without resorting to approximations such as Gumbel-SSoftaax thus paving the way for differentiable structure learning under general models and losses. Open- source code is available at https:// thub. wtt.hhc..ooo...ae.oooomoo..o..oc coom/duntaiin dagrad.
1 Introduction
Directed acyclic graphs (DAGs) the most common graphical representation for causal models [48, 60, 50], where nodes represent variables and directed edges represent cause-effect relationships among variables. We are interested in the problem of structure learning, i.e. learning DAGs from passively observed data, also known as causal discovery. Our focus will mainly be on score- based approaches DAG learning [11, 23), where the structure learning problem formulated as optimizing a given score or loss function s( B; X) that measures how well the graph, represented as an adjacency matrix B being acyclic. This combinatorial optimization problem is known to be NP-complete [10, 12).
{0, Fxx, fits the observed data X. constrained the graphical structure
Recent advances in score-based methods have introduced continuous representation of DAGs, transforming the combinatorial acyclicity constraint into a continuous constraint via a differentiable function that exactly characterizes DAGs [73]. In this case, the discrete adjacency matrix B (0. 10.1) is first relaxed to the space of real matrices, i.e., B E RPxP, and then a differentiable function h R + [0, oo) is devised that h(B) 0 if and only if B DAG [73, 4]. This results in the following optimization problem:
min s <(; X) subject to h(B)
BERPXP
Considering differentiable score function s, the differentiable program (1) facilitates the use of gradient-based optimization techniques along with the use richer models, such as neural networks,
"changleng@chicaggooottcc edu
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

for modeling the functional relationships among the variables [74, 67, 41, 29, 44, 28, 76]. One of the most attractive features of this approach is that it applies to general models, losses, and optimizers, in contrast to prior work. Moreover, it cleanly separates computational and statistical concerns, so that each can be studied in isolation, in the same spirit as the graphical lasso [36, 68, 17].
Looking back at the inception of the continuous DAG learning framework by Zheng et al. [73, 74], however, most developments this framework have focused on the design of alternative differentiable acyclicity functions/ / with better numerical/computatiie emphasis on which score function to use [41]. In fact, and unfortunately. regardless of the modeling assumptions, it has become a rather standard practice [67, 74, 4. 13, 32, 28] to simply use the least squares (LS) loss (a.k.a. "reconstruction loss") as the score by default, following the original paper by Zheng et al. [73], despite its known statistical limitations [63, 33, 1].
properties [4, 32, 72, 67], placing little
As a result, Reisach al. [55] flagged the empirical successes of continuous structure learning (CSL) methods as largely due to the high agreement between the order of marginal variances of the nodes and the topological order of the underlying simulated DAGs, a concept they describe as "varsortability.. Then, Reisach et al. [55] empirically showed that the performance in structure recovery of CSL methods drops significantly after simple data standardization. More recently, Ng et al. [42] demonstrated that this phenomenon may not be explained by varsortability, and instead pointed out that the explanations are due to the score function, albeit without proposing which score function to use. These observations motivate deeper consideration of the choice of score.
Unfortunately, despite the fact that several score functions have been proposed for learning Bayesian networks (such as BIC [23], BDeu [35], and MDL. [6]), their application to CSL methods is not well understood. This paper is precisely concerned with finding suitable and general score function with strong statistical properties for CSL methods. That is, our objective is to find a score function that is: (1) intrinsically differentiable so that it is amenable to gradient-based optimization without approximations; (2) applicable to general models; (3) scale-invariant; (4) capable of identifying the sparsest model under proper regularization; and (5) connects nicely with classical concepts from Bayesian networks such as faithfulness and Markov equivalence classes.
Contributions. The main contribution of our work is to show that a properly regularized, likelihood- based score function has the five properties outlined above. We begin with Gaussian models to convey the main ideas, and then discuss generalizations. In more detail:
1. (Section 4) Starting with Gaussian models, we show that using the log- log-likelihood with a quasi- MCP penalty (10) as the scoring function leads to optimal solutions of (1) that corre- spond to the sparsest DAG structure which is Markov to P(X) (Theorem 1). Furthermore, under the faithfulness assumption, all optimal solutions are the sparsest within the same Markov equivalence class (Theorem 2).
2. (Section 5) We provide general conditions on the log- likellhoooood under which similar results hold for general models (Theorem 4).
3. (Section 4.5) We show that for Gaussian models, the log- likeliidloood score is calle-ivvariant This means that rescaling standardizing the data does not change the DAG structure (Theorem 3), and hence is not susceptible to varsortability.
4. We conduct experiments to evaluate the advantages of using a likelihood-based score. The findings from these experiments are detailed in in Section 6. The empirical results support our theory: The likelihood-bart based score robust and scale invariant. We also release open-source code to facilitate the implementation and reproduction of our results.
2
Related work
Most methods for learning DAGs fall into two primary categories: Constraint-based algorithms, which depend on tests of conditional independence, and score-based algorithms, which aim to optimize a specific score or loss function. As our focus is on score-based methods, we only briefly mention classical constraint- based methods [59, 34, 62]. Within the umbrella of score-based methods, the linear Gaussian models is covered in works such as [1. 2, 19, 20, 37, 49], while studies on linear non-Gaussian SEMs are found in [33, 57]. Regarding nonlinear SEMs, significant contributions have been made in additive models [9, 15, 64), additive noise models [24, 49, 39], generalized linear
models [47, 46, 22], and broader nonlinear SEMs [38, 21].

---

Works that are more directly connected to our research include those developed in the continuous
structure learning (CSL) framework [e.g. 73, 74, 13, 4. 14, 29, 76, 41, 40, 28, 44]. Most of these papers focus on empirical and computational aspects, and only a few study the theoretical properties of the CSL framework in (1). These include: [65, 43] studied the optimization and convergence subtleties of problem (1); [13] studied optimality guarantees for more general types of score functions and proposed a bi-level optimization method to guarantee local minima; [14] designed an optimization scheme that converges to the global minimum of the least squares score in the bivariate case. Finally, among the few works that study score functiion under this framework, we note: [41] studied the properties of the & 1--regularized profile log-likeling likelihood, which leads to quasi-equivalent models to the ground-truth DAG; and the authors in [56] claim that a family of likelihood- based scores reduce to the least square loss, although this only holds under knowledge of the noise variances [33]. Perhaps most closely related to our work is [8], who proved a similar identifiability result under the likelihood score. However, they used &o regularizer along with the faithfulness assumption, which leads to an inherently non-differentiable optimization problem that is much simpler to analyze but requires approximations (e.g. Gumbel--oftmaa to optimize. On the other hand, they also consider interventional data, which we do not pursue in this work. Extending our results to include interventional data and interventional Markov equivalence is an important direction for future work. In contrast to the aforementioned works, we also prove that the log-likelihood has desirable properties such as being scale invariant, and when regularized by nonconvex and differentiable approximations of the lu function, it provably leads to useful solutions that are minimal models and Markov equivalent to the underlying structure, without assuming faithfulness.
3 Preliminaries
We let G = (V, E) denote a directed graph on nodes, with vertex set V = [p] == {l...., p) and edge set E c V x V. where (i.j) E E' indicates the presence of a directed edge from node to node j. We associate each node E V to random variable X,, and let X = Xpp).
Structural equation models (SEMs). An SEM (X, f, P(N)) over the random vector X is a collection of p structural equations of the form:
X, fz N3),, 0xf; = if l: & PAj,
where f f= (f)) is a collection of functions ff R - IR, here N (N Np) is vector of independent noises with distribution P(N), and PA, denotes the set of parents of node j. Here, Se denotes the partial derivative of Jj w.r.t. Xv. which is identically zero when ƒj is independent of XR. i.e. (X, [5(PA,,N The graphical structure induced by the SEM, assumed to be DAG, will be represented by the following p x xp weighted adjacency matrix B:
B B(), Bij - =100,f
(3)
and we use G(B) to denote the corresponding binary adjacency matrix. For any set B of SEMs, let
((B) (G(B()) (X, P(N)) E B);
(4)
i.e G(B) is the collection of all the DAGs implied by B. If is set of DAGs and B is a set of SEM, we also abuse notation by writing D B indicate D= G(B).
The SEM (2) is general enough to include many well-known models, such as linear SEMs [e.g.. 33, 491, generalized linear models [47, 45, 18]. and additive noise models [24, 51], nost-non models [70, 71] and general nonlinear SEM [38, 21, 26, 74]. To illustrate some of these models: In linear SEMS we have X, = (PA,) ) + N₁. where f, is a linear map; in causal additive models (CAM) we have X,= = ΣxEPA, f5 x) + Nj, where fz.x is a univariate function; in post-nonlinear models we have X,- = [ز, (PA,) + N2). In fact, essentially any distribution can be represented as an SCM of the form (2); see Proposition 7.1 in Peters et al. [50].
Faithfulness and sparsest representations. It is well-known that the DAG G is not always identifi- able from X, and there is a well-developed theory on what can be identified based on X under certain assumptions. This leads to the concepts of faithfulness and sparsest representations, which we briefly recall here; we refer the reader to [60, 48, 50] for details. Let I(P) denote the set of conditional independence relations implied by the distribution P. and let I[G) denote the set of d-separations implied by the graph G. Then
and faithful to G if I(P) C I(G). When both conditions hold, i.e. I(P) I[)), then G is called a perfect mup of P. Following common convention, we will simply call P faithful when ((G = I(P). When P is faithful to G,
Markov to G GG) C I(P),
the Markov equivalence class (MEC) of G is identifiable and can be represented by a CPDAG.
3

---

Definition 1. For any DAG G, the Markov equivalence class is M(G) = I IG)) I(G)
Since faithfulness may not always hold, there has been progress in understanding what can be identi- fied under weaker conditions. One approach which we will use is the notion of a sparsest (Markov) representation (SMR), introduced in [54]. A sparsest representation of P is Markovian DAG G that has strictly fewer edges than any other Markovian DAG G', and such sparsest representation is unique up to Markov equivalence class. Theorem 2.4 in [54] shows that if P is faithful to G, then G must be a sparsest representation of P. This notion is closely related to the notion of minimality we adopt in Definition 2 (cf. Lemma in the Appendix). These ideas can be generalized and weakened even further; see [31, 30] for details.
Parameters and the negative log-likelihood (NLL). For positive integers m, we will use & E C R and & E E C R to denote the model parameters for f = (f1....; fp) and N, respectively." Then we denote the distribution of X by P(X; tb, g). Let x R" denote one observation of X. P(X; ; v/, {), the negative log-likelihood and
Given n i.i.d. samples X = (XI,.... xn) where expected negative log- likelihoodd can be written as:
e, (v),)) () - Σ log P(
(((b,) -Bllog P(x: ()),
(5)
where the subscript n in fn is used to indicate the sample version of the log-likelihood.
Identifiability. Let v/S (resp. <d) denote the model parameters for the ground truth f (resp. N let B = B((b)) E RRxP denote the induced weighted adjacency matrix, and let G(B") E {0. 1)pxp denote the induced binary adjacency matrix. For example, the general linear Gaussian model (6), B represents the adjacency matrix, and & = Ω denotes the variance of the Gaussian noise. In another case, if f, is approximated by a multilayer perceptron (MLP), with Nj as Gaussian noise, then & includes all the parameters of the MLP, while & represents the variance of the Gaussian noise. Additionally, (B)₂ = [B()) = Ili-th column of Af. where Al")) the first hidden layer in ʃj [74]. Thus, by our definitions, P(X;e* (") is the true distribution. Here, there are two types of identifiability questions:
1. Parameter identifiability: Is it possible to uniquely determine the parameters (v0,5") based on observations from P(X; () { Formally, is there any (#. E) * (160, <0),, such that P(X, ( P(X, 1. E) almost surely?
2. Structural identifiability: Is it possible to uniquely determine the DAG G(B)) based on
observations from P(X; v0, { In other words, is there any () A (v0, (d) such that P(X, 00 500) P(X, v. E) but G(BD) + G(B(T)).
In general, parameter identifiability implics structural identifiability since the ability to uniquely determine parameter values often means that the structure they induce is also identifiable. However, the converse is not generally true, i.c. structural identifiability does not always imply parameter identifiability, as different parameter values can lead to the same structure. Classical results identifiability of SEMs include: linear SEM with equal variance [33], linear SEM with non non-Gaussian noises [57, 58], causal additive models with Gaussian noises [9], additive models with continuous noise [51], and post-nonlinear models [70, 71, 25].
In models where parameter identifiability is possible, the population NLL l E) serves as a natural choice for the score function because it attains a unique minimum at the true parameters (45%. However, this approach is not straightforward for nonidentifiable models, where multiple parameter sets can induce the same data distribution P(X; "i" . leading to ambiguities in parameter or structure estimation. In such cases, regularizing the log-likelihooe can alleviate this issue. These regularizers enforce specific characteristics like sparsity, guiding the model towards more meaningful solutions (e.g. faithful or sparsest), despite the lack of identifiability.
4 General linear Gaussian SEMs: A nonidentifiable model
Although our results apply to general models, we begin by outlining the main idea with one of the simplest nonidentifiable models, the Gaussian model. Our goal in this section is theoretically show
Given that s describes all the parameters for the functions f. we will also use B(b) to denote B() in (3).

---

how the NLL with nonconvex differentiable regularizers can lead to meaningful solutions such as minimal-edge edge models and elements of Markov equivalent classes. We also discuss and prove the scale invariance of NLI., making it amenable to CSL. approaches and addressing concerns raised in previous work [55, 42]. Then, in Section 5, we extend these results to general models.
4.1 Gaussian DAG models
A linear SEM (B, over X with independent Gaussian noises N. special case of (2), is well- known to nonidentifiable in terms of parameters and structure [see 2, for discussion). We write the model as follows:
X = B X + N,
where B E RRxP is matrix of coefficients with G(B) being DAG, and N E RP is the vector of independent noises with covariance matrix 52 = diag(w?,.... ..,,...,,)))
Given the model (6) it is easy to see that the distribution P(X) is Gaussian and is fully characterized by the pair (B, Q). That is:
X ~ (O( E), E =
= -) B)=
(7)
where Σ is the covariance matrix of X. In the sequel, we use the subscript f to refer to a function. In this case, y, denotes a function with arguments (B, S2) and returns the covariance matrix. Moreover, we use θ to denote the corresponding precision matrix (inverse of the covariance matrix):
+ (B = (I BT B)
(
Let X = (XI, Xn) be n i.i.d samples of X. Then, let the sample covariance matrix be
= 1 E x,X, The sample NLL function is given by:
log det Ω log det B) + 1 Ir(O((B0)) + const.
(B, - log II P(x:: B, = 2 1=1
The corresponding population NLL function is
((B Ex log P(X; B) = 2 log det Ω log de B) Tr TTO const.
The full derivation can be found in Appendix C.I. Here, it is important to note that the distribution of X is fully determined by either the precision matrix é or the covariance matrix L.
4.2 Equivalence and nonidentifiability
Our goal is to identify (B, 22): Unfortunately, the model is inherently nonidentifiable in terms of both parameter and structure. This means that multiple pairs (B, 22) for model (6) can induce the same data distribution P(X) given in (7), thus resulting also in the same precision matrix θ. To address this, we define the equivalence class (()) as the set of all pairs (B, such that (B, = e.
E(O)= = {(B O, Q) = θ}.
It is worth noting that the size of E()) is finite and most pl, which corresponds to the number of permutations for variables [2]. For more comprehensive details on this class, see Appendix C.2.
This ambiguity naturally leads to the question: which pair (B, (B,S) should we estimate? Since any pair would be indistinguishable based only on observational data, a natural objective is to estimate the "simplest" DAG, for example, DAG that induces the precision matrix O with the smallest number of edges. In other words, our goal is to estimate the matrix B that has the minimal number of nonzero entries in the equivalence class. Let S = (((i.)): Bij + 0}/-
Definition 2 (Minimality). (B, is called a minimal-edge I-map" in the equivalence class E()) sn S s Y(B, .)) E E()). The set of all minimal-edge I-maps in the equivalence class {()) is referred to as the minimal equivalence class Emin(())
Emin (0) = {B, 92) (B, is minimal edge I-map, (B, E ={())).
'In terms of notation given Section 3, we have parameters e B and & = (wf,..., w/:) and parameter
spaces V/ = E Rº00-
"This generalizes the classical definition for DAGs [e.g. 63] to refer to the entire model with the distribution and graph encoded by the matrix B and the emror variance f.
5

---

In the sequel, for brevity, we will often refer to such models as "minimal models".
Unlike faithfulness, which may always hold, the minimal equivalence class Emin (+) is always well-defined. Moreover, as detailed in Lemma 4 in the Appendix, Definition 2 is closely related to the SMR assumption [54]: Under the SMR assumption (and hence also faithfulness) for G, have M(G) min (+), i.c., Emin (9) is the Markov equivalence class of G. However, there could be multiple pairs [B, 52) within Emin (e). Nevertheless, our goal to recover one element from Emin(O) The elements in Emin (+) not only represent the "simplest" DAG model for X in terms of edge count, but also bear a deep connection to classical notions such as Markov equivalence. For example, under faithfulness, all these elements describe the same independence statements.
Lemma 1. Let X follow model (6) with (B", and @ = O ( Assume that P(X) is faithful to GD == GB.). Then M(G") = Emin )
Recall our convention that this means that M(G²) G(Emin i... the DAG structures contained in Emin(9)) coincide with M(G)). Thus, under the faithfulness assumption, recovering Emin(")) is the same as recovering the MEC, which is the usual goal in causal discovery. Morcover, emphasize that these apply generally: For non-Gaussian X following the model specified in (2), the same conclusion can be made; see Lemma 3 in the Appendix. Finally, we note that the commonly used LS loss does not have the same minimizers as the log- bog-likelihood when the noise variances are different; see Appendix C.3.
4.3 Regularization
In order to distinguish elements in (e) from the minimal elements in Emin we need to somehow account for the number of edges when evaluating the score function. The common approach to this is to use BIC, or cquivalently the to penalty. Although both approaches effectively penalize the number of nonzero entries in B, their non-differentiability makes them unsuitable for differentiable structure learning. The & penalty, while amenable to differentiable approaches, 5 is not effective in precisely counting the number of edges, and also biased in parameter estimation To mitigate these shortcomings alternatives such as the smoothly clipped absolute deviation (SCAD) penalty [16] and the minimax concave penalty (MCP) [69] have been proposed. We choose use reparametrized version of MCP, termed quasi-MCP, defined as follows:
quasi-MCP: P) = AM14
28 II( < 5) 2 1(4 > 6)]
(10)
Here, 1(.) is the indicator function; corresponding plot can be found in Appendix C.5. Similar to
MCP, quasi-MCP is a symmetric function that takes on a quadratic form between [0, 8) and remains constant for values greater than 8. The function is smooth, and for values |t] > 5, it approximates the behavior of the lo penalty, thus serving to penalize the number of non-zero coefficients in B.
The score function in (1) can be naturally written as
8 (B, S2 X, 6, X) = (B, + (A,)(B)
(11)
where PA,s (B) En px, sB,)) (B,y). Then, the optimization problem can be written
min
( 1, 8, subject to h(B) = 0, L > 0. B,52
(12)
It is worth noting that for any B, the corresponding optimal Ω that minimizes sBB X, 8, X) can be easily be expressed in terms of B as (B) (see Appendix C.I). Therefore, we can always plug 22(B) into (12) to profile out f
4.4 Provably recovering minimal models
Solving problem (12) requires minimizing é,, (B, and px.s simultaneously. To study the behavior of these minimizers, let us define the of global minimizers,
On,2.5 ((B", .*) (B*, .2") is minimizer of (12)).
(13)
5 Althoughh l is nonsmooth, standard smoothing techniques can be applied to 4 regularizers as [73, 4].
SSee Appendix C.4 for examples.

---

Ideally, we would like Nn,,55 = E (e"), however, it is unclear whether there exist values of À and ò such that any optimal solution (B*, (2*) lies within min (9). The following theorem provides an affirmative answer to this question. In the sequel, we say that property S(x:) holds for all sufficiently small x > O if there is some fixed 6 > 0 such that for every x S e, the property ((x) holds.
Theorem 1. Let X follow model (6) with from P(X), and On,,, be defined as in (13). Then, for all sufficiently small X, 8 >> 0 (independent n), it holds that P(On,,,5 = Emin (0")) + 1 as n2 + DC.
and 6 = ,, ).). Let X be i.i.d. samples
In other words, we can always guarantee that On.2.5 = Emin (+") by taking X, å sufficiently small, which is easily accomplished in practice. In the following, we use the superseript O to denote ground truth parameters. Additionally, we can assume that Bd always belongs to Emiu { ensuring that our reference to the ground truth aligns with the simplest or minimal representation within the equivalence class. Moreover, by Lemma I, under the faithfulness assumption, Theorem can be interpreted as recovering the Markov equivalence class M(G"):
Theorem 2. Consider the setup in Theorem and assume additionally that P(X) is faithiful to G = G(B"). Then, for all sufficiently small X, (independent of ni), it holds that P(On. M(G)) + 1 as
xc.
Theorem indicates with properly chosen hyperparameters, the optimal solution from optimization (12) will produce a graph that adheres to the same independence statements as GO. This implies that the structure learned through the optimization process accurately reflect the underlying causal or conditional independence structure of underlying data generating process.
Remark 1. Although we use quasi-MCP (mainly for its simplicity), it turns out MCP or SCAD can also be used. See Corollary Appendix for details.
4.5 Scale invariance and standardization
is known that the LS loss is not scale-invariant, i.e. re-scaling the data (and in particular, standard- izing it) can drastically change the structure [33]. fact which Reisach et al. [55] use argue that differentiable DAG learning with the LS Loss is also not scale-invariant. Here we show that by using different score-aie in this case the log- bog-likelihood-fii fixes this and results in (provable) scale invariance. Thus, the choice of score function is crucial if certain properties such as scale invariance are desired. The following result restates the well-known fact that Gaussian DAGs invariant re-scaling (i.e. re-scaling does not change the support for any (B, Q) E (()) using our notation:
Lemma 2. Let X (0, E), suppose E is positive definite covariance matrix and let θ = Σ-1 suppose D is a diagonal matrix with positive diagonal entries. Then G(E(O)) = G(E(DOD)).
Lemma has appealing consequences for standardization. Given raw data X. denote its standardized version by Z (cf. Appendix C.6). Ideally, structure learning algorithms will output the same structure used input, and Lemma suggests that re-scaling X will not alter the structure of the DAG that is recovered from optimizing (12). The following theorem formalizes this:
whether X or
Theorem 3. Under the same setting as Theorem I, the solutions to (12) are grale-invarian That is, for any 0, let
O, ( (B" (B", is a minimizer of (12) with data X): On,,.5 (Z)=
is minimizer of (12) with data Z),
where Z is the standardized version of X. Then, for all sufficiently small X, 8 2 0 and all n, we have g(On.x,((X) = G(O ((On,,,s(Z)) (Z)). Moreover, for all sufficiently small X, 8 > 0 we have P (X)) G(O,,s, (Z)) G(Emii (O,(B) (2)))]]
00
Thus, even on finite samples, the set of DAG structures G(On.,,s (X)) derived from the raw (un- standardized) data X will always be the same as data Z. As a result, standardizing Gaussian data does not affect the recovered DAG structure if the optimization problem (12) can be solved exactly.
which is derived from standardized
Remark 2. Theorem 3 applies to global optimization of the objective (12). Of course, in pructice, algorithms can get stuck in local optima, but the global solutions (even for finite samples n) will always be scale invariant.
7

---

Nonconvex regularized log-likelihood for general models
The results in the previous section are not specific Gaussian models, although this helps with interpretability in a familiar setting. We now extend these results from linear Gaussian SEMs to more general SEMs. Here, we assume that follows model (2) and the induced distribution denoted by ) Let us define the equivalence class 2(00.5")
EwV. = ((,,E) E): P(x; () = P(x; ,"),,, Va [RP).
That is, E(bb, < is set of pairs (v).E) that induce the same distribution P(X; .0,55) -0). As result, any pair (1/,, {) within this equivalence class will be a minimizer of the NI. ((tb, E) Analogously to Definition 2, we also define the collection of minimal elements in the equivalence class ((b)), (0)) Definition 3. (w,,), E) is called minimal-edae I-map in the equivalence class Evv",5) if iB8) S
BB V(P. 5) E E((b, ) We further define
(vb ( = (1b,)) (+),) is minimal-edge I-map, (v,E) E ((b)5)))
Here, it is crucial that our concept of minimality concerns BB() which is the number of nonzero entries in the weighted adjacency matrix ( rather than the number of nonzero entries in the parameter & itself. Therefore, (n()) essentially counts the number of edges in the adjacency matrix.
Assumption A. (1) (16", (0)) is finite. (2) B(+) is L-Liscchitz w.r.. v, i.e. 1| - -aalla
I L.
Assumption B. For any or such that E(bb <d) a, the level set {(bb,) E(bb <) a) is bounded, where f(b,E) () is the expected NLL. defined in (5).
Assumption A(1) is relatively mild; it requires that the equivalence class contains only finitely many points. This assumption is satisfied by Gaussian models, generalized linear models with continuous output [66]. binary output [74, 13], and most exponential families. It is also obviously satisfied by any identifiable model since ( _ 1. Assumption A(2) is mild continuity requirement on (()). Assumption simply guarantees that the optimization problem has a minimizer, and is standard [7]. More discussions about the assumptions are included in Appendix C.7. Without this type of assumption, score-based learning is not even well-defined.
Similar in spirit to Theorem 1, we can show that by combining the NLL with quasi-MCP for appropriate X, 6, solving the following problem, we recover elements of Emin (
min WEW.EEE
(vb. + PA.s (B(+)) subject to h M(B()) = 0,
(14)
where px,o is quasi-MCP defined in (10). Next, define its set of global minimizers.
On,,, {(v", ,²"= (v", .e") is minimizer of (14)).
Theorem 4. Let X follow model (2) with parameters (06 and let X n. i.i.d. samples from P(X; t0, )²). Under Assumptions A-B, for all sufficiently small X, 8 >> 0 (independent of n), it holds that P(On.,.5 = Emin (v6,, ,)")
as n + co.
Theorem 5. Under the setting in Theorem 4 and assuming that p" (.:) is faithful with respect to G G(B(1,)) Then, for all sufficiently small X, 5 O (independent of n), it holds that P(O = M(G)) as oc.
Experiments
To solve (12) and (14), we employ the augmented Lagrangian algorithm [5] from NOTEARS [73, 74], modifying their least squares score with & penalty into the log-likelihood with MCP (10). We compare our approach
relevant baselines, e.g. NOTEARS [73], GOLEM [41]. DAGMA [4], VarSort [55], FGES [52] and PC [59]. For our variation of NOTEARS that employs score function based on the NLL. with MCP, we name it as LOGLL-NOTEARS. The suffixes 'POPULATION' and 'SAMPLEE denote the use of the population and sample covariance matrix, respectively. Full details of the experiments are given Appendix Open-source code implementing these methods available at https: gigt.u........ com / wwwpd/unnnnnnnnee daggad...
8

---

Our primary empirical results are shown in Figures 1 and 2. We use the structural Hamming distance (SHD) as the main metric to evaluate the difference between the estimated graph and the ground truth graph. Lower SHD values indicate better estimation accuracy. Given that the model specified in (6) nonidentifiable, we compare the CPDAGs of the estimated graph and the ground truth graph.
In Figure I(a), we observe that using the NLL+MCP achieves the best performance for the different types of graphs and ranks second best for sparse graphs (ERI, JF In Figure 1(b), standardizing X significantly impacts the performance of GOLEM, NOTEARS, and DAGMA; the SHD values are not any better than an empty graph, exactly as predicted by prior theory. The performance of LOGLL-NOTEARS-S SAMPLE and LOGLL-NOTEARE LALL-NNTTARR-OPUULAAEAA POPULATION are also affected by standardization, but these methods remain robust and continue to make meaningful discoveries. It is important to note that this observation does not contradict our Lemma 2. The challenges arise because solving the optimization problems (12) and (14) to find global solutions becomes inherently difficult as p increases. To verify the scale invariance property in Theorem 3, we also conduct experiments on small graphs and include exact method that solve (12) and (14) to global optimal, see Figure 5.
In Figure 2, we replicate the Figure I in [55], providing a more direct comparison between various methods applied to raw data (X) and standardized data (X standardized). We include VarSort (referred to as sortnregress in [55]) as a baseline. Notably, for smaller graphs (p = 10), both LOGLL(- exhibit the scale-invariant property alongside PC and FGES, in aligmment with Lemma 2. This contrasts sharply with other methods, which completely deteriorate. For larger graphs (p- 50), standardizing the data mildly degrades the performance of L.G....-NOTEARSS SSMPLE FF and L-GIIL.-NOTAATT (ERTTAARSS-PPUUTT POPUUL.TTION. This can be attributed to the increased complexity of optimization as the size of the graph grows.
NOTEARS)- -OTEAAR)-SAMPLE SAMPLE and LOGLL(-NOTEARS)) PPPTTEE
In Figure 3, we use a concrete toy example to investigate two key factors in the implementation: (1) the impact of random initialization, and (2) the upper limit for do that can be applied according to Theorem 1. We generate 10s initializations Brandom with weight for each edge uniformly sampled within [-5, 5]. and perform optimization using LOGLL-NOTEARS starting from these points. The "maximal 8" is the theoretical maximum dy that ensures the validity of Theorem 1. We computed the SHD and the distances between the estimated M eaa) and M(B"). The red line in Figure represents the average SHD and distances. The distribution of these 10y estimated SHD and distances visualized using dots of varying sizes, where larger dots indicate a higher frequency of points. In some cases where SHD takes a value of 11, this value is used to indicate that the estimated Bes does not form valid DAG, which is an artifact of thresholding and affects < 0.5% of models. For the remaining models, the optimization (12) can typically be solved very close to a globally optimal, and according to Theorem 2, the SHD should ideally be zero, which is consistent with the figure.
Our results are not limited to the linear model with Gaussian noise. In Appendix E.3, we provide additional experiments on a logistic model (binary X,) and neural networks. Further details on the experimental settings and additional experiments can be found in Appendix D and E.
7 Conclusion
Continuous score- boosed structure learning is a relative newcomer to the literature on causal structure learning, which goes back several decades. It has attracted significant attention due to its simplicity and generality, however, its theoretical properties are often misunderstood. We have sought to fill in this gap by studying its statistical aspects (to complement ongoing computational studies, e.g. [41, 65, 4, 13, 14, 43)). To this end, we proposed a fully differentiable score function for structure learning. composed of log-likelihood and MCas We demonstrated that the global solution corresponds to the sparsest DAG structure that is Markov to the data distribution. Under mild assumptions, we conclude that all optimal solutions are the sparsest within the same Markov cquivalence class. Additionally, the proposed score is iaal-iinvriiant, producing the same structure regardless of the data scale under the linear Gaussian model. Experimental results validate our theory, showing that our score provides better and more robust structure recovery compared to other scores.
We hope that this work stimulates further statistical inquiry into the properties of CSL. For example, we have focused on parametric models, and left extensions to nonparametric models to future work. Certain assumptions such as the finiteness of the equivalence class and the boundedness of the level set of the likeeiilldoo become more interesting this regime. We have mentioned already that extensions to richer data types including interventions is important direction. It would be of great

---

interest to explore ways to relax our assumptions to expand our statistical understanding of CSL in broader scenarios.
150
205
030
105
0.5
180
480
185
(a) X (generated by Equation (6))
(b) standardization of (generated by Equation (6))
Figure Results terms of SHD between MECs of estimated graph and ground truth. Lower
better. Column: &o- f1, 2, 4}. Row: random graph types. {ER,SF)-k (Scale-Free, Free,Erds-Rény
graphs with kd expected edges. Here = {10, 20, 50, 70, 100), n = 1000.
PC
PC
NOTEANS
NOTEARS
Merhod
GOLEN
GOLEN
VarSert
VeeSort
Loght
LogiL sample
Loglt_papuuttton
125
(a) p 10, graph =ER"", &
p graph "ER" 2
Figure 2: Comparison of raw (orange) standardized (green) data. SHD (lower is better) between Markov equivalence classes (MEC) of recovered and ground truth graphs for ER-2 graphs with 10 (left) or 50 (right) nodes. In (b), SHD for VarSort with standardized data is omitted due to its average exceeding 300.
MEC
MATCHE PATTS
0.
3.2
0.0
0.6
3.0
masirol
madmalo
Figure 3: Graph: fork structure Xo X1 and Xo + X2. (Bea, Dest) E Emin(O") because SHD and distance are close to O.
For 0 < 8 < do. the estimated
10

---

References
[1] Aragam, B., Amini, A. and Zhou, Q. [2019], 'Globally optimal score-based learning of directed acyclic graphs in high-dimensions' Advances Neural Information Processing Systems 32.
[2] Aragam, B. and Zhou, Q. [2015], 'Concave penalized estimation of sparse Gaussian Bayesian networks', The Journal of Machine Learning Research 16(1), 1273-2328
[3] Barabási, A.-L. and Albert, R. [1999], 'Emergence of scaling in random networks', science 286(5439), 509-512.
[4] Bello, K., Aragam, B. and Ravikumar, P. [2022), 'Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization' Advances in Neural Information Processing Systems 35, 8226-8239.
[5] Bertsekas, D. P. [1997]. 'Nonlinear programming' JJuunnal of the Operational Research Society 48(3), 334-334.
[6] Bouckaert, R. R. [1993], Probabilistic network construction using the minimum description length principle, in 'European conference on symbolic and quantitative approaches to reasoning and uncertainty', Springer, pp. 41-48.
[7] Boyd, S., Boyd, S. P. and Vandenberghe, I.. [2004], Convex optimization, Cambridge university press.
[8] Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien, S. and Drouin, A. [2020]. 'Differen- tiable causal discovery from interventional data', Advances in Neural Information Processing Systems 33, 21865-21877
[9] Bühlmann, P., Peters, J. and Ernest, J. [2014], 'Cam: Causal additive models, high-dimensional order search and penalized regression', The Annals of Statistics 42(6), 2526-2556.
[10] Chickering, D. M. [1996), Learning bayesian networks is np-complete, in 'Learning from data', Springer, pp. 121-130.
[11] Chickering, D. M. [2003], 'Optimal structure identification with greedy search', JMLR 507554.
[12] Chickering, D. M., Heckerman, D. and Meek, C. [2004), 'Large-sample learning of Bayesian networks is NP-hard', Journal of Machine Learning Research 5, 1287-1330.
[13] Deng, C., Bello, K., Aragam, B. and Ravikumar, P. K. [2023]. Optimizing notears objectives via topological swaps, in 'International Conference on Machine Learning', PMLR, pp. 7563-7595.
[14] Deng, C., Bello, K., Ravikumar, P. and Aragam, B. [2023], 'Global optimality in bivariate gradient-based dag learning', Advances in Neural Information Processing Systems 36.
[15] Ernest, J., Rothenhäusler, D. and Bühlmann, [2016], 'Causal inference in partially linear structural equation models: identifiability and estimation', arXiv preprint arXiv: 1607.05980.
[16] Fan, and Li, R. [2001], 'Variable selection nonconcave penalized likelihood and its oracle properties'. Journal of the American statistical Association 96(456), 1348-1360.
[17] Friedman, J., Hastie, T. and Tibshirani, R. [2008], 'Sparse inverse covariance estimation with the graphical lasso', Biostatistics 93), 432-441.
[18] Gao, M., Ding, Y. and Aragam, B. [2020]. 'A polynomial-tim -ime algorithm for learning nonpara-
metric causal graphs', Advances in Neural Information Processing Systems 33, 11199-11611.
[19] Ghoshal, A. and Honorio, J. [2017], Learning identifiable gaussian bayesian networks in polynomial time and sample complexity, in 'Proceedings of the 31 International Conference on Neural Information Processing Systems', pp. 6460-6469.

---

[20] Ghoshal, A. and Honorio. J. [2018], Learning linear structural equation models in polynomial time and sample complexity, in "Proceddings of the Twenty-Fiig International Conference on Artificial Intelligence and Statistics', Vol. 84 of Proceedings of Machine Learning Research, PMIR, pp. 1466-1475
[21] Goudet, O., Kalainathan, D., Caillou, P., Guyon, I., Lopez-Paze D. and Sebag, M. [2018], 'Learn- ing functional causal models with generative neural networks', Explainable and interpretable models in computer vision and machine learning pp. 39-80.
[22] Gu, J., Fu, F. and Zhou, Q. [2019]. 'Penalized estimation of directed acyclic graphs from discrete data', Statistics and Computing 19(1), 161-176.
[23] Heckerman, D... Geiger, D. and Chickering, D. M. [1995]. 'Learning bayesian networks: The combination of knowledge and statistical data', Machine learning 20(3), 197-243.
[24] Hoyer, P., Janzing, D., Mooij, J. M., Peters, J. and Schölkopf, B. [2008], 'Nonlinear causal discovery with additive noise models', Advances in neurul information processing systems 21.
[25] Immer, A., Schultheiss, C., Vogt, J. E.. Schölkopf, B., Bühlmann, P. and Marx, A. [2023], On the identifiability and estimation of causal location-scale noise models, 'International Conference on Machine Learning'. PMLR, pp. 14316-14332.
[26] Kalainathan, D., Goudet, O., Guyon, I., Lopez-Paz, D. and Sebag, M. [2022], 'Structural agnostic modeling: Adversarial learning of causal graphs', Journal of Machine Learning Research 2319)) 1-62.
[27] Kingma, D. P. and Ba, J. [2014), 'Adam: A method for stochastic optimization', arXiv preprint 1rXiv:1412.6980
[28] Kyono, T., Zhang, Y. and van der Schaar, M. [2020], 'Castle: Regularization via auxiliary causal graph discovery' Advances in Neural Information Processing Systems 33, 1501-1512.
[29] Lachapelle, S., Brouillard, P., Deleu, T. and Lacoste-Julien, S. [2020], Gradient-based neural dag learning, in "International Conference on Learning Representations'
[30] Lam, W. [2023], Causal Razors and Causal Search Algorithms, PhD thesis, Carnegie Mellon University.
[31] Lam, W.-Y., Andrews, B. and Ramsey, J. [2022], Greedy relaxations of the sparsest permutation algorithm, in 'Uncertainty in Artificial Intelligence', PMLR, pp. 1052-1062.
[32] Lee, H.-C., Danieletto, M., Miotto, R., Cherng, S. T. and Dudley, J. T. [2019], Scaling structural learning with no-bears to infer causal transcriptome networks, in 'Pacific Symposium on Biocomputing 2020', World Scientific, pp. 392-402.
[33] Loh, P.-L. and Bühlmann, P. [2014], 'High- #High-dimensional learning of linear causal networks via inverse covariance estimation', The Journal of Machine Learning Research 15(1), 3065-3105.
[34] Margaritis, D. and Thrun, S. [1999], Bayesian network induction via local neighborhoods, in 'Proceedings of the 12th International Conference on Neural Information Processing Systems', pp. 505-511.
[35] Maxwell Chickering, D. and Heckerman, D. [1997], 'Efficient approximations for the marginal likelihood of bayesian networks with hidden variables', Machine learning 29(2),
[36] Meinshausen, N. and Bühlmann, P. [2006), 'High-dimensional graphs and variable selection with the lasso'.
[37] Meinshausen, N. and Bühlmann, P. [2006], 'High High--dimensional graphs and variable selection with the Lasso', The Annals of Statistics 34(3).
[38] Monti, R. P., Zhang, and Hyvärinen, A. [2020], Causal discovery with general non-linear relationships using non-linear ica, in 'Uncertainty in artificial intelligence', PMLR, pp. 186-195.
12

---

[39] Mooij. J. M., Peters, J., Janzing, D.. Zscheischler, J. and Schölkopf, B. [2016], 'Distingusshin cause from effect using observational data: methods and benchmarks' The Journal of Machine Learning Research 17(1), 1103-1204.
[40] Moraffah, R., Moraffah, B., Karami, M., Raglin, A. and Liu, H. [2020], 'Causal adversarial network for learning conditional and interventional distributions', TrXXi:2008.11377
[41] Ng. I., Ghassami, A. and Zhang, K. [2020], 'On the role of sparsity and dag constraints for learning linear dags', Advances in Neural Information Processing Systems 33, 17943-17954.
[42] Ng, I., Huang, B. and Zhang. K. [2024], Structure Learning with Continuous Optimization: Sober Look and Beyond, in "Proceedings of the Third Conference on Causal Learning and Reasoning', PMLR, pp. 71-105.
[43] Ng, I., Lachapelle, S., Ke, N. R., Lacoste-Juile Julien, S. and Zhang, K. [2022], On the convergence of continuous constrained optimization for structure learning, in 'International Conference on Artificial Intelligence and Statistics', PMLR, pp. 8176-8198.
[44] Pamfil, R., Sriwattanaworachai, N., Desai, S., Pilgerstorfer, P., Georgatzis, K., Beaumont, P. and Aragam, B. [2020], Dynotears: Structure learning from time-series data, 'International Conference Artificial Intelligence and Statistics', PMLR, pp. 1595-1605.
[45] Park, G. and Park, H. [2019a], Identifiability of generalized hypergeometric distribution (ghd) di- rected acyclic graphical models, 'The 22nd International Conference on Artificial Intelligence and Statistics', PMLR. pp. 158-166.
[46] Park, G. and Park, S. [2019b], High-ddiennioonll poisson structural equation model learning via \ell_ regression." J. Mach. Learn, Res. 20, 95-1.
[47] Park, and Raskutti, G. [2017], 'Learning quadratic variance function (qvf) dag models via overdispersion scoring (ods)., J. Mach. Learn. Res. 18, 224-1.
[48] Pearl, J. [2009], Causality, Cambridge university press.
[49] Peters, J. and Bühhlann, P. [2014], 'Identifiability of gaussian structural equation models with equal error variances', Biometrika 101(1), 219-228.
[50] Peters, J., Janzing, D. and Schölkopf. B. [2017], Elements of causal inference: foundations and learning algorithms,, MIT press.
[51] Peters, J., Mooij. J. M., Janzing, D. and Schölkopf, B. [2014), 'Causal discovery with continuous additive noise models', JMLR.
[52] Ramsey, J., Glymour, M., Sanchez-Romero, R. and Glymour, C. [2017]. 'A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application functional magnetic resonance images', International journal of data science and analytics 3, 121-129
[53] Ramsey, J., Zhang, J. and Spirtes, P. I.. [2012], "Adjaceney-faithuulness and conservative causal inference", arXiv preprint arXiv: 1206.6843.
[54] Raskutti, G. and Uhler, C. [2018], 'Learning directed acyelic graph models based on sparsest permutations'. Stat 7(1), e183.
[55] Reisach, A., Seiler, C. and Weichwald, S. [2021], 'Beware of the simulated dag! causal discovery benchmarks may be easy to game', Advances Neural Information Processing Systems 34, 27772-27784.
[56) Seng, J., Zečević, M., Dhami, D. S. and Kersting, K. [2023], Learning large dags is harder than you think: Many losses are minimal for the wrong dag, in "The Twelfth International Conference Learning Representations'.
[57] Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A. and Jordan, M. [2006), 'A linear non-gause non-gaussian acyelic model for causal discovery.', Journal of Machine Learning Research 7(10).
13

---

[58] Shimizu, S., Inazumi, T., Sogawa. Y. Hyvarinen, A., Kawahara, Y., Washio, T., Hoyer, O., Bollen, K. and Hoyer, P. [2011], 'Directlingam: A direct method for learning lin- ear non-gaussian structural equation model', Journal of Machine Learning Research-IMI.R 12()))), 1225-1248.
[59] Spirtes, P. and Glymour, C. [1991]. 'An algorithm for fast recovery of sparse causal graphs', Social science computer review 9(). 62-72.
[60] Spirtes, P., Glymour, C. Scheines, R. and Heckerman, D. [2000], Causation, prediction, and search, MIT press.
[61] Tibshirani, R. [1996], 'Regression shrinkage and selection via the lasso', Journal of the Royal Statistical Society Series B: Statistical Methodology 58(1), 267-288.
[62] Tsamardinos, I,, Aliferis, C. F., Statnikov, A. R. and Statnikov, E. [2003], Algorithms for large scale markov blanket discovery, in 'FLAIRS conference', Vol. 2, pp. 376-380.
[63] Van de Geer, S. and Bühlmann, P. [2013], 'Lo Eo-pennlized maximum likelihood for sparse directed acyclic graphs'. The Annals of Statistics 41(,), 536-567.
[64] Voorman, A., Shojaie, A. and Witten, D. [2014], 'Graph estimation with joint additive models', Biometrika 101(1), 85-101.
[65] Wei, D., Gao, T. and Yu, Y. [2020], DAGs with no fears: A closer look at continuous optimization for learning bayesian networks, in 'Advances in Neural Information Processing Systems'.
[66] Ye, Q. Amini, A. A. and Zhou, Q. [2024], 'Federated learning of generalized linear causal networks', IEEE Transactions on Pattern Analysis and Machine Intelligence.
[67] Yu, Y., Chen, J., Gao, T. and Yu, M. [2019], Dag-gnn: Dag structure learning with graph neural networks, in 'International Conference on Machine Learning'. PMLR. pp. 7154-7163.
[68] Yuan, M. and Lin, Y. [2007], 'Model selection and estimation in the gaussian graphical model', Biometrika 94(1), 19-35.
[69] Zhang, C.-H. [2010], 'Nearly unbiased variable selection under minimax concave penalty', The Annals of Statistics 38(2), 894 942.
[70] Zhang, K. and Hyvarinen, A. [2012], 'On the identifiability of the post-nonlinear causal model", arXiv preprint arXiv: 2505.2599
[71] Zhang, K., Wang, Z., Zhang, J. and Schölkopf, B. [2015], 'On estimation of functional causal models: general results and application to the post-nonlinear causal model', ACM Transactions on Intelligent Systems and Technology (TIST) 7(2), 1-22.
[72] Zhang, Z., Ng. I., Gong, D., Liu, Y., Abbasnejad, E., Gong, M., Zhang, K. and Shi, Q [2022], 'Truncated matrix power iteration for differentiable dag learning', Advances in Neural Information Processing Systems 35, 18390-18402.
[73] Zheng, X., Aragam, B., Ravikumar, P. and Xing, P. [2018], 'Dags with no tears: Continu- ous optimization for structure learning', Advances in neural information processing systems 31.
[74] Zheng, X., Dan, C., Aragam, B., Ravikumar, P. and Xing, E. [2020], Learning sparse non- parametric dags, in IIternational Conference on Artificial Intelligence and Statistics', Pmlr, pp. 3414-3425.
[75] Zheng, Y., Huang, B., Chen, W., Ramsey, J., Gong, M., Cai, R., Shimizu, S., Spirtes, P. and Zhang, K. [2024], -Causal-learn: Causal discovery in python', Journal of Machine Learning Research 25())), 1-8.
[76] Zhu, S.. Ng. I. and Chen, Z. [2020], Causal discovery with reinforcement learning, in 'Interna- tional Conference on Learning Representations'
14

---

SUPPLEMENTARY MATERIAL Markov Equivalence and Consistency in Differentiable Structure Learning
A Preliminary Technical Results
In this appendix, we include various technical results used to prove the main theorems of the paper. Proofs can be found in Appendix B.
The following corollary supports Remark 1. In the main paper, we use quasi-MCP (10) penalty in the optimization problems (12) and (14) for simplicity. However, similar conclusions hold when MCP or SCAD is used as the penalty term.
Corollary (MCPISCAD). Under the same setting as Theorem 1. Let optimal solutions collection be
On
(B*,2") is minimizer of (12) with px, replaced by p>.a PXCP(()) (t) or PA,a
(CCAD())
Then, for all sufficiently small X, > (independent of n), it holds that Oa,5a Simin(09 as + 0o, where MCP px,a (.) and SCAD p (-) are defined in Appendix C.5.
The following lemma is generalization of Lemma 1. Even for the general model, under the faithfulness assumption, all elements in the minimal equivalence class Emin (v _0)) belong to the same Markov equivalence class, as is the case in the general linear Gaussian model (6).
Lemma 3. Consider that X is generated by (2) with ()) Assume that P , v00) is faithful to G = G(B()))). Then
M(G²) (Emiin(v)) (( (G(B()) (v),) Emin( (v0,50))
where B(b) is the aljacency matrix implied by the parameterization (w) E), see (3). M(G") the Markov equivalence class of Go, see Definition 1.
Under the Sparsest Markov representation assumption, elements in the minimal equivalence class also in the same Markov equivalence class. It is important to note that the faithfulness assumption is stronger than the Sparsest Markov representation assumption. Specifically, if P is faithful with respect to G, then the pair (G, P) satisties the Sparsest Markov representation assumption.
Lemma 4. If a pair (G, P(X; B, (2) satisfies Sparsest Markov representation (SMR) (see Definition 4), then M(G) ((Eiin (e)) (G(B): (B, Emin( where M(G) is Markov equivalence class of G (see Definition! l).
The following lemma provides the formulation for the standardization of X, along with its covariance and precision matrices.
Lemma 5 (standardization). Let X ~ N(O, of = =Vr((X,) and D == diag(o1......0 Then the standardization of X, corresponding covariance matrix and precision matrix can be expressed as
D- - -EX), Cov(Xxta)= = D- [Cov(Xeta)]- = DOD
The following lemma establishes an useful identity that holds for any adjacency matrix of DAG, which is used in the derivation of the log- logllikeldoood function for the model Equation (6).
Lemma 6. If B is adjacency matrix of a DAG, then log det(I - B)
O.
The following lemma provides condition under which the optimization problem (12) is well-defined, ensuring that {BB, 52) > -o0 for any (B, ?). Lemma 7. For any (B, if Ω 0, then Σ = (B) is positive definite. Moreover, if X is generated by Equation (6) with (B", (2)), then ((B, Q)) for any (B,

---

The following lemma is used in the proof of Theorem 1. It justifies that the loss of every element in A3 is strictly greater than the loss of the ground truth, i.e. (B", ,
Lemma 8. Under the same setling and notation in the proof of Theorem I, see Section B. for any (B, 92) E E, it holds that dist B, As) > 0. there exists a > 0 such that {(B, (2) (B (2) a for all (B,1)= A3.
B Detailed Proofs
B.1 Proof of Theorem 1
Proof. It suffices to consider the population case, i... (B, is replaced by its population counter- part (B, By Lemma 7, we have
{(B -x
Also, px,o(B) 2 0 for any B. Consequently, optimization problem (12) is well-defined.
By convention, we assume that (B", 529) E Emin (0%) Now, consider the case where ox.o B3) = which equivalent Bd 0, since ((B, 92) 2 { and px,s(B) 2 P2.5 (B))( for any B. Therefore, for all À > 0 and 8 >> O, BD is the unique optimal solution to optimization problem (12), proving the conclusion.
In the subsequent proof, assume that (+") = 1. that is, Emin( {(B 52 This assumption simplifies the proof because any element of Emin (@") is indistinguishable based on the value of {BB, Q2) and the penalty for the chosen (6, A), as shown below. Our goal is to identify one element via optimization problem (12), which significantly simplifies the argument.
First, let us define
60 I+
where
= min min (B))) ((
[iy min
min ("
[ (T))
with any A 0. (a) is due to the fact that each element in equivalence class E(Oº) is one-to-ons associated with (T), see Section C.2 or [2] for detailed discussion. Then, for any À 0 and
0 < 8 < da. consider the set A {(B,5)) PP.a(B) PB))(B))) For any (B, 52) E A1, we have (B, & E(O") \ ((B", ,,)).. This follows from the fact that
Px. (B) = A6 SB
Ad
Sga = PA.s(Ba)
V(B E E(0)) ((B922))
2
2
As consequence, this implies that {(B", < ((B (( E A1. Therefore,
((B", PA.5 )B)) (B, + (B)
Y(B A
Next, we define A2 = ((B.2) sA.s (B) > px.s (B))). Since (BB,, 2 < ((B, (2), follows that for all (B, 22) E A2, the following inequality holds:
((B", (2) Px.s ) ((, + (A)
Y(B A2-
Therefore, we need examine the set A3 {(B, (B) < px,s (X,)(B)). For (B" to achieve the minimum value of the score function, it is crucial that the following condition is satisfied:
((B", Ω2) + (B*) < (BB Ω) + (B)
V(BΩ)) E Ag.
This condition guarantees that the ground truth parameters (D", Ω2) correspond to the optimal solution by comparing their score with any other parameters in the subset Ag.
It is important to note that (t) = Xp1,5(t), Vt. Thus, necessary and sufficient condition for this hold is:
X < min (B, - ((B",22)
(B, P1 (B)) P1.5 (B)
16

---

2 Su, with equality achieved when B O. Therefore, the denominator on the RHS cannot be arbitrarily large. Moreover, since (B, E A3. it follows that (B,S)) & E(("), as A3 n E((") 0.
Note that for all (B, A3, we have PI.) P1,s(B) S 6
We define the distance from B to the Ag
dist( B, A3)
inf (B,0)EAs
|
B
For all (B, E( it turns out that dist( B, As) must be positive due to the design of óo, giving:
di Aa)
min (6))
min
Bisl - 60
((..)]...,
A
1+A 1+ > 0.
By Lemma 8, there exists some a > such that ((B {(Bo, > a for all (B, E Ag. Consequently, have:
(B,0)) - ((BD, (B PI,5( - (B) 0.
inf
Thus, we can define
(B,0)) ,(")
A inf
(B Pi (8) P1 >
In summary, for all 0 A A and 0 5 < 6n. for any (B, Emiu (89), and for all
(B, ()) & Emin (0%), the following holds:
{(B PA.s (B) ((B PB).
This concludes the proof.
B.2
Proof of Theorem 2
Proof. Here, O, 22) = O". From Theorem 1, we know that when n
On, = Emin (0*).
Given the additional assumption that p(X) is faithful with respect G == G(B)), by Lemma 1, have
M(G)) (09) (G(B) (B, E Emin (0")
Note that
(G(B) : (B,S2 Emin (09) = (G(B): (B,9) E On,.5). Emin
Thus, conclude that
M(G²) = On,266
as m2 1
This completes the proof.
B.3 Proof of Theorem 3
Proof. In this proof, we use the notation introduced in Section C.6. Note that when X is used in (12),
we essentially compute the sample covariance matrix E based on X as follows:
-= 1 x (f.....,,pp)) 7b)) X - I, (T...... (P......74 Pp)
n
and plug it into the negative sample logelikeldoood function. The same procedure applies to Z:
Eyy
D X 1, .......... Rp)i1 X (p,..... ...h................. .P)) D-1
-
17

---

Denote ② = (E)-1 and Osu (Ema)- = DÔD. For θ by applying Theorem 1, there exist 10 Ap
and 00 > such that for any 0
1 Apwrn and 0 8 % di we have
OnA.s () = Emin (0). For DÔD, we apply Theorem again, and there exist Ao Aud. > 0 and
50d.n O such that for any 0 < A < Age. and 0 8 < dgd.n, we have O,,,,6 Emin (DOD).
>oul. } in optimization (12) ensure We select on min ( i 00 Soul." and An min (1) (Apm. Ao
that O,,,,5(X) (X) = Emin ()) and O,,x,(Z) = Emin (DOD) hold simultancously.
By applying Lemma 2, we conclude that:
g(Om,A, G(Emin (e)) ( (DGD) ((OnnA,,
Furthermore, as n + 20, we have Σ + E and θ + O. Therefore, Emiu (e) + Emiu (8), and Emin (Qsu) + Emin ( (Qau) as + oc. Thus,
where n - oo Note that we use - to indicate we consider the result in population level.
G(O,,,, (X)) = (mmin (0)) = g(On,
B.4 Proof of Theorem 4
Proof. This proof shares many similarities with the proof of Theorem 1. First, when n -> 00, we consider the result at the population level. Thus, (v). g) ((v,, E), and we will focus E(tb {) in the following. As a result, we only work with E((), instead of (xb, E)-
By convention, we can always assume that (w/ ) E (min(
Now, consider the case where Px,s p P(X, 60, (,,,)), we have
which implies that B(b))) 0. Since X
f(vb E((,E),
and thus
1(60, + PA. S
A(b E V E, YA 0, 0.
Therefore, (w0,55) is the optimal solution to optimization (12).
As we iterate, we can assume that ( assumption simplifies the proof because any element in Emin (v/0 5 c0) is indistinguishable based on the value of f {) and the penalty for the chosen parameters (6, A). Our goal is to find one element by solving the optimization problem (14), and this assumption simplifies the argument significantly.
1, meaning E Emin (v (V),<") This
First, we define
80
min [Bb)) (W. ,0) B(O), #0}
min
A
is important to note that, under Assumption A (1), since 18(,0,, {0) is finite, we have > 0.
Then, for any À >> and 0 < 5 < dg, consider the set
A1 = { PA.6 (B())) = Px, (B(b5))
It clear that for any (w.E) E A,, we must have (v),E) & E(W,, , since for any (v,E) E E(60.
we have sx,s( B (v+) Px, PB,((B( (B(w))). As result,
(B.
((b0,59) (()))) v(v,E) E A
Therefore,
PB.s(()))))) < ((v,E) + PPo(((())) v(W,E) A. Next, we consider the set A2 ((v,E) px,(((()) > (B)))()))),, and we know that E1150 ( s ((),E). Therefore,
((v0.c") PA,s (B(x)) < ((},E) + PA.((B())
v(v,E) E A2.
18

---

Consequently, we need to check As = {(tb.E) (B,)( < PA. (B,) For (w0.E) to be the minimizer of (14), we require that the following condition holds:
l(v 5 PA, (B((
((),E) PA (B(())
((V, A3.
This condition ensures that the ground truth parameters (5)) correspond to the optimal solution by comparing their score with that of any other parameters in the subset A3.
It is also worth noting that px,s(t) (t), for all t. Therefore, necessary and sufficient condition for this to hold is:
((1,E) (0)
inf
(ViE)EAs P1, (B,(()))))) (1,)(
Note that for (w) E) E A3, we have (B( (B,)( S 8 ( Therefore, the de- nominator the RHS cannot be arbitrarily large. Moreover, for any 0 6o, the following holds:
A & B((") B( & LIb 11l2
V(b,E) E A3.
1
The second inequality follows from Assumption A (b). As a consequence,
11 E) (eb", 1l) 2 1 2112 z
L(1 + A)
Thus, we obtain
A
Ag S ((b,E) (1(,E)) - (30 2
= A4.
L(1 A)
First, note that A3 a nonempty set. Otherwise, the conclusion would hold immediately. Let us select any E) E As + 0, and define As, = {(v+E) ((bb <) - {) Then,
E(w), ()),)) 2 inf inf inf E((bE)) (w.D)A,nAs (+tE)EAs
((vbE).
(A.
It is important note that A, As is nonempty, since (15,) E A₂ Ag- By Assumption B and the properties of f(b, E). we know that As is a bounded and closed set, and A₁ is closed set. Consequently, A4sAs compact. Furthermore, for all (v)&) E E(, _²), we have (v, E) & A,S All of this leads to the following conclusion:
inf E(bbE) Z
min (nAA)
E(WiE) inf (a
((b,E) > E(w.
(w.E)
As a result, we define An as follows:
inf ((vbE) 0(160, ( (wE))A PI,(((()))) - P1,s (I,s(((b))
Ao
>0.
B.5 Proof of Theorem 5
Proof.
The proof is combination of Lemma 3 and Theorem 4, similar to Proof of Theorem 2.
B.6 Proof of Lemma
Before proving the result, we introduce the definitions of the Sparsest Markov representation assump- tion and restricted faithfulness, along with few useful theorems.
Definition 4 (Sparsest Markov representation [54]). A pair (G", P) satisfies the Sparsest Markov Representation (SMR) assumption if (G P)) satisfies the Markov property and |G] |G for every DAG G such that (G, P) satisfies the Markov property and G & M(G)).
In other words, the SMR assumption asserts that the true DAG G is the (unique up to Markov equivalence) sparsest DAG satisfying the Markov property.
19

---

Definition 5
[53, 54D. A distribution P satisfies the restricted faaricteddfathfuulnees assumption with respect DAG if is Markov G and following two conditions hold:
Adjucency-faithuuues:: for all ().)) E and all subsets S C [p] it holds that Xj IL Xx Xs
for all triples k,, with skeleton - k and all subsets S C. [pl' k) such that j is d-connected to ls given S it holds that X; J. X2 Xs
Orientation-faittffllee
Theorem 6 (153)). If distribution P is faithful to G, then such distribution also satisfies the instrcced-faittfulhee assumption with respect to G.
Theorem 7 (Theorem 2.4 in [54]). Let (G, P) satisfy the Markov property. Then the restricted- faithfulness assumption implies the SMR assumption.
Proof. First, by Theorem 6, the faithfulness assumption implies the restricted faithfulness assump- tion. Second, by Theorem 7, the restricted faithfulness assumption implies the Sparsest Markov representation assumption. Furthermore, note that for any (B, 92) E G(Emin (@")). the distribution P(X) is Markov to G (B)), since (B, .2) E ((E()))). According
the definition of the Sparsest
Markov representation assumption, all sparsest DAGs that satisfy the Markov property must belong the same Markov equivalence class. In our case, this means ((miin M(G').
B.7 Proof of Lemma 2
Proof. For ~ ,_,N0,,)) Σ). let θ = DeD, and denote the inverse of as E = - It follows that X = D X (0, Σ). Now, consider the following least squares regression for j E (I,...,,p and 8 C {1..... p} (i). Let 3 e gl Then the following relationships hold:
Bsj = arg min B||x, 3 x 2
- 8sj Egg
Bs, = arg min B x 3 Xs Beg EigKis
(ss) Dss Dss
Egj -DD-' 1sj = D5 Djj
Σsj = -5g = Dss Deg D3 Ls; DJj Dss
As a consequence, supp(3ss) supp(3s;) Note that for all (B, E E()), we know from Section such that B = B(x). Moreover, B can be recovered by least squares
C.2 that there exists a regression using X with its topological sort [2, 13] that is consistent with æ. For such a æ, we can find a pair (B, E (DDD), where B has the same topological sort as n, and it can be recovered by least squares regression on X. We have shown that, for the same S, supp(3ss) = supp(3s:). Therefore, supp(B) (B) = supp(B).
B.8 Proof of Lemma
Proof.
The proof is the same as Lemma 1.
B.9
Proof Lemma 4
Proof. This follows directly from the definition of the Sparsest Markov Representation (SMR) assumption. Since for all (B, E (O), G(B) is Markovian to P and G(B) is the Sparsest, by Definition 2, all G(B) must belong the same Markov equivalence class by the definition of SMR.
B.10 Proof of Lemma 5
Proof. Xsu D EX) is based on definition of standardization.
= Cov(Xaa) = Cov( Ex))= D= - E-))
D-1ED-1
(Cov(X = - - DOD.
20

---

B.11
Proof of Lemma 6
Proof. Detailed proof can be found in [41], Appendix Section
B.12
Proof of Lemma 7
Proof. From the definition of (B,
Ey = B)- -B) - -B)-1 (T - B)
where 12//2 diag(w)...... wp). It is clear that E( is positive semidefinite, as E(B, ||2|/2 -) 21 =1:20 0,
E
Next, we just need show that /2- B3
for all x B) + + B) B
11/2 B)
Here, w > 0 for all j. so //2 is invertible. As (1 - B) is full rank matrix, then (1 B) is also a full rank matrix, it indicates that Σ(B, [2) is positive definite matrix.
Since Ω2 0, follows that Σd is positive definite. By Lemma have:
8(B,2)
2 log det - log det(1 - B) + 2 TTE" + const.
2 log det 1 Tr(E"((B,1)) const.
1
2
BB (B)) ((B)
where (B) and ((B) are defined in Equations (15) and (17), respectively, The last inequality follows from Section C.I. Next, we need to prove that ((B) -00.
((B) 2 log det diag((1 (I- B) ET B)) const.
log E |X, - -BJ XI2 const.
2
P
Σ
log E Ei((e, - B₁) x + const.
2
i=1
Σ
log(e; -Bi) T2 (0) B₁) + const.
2
j=1
1
log [I(e; BBy) (2º) const.
2 j=1
1
P
Σ
log A Amii (En)
2
Here, E RP a unit vector with the j-th position equal to and all other positions being zero, and
Amiul (Σº) ) is the minimum eigenvalue of Σº Since n is positive definite, we have Amiu (Σo) 0.
Because B is the adjacency matrix of DAG, it follows that B,N = O, which implies ||ez By 2
||e,|| 1. As result,
8 BB, Ω2) 2 ((B) -O.
B.13 Proof of Lemma 8
Proof. Note that for a fixed B, the corresponding optimal (B) = diag(( TT B)) is the solution with respect to ((, Therefore, without causing confusion, take 22 (B)
21

---

diag((I - T B)) and consider the log-likeliban leellike as function of B only, i.e., {BB), for simpler representation. See Equation (17) in Section C.I for details. It is clear that Ag, so define A4 = {B (BB) z {(0)). Note that ((0) is finite.
((0) 2 (B)
log(e, - [ B3)
2
Σ j=1
log ||e, - Amin
p II
log
(Emin(""))²
lle Byll
==1
This indicates that
II lej Byll2 5 exp(e(0))
(Amiu (Eo)
3=1
Moreover,
le, - B,1 < exp(((0)) A
lee Bxll S
II 1=1
Vk E {l.....p}
This implies that BN must be bounded, and therefore every B in A4 is bounded. It is clear that arg mill E(B) A4. Thus, we need to show that minne ((B) minneegnA, ((B) > ((B)). Define
A5 = (B dist B, (,,(()))) 2
min dist( B, As) 2 BEE(EM)
is easy to see that Ag Ag. Then.
min BEAsNAa
{BB)
min
{(B).
BEAdNAs
Note that A4 n A5 is closed, bounded, and nonempty (0 E A4 n As), and L(B) is continuous function of B. Consequently, there exists at least one minimizer ((B) in Aq A.Ag Combining this with the fact that B A4 A, for all B E(O"), conclude that:
min {(B) 2 BEAs
min BEAdnAs
((B) > E(B))
B.14 Proof of Corollary
Proof. By Theorem 1, we know there exists An > and 6n > 0. For MCP, it can be transformed into quasi-MCP, by reparameterization from Section C.5. Then, combining these results together.
80
Amop Ag,
admep < 60
Amep < do. <a<<<<
Amcp
5u and (Amop)o = Ao. For SCAD, we just requires the following is satisfied We could simple set ag Xu to satisfies the pattern the proof Theorem 1.
12 (a 1) Au
0 < adsend < So, 0< <
One simple choice is to let
ag = (Dueado < minff Von, Voo, 1}
This completes the proof of Corollary 1.
22

---

C Additional Examples and Details
In this appendix, we provide the additional details of derivations, examples, concepts, and discussions referenced in the main paper. These include:
The derivation of the log- fog-likelihood function for the model in Equation (6) (Appendix C.I).
A brief introduction to the characterization of the equivalence class {()) (Appendix C.2).
Examples demonstrating that the optimal solution for the least squares loss differs from the optimal solution of the logllikelihood (Appendix C.3).
An example illustrating the estimation bias when the & penalty is applied (Appendix C.4).
The formulations for quasi-MCP, MCP, and SCAD (Appendix C.5).
The standardization the random variable and the dataset X (Appendix C.6).
A detailed discussion of Assumptions and B (Appendix C.7).
C.1 Log-likelihood of Model (6)
In this subsection, we detail the negative log- bog-likelihood of the model in Equation (6).
e, (B, = 71 log (x:: B,
x ( 2)x;
log
(2 227 E(B,2)) exp
71
P log 2T + logdet det E(B
1
x @(B
2
2n
1 log det (I -T B) + Tr()(B,
E x + const.
2
1
log det log det(I - B) Tr(Se(B,0)) const.
2
2
1
P
log w/
Tr(O (I- B)T TĀ[I B)) log det(I - B) const.
2
f=1
1
P
p ((I- - T) B)); log det(I - B) + const.
Σ
log w
2
2 j=1
w/
j=1
It is easy to know that for any fix B, the optimal solution of (wz)² can be written as:
(w:)² = [(T -BI-BB)I---B B)];
Therefore, optimal solution (B) for any fixed B can be written as:
(B) = diag((I - B) B) E(I B))
(15)
Let us define profile sample log- likelihooho (B) as function of B with such optimal 2(B) plugged in
1
(B) log det diag((I -B) T B)) log det(I B) const.
2
(16)
1
log - |X XB, log det B) + const.
2
Where X (X1,....
X,) (x1; (x)...., xn) and corresponding profile population log-likelihood
{(B)
log det diag((T (TT B) T - B)) - log det((- - B) + const.
2
(17)
log E BJ XI| log det([ I B) + const.
2
23

---

C.2 Equivalence class E())
We provide brief introduction to the equivalence class {()), which has been extensively studied in [2]. We adopt the notation from [2], and further details can be found in that work.
Definition 6 (topological sort). a topological sort of a directed graph an ordering on the nodes, often denoted by -,, such that the existence of a directed edge Xx + Xj implies that Xx. A X2 in the ordering.
Let P denote the collection of all permutations of the indices f1...., p}. For an arbitrary matrix A and any P. let P. represent the matrix obtained by permuting the rows and columns of A according to such that (P=
(x(())())-
DAG B is said be compatible with permutation if B is lower-triangular matrix, which is equivalent to saying that X - Xy in B implies that π (k) > π (j). Similarly, π is also called compatible with B.
For any positive definite matrix & and π E P, the matrix P-O represents the same covariance structure as O, up to reordering of the variables. The Cholesky decomposition of (+) can be uniquely written as:
P -I- (I-L)) L)
(),
where L is strictly lower triangular and D is diagonal. By Lemma 8 in [2], the following holds: E P.
P-O(L, =
Therefore,
D)-1D).
For each æ, we define:
B(τ) Ps LL, ((x)) = Ps 1D.
This suggests that for any E P, there exists pair (B(r), M(()) E(O), where B(x) can be uniquely determined based on the permutation π and θ [2]. is important to emphasize that different permutations, T1 * π2, can still result in the same pairs, i.e., (B(x1), ((T)) (B())),(TT) = (B(rg), δ (mT2)). Furthermore, this indicates that for any (B, E E()), there exists at least one permutation i such that (B, .2) = (T),, ((T)). Moreover, it turns out that the collection of pair of (()),((T))) forms the entire equivalence class E()).
Lemma (Lemma 1, 2). Suppose E is a positive definite covariance matrix and θ = Σ-1 Then,
E()) -P :) P.O =
-((B()),)(m)): : PP
This result indicates that the size of E()) is at most pl, which large but finite.
C.3
LS loss vs. log-likelihood
The following examples show that when the variances are unequal, the LS loss will not in general have the same minimizers as the log- likelihood. The first example just Example [33]. Suppose (X1, X2) is distributed according to the following linear SEM with unequal variances:
X1 2
X1 C1,
X₂
+
E2,
E1
N(,,
62 N(0,1/4)
Thus
0 0
-1/2) 0
11
0
-1/2 (-11/2
B0
Ω
E = Ey(B".
0
1/4)
-1/2 1/2
and also
E(e0) = Emin (09) = (B", BB}
24

---

where
B1 = (- Q ( 9)
Moreover, E(B", (2 I, Q1), since both SEM represent the same covariance.
But it turns out that B|||X BJ x| B[||X - (B9) T x(1): More precisely, it is easy to check that
B||| - B x1 = TT((T - B1) T - B1)) = 1, B||| (B") T x Tr((T- (B) ET B))) 5/4,
and moreover Bi is the global minimizer of the LS loss B[||X - B X ]. It follows that when the variances are different, the log-likelihood and LS loss have different global minimizers.
Similar calculations can be carried out for d > 2. but are tedious owing to the size E()). For example, here is an example of SEM over 3 nodes such that the LS loss has a different set of global minimizers, but also the LS-global minimizer has more edges than the sparsest Markov representation:
0 0
-
3/110
7 0
00
7
0
-2
Bd
)
0 0000
2
Ω
0
Σ ===E,h EB .29)
3
-5
0
0
0
2
2
- 10
For this model, LS loss selects the following SEM with 3 edges:
0
B
- 1.199
0 1.589
0.7733
We have 1a E Emin (89), but Bi & (0%)
C.4 Estimation bias under &
We provide example showing that when the l penalty is applied, the estimation becomes biased. Therefore, l should not be used. Consider the following linear Structural Equation Model (SEM):
Xi = NNO, 1) X2, = X1 X1+N
If the topological sort is known, i.e., X1 X2, and & penalty is used for minimizing the negative log-likelihood:
min log E| ax log | | + A/a|
Ideally, we would expect to be the minimal solution to the loss function. However, the problem is equivalent to:
log((1 a)² + o²) + A/a
is clear that nonzero for any À > 0, indicating that the é penalty leads to biased estimator. This bias does not occur when using MCP or SCAD with appropriate hyperparameters.
not the minimal solution to the loss function, as the derivative at a is
C.5 quasi-MCP, MCP and SCAD
We present the formulas for quasi-MCP, MCP, and SCAD, and demonstrate that quasi-MCP and MCP are equivalent.
pI,(()= 1,.(0==1(11))
1(4 < 6) + 1(4 > 0)
quasi-MCP
28
12
2aa 2
MCP
Pxa
MCP(() = 1(() < aA) ||tl|
+ +((() 2 Xa)
2a
2aa/1 - 12 X2
2² +
SCAD
Pxa
SCAD (t) = 1|//1(1| < 2) + 1(A 18
< aa)
1(t > Ao)
2(a 1)
25

---

Plot of Ps.
Figure 4: The plot of pt,)(t) with A = 2,8 = 1
is worth noting that we set as aA in quasi-MCP, then px,x = Px, (t). In another way, in MCP, then MCP PX. (t) Px,a(t). (f). Thus, quasi-MCP and MCP equivalent each
MCP
other.
C.6 Standardization of X and X
We present the formulas for the standardization of X and the standardization of the corresponding dataset X.
Let of == Var(X:) and D = dieg(a)....,op) Denote the standardized version of X as X which can be expressed as:
X D 1 IEX)
For X E R2xp, we can write X = (Xiz) = (xI,..., xn) = 1X1,... (1,,,, and define sample average for node Pij
VjE bll
i=1
Next, we define the sample variance for node j as
[x. 22
X
n
The diagonal matrix of sample standard deviations then
D diag ( a,)
Finally, we standardize X by subtracting the sample means and scaling by the inverse of D:
Z = X 1, (f....., D-1
where 1, E IR" an n-dimensional vector with all entries equal to 1.
C.7 Discussion of Assumption A,
In this subsection, provide more detailed discussion of Assumptions A and B.
26

---

First, it is important to emphasize that if px.5 is replaced by the to penalty in (12) or (14), then Assumptions A and can be omitted, and all the results still hold. In this case, the proof would be significantly simplified. However, the use of the differentiable quasi-MCP, in contrast to the to penalty, introduces substantial complications, necessitating some additional assumptions that are fundamentally different from the results in [8]. Specifically, Assumptions A and B are exactly what is required to make the problem suitable for gradient- based optimization.
Assumption A is not restrictive and is satisfied by all identifiable models, including linear Gaussian models, generalized linear models with continuous output [66]. binary output [74, 4], and most exponential families. However, the requirement for the finiteness of the equivalence class can be relaxed. What truly needed is that the minimal nonzero edge has sufficient "signal," i.e.,
min (W.E)EE(W9 (0) ((.)):(())
min
[B()) > 0
This is trivially true when [E(b". is finite. When (((.%.0) is infinite, each [B could be positive, but it is possible lim (V. ga) min(( (( can be arbitrarily small. The to penalty deals with this with its discontinuity at zero, whereas the continuity of quasi-MPC makes this more challenging. This is the cost of differentiability, which we argue is worthwhile
[B()) 0, because [Bg))
standard assumption in the optimization literature [7] and is generally quite weak. Moreover, it is almost necessary because quasi-MCP cannot exactly count the number of edges in (()). The magnitude of the quasi-MCT penalty does not directly reveal the number of edges. This is sparsity-inuccing the trade-off for replacing the to penalty with fully differentiable sparsity- inducing penalty. Finally, it is worth noting that this assumption can also be relaxed: what is truly required is that for any > there exists d' > 0 such that
Assumption
(( ((V. 50) >
for all ((v,E) dist((ib, <), E. (0)))) > c).
In other words, we require loss gap when () is not in ((bb,,))) This can be inferred from Assumption B.
The following is the proof when Px.s is replaced with the lo penalty.
Proof when Px,5 is replaced with bo:
{0) other words, there is unique element in the minimal equivalence class. This is because any element in (1/ ( is indistinguishahle based on the score function, i.e., the value of E(w,, E) and the penalty for the number of edges in (1)). Our objective is to find this unique element by solving equation (12) in the paper, which simplifies the proof.
Proof. We can also assume that (v (0) = 1, meaning Emin 59)
When 8 BB)") 0, the result is straightffrward:
ll (0)) + SB((O) S ((),,) + 8B())-
Now, let us consider the more general case where 8 BB(0) > 0 and divide the parameter space into three regions:
A1 = {()+E) (n()) > )n(wn) Ag = {(b.g) (n()) = {n As = {(b)g) sn(v)) < $n(w))}
Case 1: Consider A1. Since t(150 ()) s ((v,,)) the following holds for any À >
(0 + \ < E(1, { + A8 V( E A1.
Case 2: Consider Ag. Since |Emiu (v", " 1, it follows that for all (v) E(wo, gd) and (1, E) + (160,9), we have SH(W) > ((())) Therefore, for all (w.E) E Ag. it holds that ( 50) < E(1b, Consequently, for any > 0:
1(00 E0) + As < ((b,E) + ASB(y) V(W,E) E Ag.
Case 3: Consider Ag. We need to prove that:
((00 + (SB(Do) < E(v,,) + ASB(v) V(A,E) E As.
27

---

This is equivalent to showing that there exists positive À such that:
X (( E0
v(thE) As.
BB())) BB(e)
Since Sm(0o) is the minimal number of edges the equivalence class, any (v, () E A3 corresponds (1)) with a number of edges strictly less than SB(0o)- This implies that:
{( () El
v(O,E) A3.
Furthermore, we have S BB()0) BB() S BB(a)), which implies that there exists small but positive À that satisfies the inequality.
D Experiment Details
In this section, we provide all the details about the experiments. These include: (1) the types of graphs used. (2) the process for generating the samples, (3) the baseline methods we compare against and where to find the code for these methods, (4) the implementation details of our method and how replicate the results, and (5) the metries used evaluate the estimation.
D.1 Experimental Setting
In this section, we outline the process for generating graphs and data for Structural Equation Models (SEMs) (2). For each model, random graph generated using one of two types of random graph models: Erds-Rényi (ER) or Scale-Free (SF). The models are specified to have, on average, kp edges, where & E 21, 4]. These configurations are denoted as ERK or SFk, respectively.
Erds-Rényi (ER), Random graphs whose edges are add independently with equal probability, We simulated models with p, 2p and 4y edges (in expectation) each, denoted by ERI, ER2, and EA respectively.
Scale-free netw Network simulated according to the preferential attachment process [3]. We simulated scale-free network with p, 2p and 4p edges and where is the exponent used the preferential attachment process.
Linear SEMs. Given random DAG B (0, 1 from of these two graph models, edge weights were assigned independently from Unif([-1.5, -0.5] U [0.5, 1.5)) to obtain a weight matrix B E pp Given B, we sampled X B X + e R" according to:
Gaussian noise with unequal variance (Gauss-NV): 2 N(0,,7),, = I,...,p where 0, Unif(0.1,0.7]
We chose set o, the noise variances in our models, be relatively smaller compared to the settings used in previous studies such as [73], [41], and [4]. This decision aims to mitigate the potential exploitation of accumulated variance along the topological sort, highlighted in [55].
Generalized Linear Model with Binary Output Given a random DAG B E (0. 1 from one of these two graph models, edge weights were assigned independently from Unif([-1.5, -055 [0.5, 1.5]) to obtain weight matrix B E RRx,, Given B, we sample X, according to the following
X, Bernoull(cxxp(B] x)/(1 cxp(BB X)))
1=
where Bj j-th column of B. The corresponding negative log- ligeliild function:
p I (log(1,, (1, + exp(XB)) X, (XBB)
sB:
n i=1
where x= =Xx;) = (x1, (Xi...., Xn) (X
X,)
28

---

Nonlinear Models with Neural Networks. We primarily follow the nonlinear setting described in Zheng al. [74]. Given G, we simulate the SEM as follows:
Xj = (X()))))) + N2
Vj E [p].
where Nj ~ _N(0, o and o; 1n, Here, fj is a randomly initialized MLP with one hidden layer of size 100 and sigmoid activation. It is worth noting that the score function used in nonlinearNNOTEARG [74] is least square loss:
s sf. X
=
x
F(
2n
1=1
where each S is an MLP with one hidden layer of size 30 and sigmoid activation.
We generated random datasets X E xR* by sampling rows i.i.d. from the models deseribed above. For each simulation, produced datasets with n samples across graphs with p nodes.
Simulation
Linear Model: p {10, 20, 50, 70, 100). k {1,2, 4}, 1000 and graph types (ER, SF).
Generalized Linear Model: p p= {10, 20, 40]. k = [1,2], n = 10000 and graph types (ER, SF).
Nonlinear Model: p {10, 20, 403, {1, 2). n 1000 and graph types (ER, SF).
For each dataset, applicd several structural learning algorithms, including fast greedy equivalence search (FGES [52]), constraint based methods (PC [60]), NOTEARS [73, 74] (using least squares loss), GOLEM [41] (using NLL with & penalty), VarSort [55], causal additive models (CAM [9]), LOGGLL(-NOTAARS/DAGGO SSAMLL (utilizing the sample covariance matrix E), LOGLL(- (using the population covariance matrix Σ) and exact method (EXACT-SEARCH) Implementation details are provided in the following paragraph. After running the algorithms, post-processing threshold of 0.3 was applied to the estimated matrix Bea to prune small values, following the methodology in [73, 74].
Implementation
The implementation details of baseline are listed below:
Fast Greedy Equivalence Search (FGES [52]) is based on greedy search and assumes linear dependency between variables. The implementation is based the py-tetrad package, available at score function with default parameters.
com/ http://mmplpyyyttteeeeeee
We use BIC as the
PC [60] is constraint-based method and based uses conditional independence induced by causal relationships to learn those causal relationships. The implementation is based on the py-tetrad package, available at bettps://gttae com/cmu--piilii patnuuuhhillppyttteee We use Fisher-Z test with a 0.5.
NOTEARS(73, 74] is the continuous DAG learning algorithm using least square loss with I regularization. implemented in python: /ttps://githie com/ wwwzz/ngnnnoneeeeeeeeeeee
GOLEM [41] is implemented using Python and TensorFlow. The code is available https://gits thtb.b.t.huuuue/ com/ http://nnaainnng/geeeeee wereeegogggggeeeeeeee golem.
VarSort [55] is based on the observation that variances tend to accumulate along the topolog- ical sort. It uses Lasso [61] to recover the coefficients. The code implemented in Python and is available https: : / /////ttuu
DAGMA[4] is continues DAG learning algorithm with better accuracy and faster computa- tional speed. It also use least square loss with & penalty as NOTEARS. The implementation is available at https:: // gi thub. com com/kevvnsbelll// dagma.
Causal additive model (CAM [9]) learns an addititve SEM by leveraging efficient nonpara- metric regression techiques and greedy search over edges. The code implemented in R, and avaiable at https: //rr/ io/ ctto/rra///e Cat////AA/////AC/ html
29

---

is our approach, which modifies the original NOTEARS or DAGMA algorithm by replacing its scoring function. Instead of using the least squares loss with an & penalty, it employs a log-likelihood function that includes quasi-MCP penalty, as defined in (10). For (OGILL
SAMPLE, we use the sample covariance matrix E in the score function. In contrast, LOGLL(- COTTARRSDAAMMAA-PPPUAAATT
uses the true covariance matrix Σ as a baseline approach. In this paper, LOGLL-NOTEA refers to solving (12) using NOTEARS with NLL and quasi- MCP as the score function, while LOGLL-DAGAN refers to solving (12) using DAGMA with NLL and quasi-MCP as the score function.
EXACT-SEARCH is used to indicate that the optimization problem (12) is solved exactly. This approach is feasible only for small graphs, where we attempt to calculate all possible configurations (TT) defined in Section C.2. These calculations can be performed using Cholesky Decomposition or Ordinary Least Squares (OLS). The label POPULATION signifies that the operation is based on the population covariance matrix Σ, while SAMPLE denotes that it is based on the sample covariance matrix Σ. The Structural Hamming Distance (SHD) for EXACT-SEARCH is calculated on average basis. This involves identifying the set .Mmiu (+) Muui (+), calculating the SHD for each DAG within this set, and then computing the average SHD.
D.2 Implementation of LugeLL/AAOTEAAASSSSAAAAAApppappn
Linear model There are two main challenges solving (14). The first challenge is that (12) is a highly nonconvex optimization problem and is sensitive to initialization. If we randomly initialize or set the initialization to zero, as done in [73, 74, 41], LOGLL-NOTEARS/DAGAAA
often gets stuck at a local optimal solution. The second challenge arises from Theorem 1, where we advised to select and such that < < Ag and 0 < 8 < dg. Theoretically, smaller values for À and 8 should be used adhere the theorem's 's guidelines. however, in practice, solving the optimization problem (12) to global optimality is not always feasible.
To address the first challenge, we adopt the approach from Ng et al. [41]. We first run NOTEARS (with least squares loss and l penalty) or DAGMA (with least squares loss and & penalty) to obtain "good" initialization point. Then, we apply (OOLLLLONTAARS/AAGMM.. POPUAATTTOOAAAAAAAPPL.L. to obtain the final output.
To address second challenge, we use warm starts. We begin with larger values for À and d and solve (12) using LOGLL-NOTEARS/DAGMA to obtain an initial Beg. We then reduce À and 8 by factor of < 1 and use the previous output as the starting point for the next iteration of LOGLL- NOTEARS/DAGMA. This process is repeated until the negative log- likelihood fn (Bea Ress) begins to increase. This iterative approach helps to refine the solutions gradually, ensuring that each step starts from a potentially better approximation, as formally outlined in Algorithm
In Algorithm 1, we detail the complete implementation of replacing E with Σ and substituting f, with e, this algorithm is adapted to the full implementation of DOGLIL-XOTEARSIIAAAAAAPPPPU..L
By
It turns out that LOGLL-DAGMA outperforms LOGLL-NOTEARS in our experiments. Therefore, we present the results from LOGLL-DAGMA.
Nonlinear model
We utilize LOGLL-NOTEARS which uses the same optimization framework in [74] with replacement of least square loss to negative log-likelihood loss, we keep other hyperparam- eter unchanged. The score function (negative log- likeliheeeod use
SNLL(f,X)
log
[|x - F(())))
2n
Here f1 is i-th MLP with one hidden layer of size 40 and sigmoid activation.
Standardized data Z Although it has been shown that the log-likelihood score call--ivarriant for the linear model with Gaussian noise (see Theorem 3), it was observed that using standardized
30

---

Algorithm 1: Full Implementation
Input:
Sample covariance Σ, decay factor (0, 1), and X, 6, initial point (Bm: Pmn), initial loss ln (typically very large)
// initial point (1)) is obtained from NOTEARS DAGMA Output: (Beat Qes))
while True do Solve LOGIL-NOTEARS or LOGLL-DAGMA with input (Bi, Pnn), and get output (Bout, Pour) Calculate fn (Bout: Sout)
if lm (Bur, Dow)) then
Ln t (Bout: Qow))
À t 7X
8
7d
(Bmn lmn) + (Bout Poat)
else
10
return (Bin Din))
11
end
12 end
data Z makes solving the optimization problem (12) significantly more challenging. For the LOGLL- NOTEARS implementation, the LBFGS-B algorithm fails to produce meaningful solutions. As result, we replaced LBFGS-B with ADAM [27], an optimizer better suited for handling the difficulties of standardized data, to solve the subproblem in NOTEARS. Alternatively, directly using LOGLL- DAGMA another effective way to address this challenge. Empirically, find that setting 7 0.8, À = 0.4, and 8 0.2 usually serves as good choice for the parameters in our optimization procedures.
D.3 Metries
We evaluate the performance of each algorithm with the following three metrics:
Structure Hamming distance (SHD): standard benchmark in the structure learning literature that counts the total number of edges additions, deletions, and reversals needed convert the estimated graph into the true graph. Since our model specified in (6) is unidentifiable, the Structural Hamming Distance (SHD) is calculated with respect to the completed partially directed acyelic graph (CPDAG) of the ground truth and Bes- We utilize the code from Zheng et al. [75].
Times:
The amount of time the algorithm takes to run, measured in seconds. This metric is
used to evaluate the speed of the algorithms.
31

---

E Additional Results
E.1
Linear Model (SHD)
PC
Counthuuus
Combinatorial
FGES
raw: standardizzd
NOTEARS
GOLEM
Method
DAGMA
VarSort
LogLl_sample
LogLL_population
Exact_sample
Exact_population
5
10 15
20
25
Structural Hamming Distance (SHD)
Figure Structural Hamming Distance (SHD, with lower values indicating better performance) between Markov equivalence classes (MEC) of recovered and ground truth graphs for ER-2 graphs with 8 nodes. Here Exact-search is added to illustrate Theorem 3. Standardization does not affect the DAG structure if the optimization (10) can be solved globally. Both Exact-sample and Exact- population produce the same DAG structure for raw data X and standardized data Z. When the population covariance matrix is known, Emia (89) = -M(G)), resulting in SHD of zero. The poor performance of Exact-sample can be attributed to the lack thresholding applied to the coefficients recovered from Ordinary Least Squares (OLS). Since Σ is only an approximation of Σ, coefficients derived from OLS based on different permutations may shift from zero to nonzero, even though such coefficients might be very small. However, since Exact is impractical for real-world applications, we use this example primarily for illustrative purposes, and thus no threshold is applied this method.
32

---

PC
FGES
standard
NOTEARS
Method
GOLEM
DAGMA
VarSort
LoglL_sampll
2
Structural Hamming Distance (SHD)
Figure 6: Comparison of raw (orange) vs. standardized (green) data. Structural Hamming Distance (SHD, with lower values indicating better performance) between Markov equivalence classes (MEC) of recovered and ground truth graphs for ER-2 graphs with 5 nodes
Logei
Buseline
Countiouus
Combinatoo
FGES
standarrtired
NOTEARS
Method
GOLEM
DAGMA
VarSort
LoglL_sample
populaoopplation
20
40
60
80
100
120
Structural Hamming Distance (SHD)
Figure 7: Comparison of raw (orange) standardized (green) data. Structural Hamming Distance (SHD, with lower values indicating better performance) between Markov equivalence classes (MEC) of recovered and ground truth graphs for ER-2 graphs with 20 nodes
33

---

0=-100-2.1-93---5
40000
SHD builwenn MLC
-00000
20000
10000
0.2
0.4
0.5
0.>110-a2.1-9332-55
50000
400000
Distnnce
20000
10000
0.2
0.4
0.6
0.3
1.0
maximal
Figure 8: Graph: structure Xo - Xo X2, X1 estimated (Bea: Dea)) E Emin(() because SHD and distance are closed to O.
X2 + X3. For 0 < 8o, the
Structure: 8.2-17
100000
SHD belwenn MEC
-60000
400000
-3.5
-1.0
maximal + 2hd
Structure: 0-> 1-2->1>>1
100000
-60000
of
20000
0.8
Figure 9: Graph: structure X₂ - X1, X1. For 0 < 8 < do. the estimated (Bes, Rea)
Emin (0") because SHD and distance are closed to O.
34

---

E.2 Linear Model (Time)
500
400
400
600
300
300
400
200
200
Time (seconds)
100
100
0-
500
400
400
400
300
200
200
100
0-
0-
100
nodes)
nodes)
[Number nodes]
method
DAGMA
NOTEARS POPULATION
LOTEARS NOTEARS SAMPLE
NOTEARS
Figure 10: Results in term of Time. Lower is better. Column: k {1, (1. Row: random graph types. JER,SFi-k = Frade--FFee,rss----e Rényi } graphs with kd expected edges. Here d {10, 20, 50, 70, 100), n 1000. Standard error is removed for better visualization. is for different methods on raw data X
600
400
200
200
200
Tine (secands)
0
600
400
400
200
200
200
0
50
100
25
100
50
75
100
(Number of nodes)
(Number of nodes)
NNumber of nodes)
LOGLL, COTTARSSPOPPULTTOO
LOTEL_ SSANPLE
Figure 11: Results in term of Time. Lower is better. Column: = (1, Row: random graph types. |ER,SF)-K = Seaale-Free,Erds-PPP {10, 20, 50, 70, 100), n 1000. Standard error is removed for better visualization. It is for different methods standardized data Z
} graphs with kd expected edges. Here d
35

---

E.3 Nonlinear Model (SHD)
E.3.1
Neural Network
MEC Structural Hamming Distance Method
pushh
gruph guahh
12, gegr
gogth Sk-k
26 Bki
goth whhhshkk2
20
CAM
gosh
guzh
Figure 12: Structural Hamming distance (SHD) between Markov equivalence classes (MEC) of recovered and ground truth graphs. LOGLL (i.c. LOGLL-NOTEARS) stands for NOTEARS method with log-likelihood and quasi-MCP, L2 (i.e. NOTEARS) stands for NOTEARS method with least square and f1-
E.3.2 General Linear Model with Binary Output (Logistic Model)
General Linear Model with Binary Output Model)
83
mn Distine (5+D)
10
(Number of modes)
(Number of nodes)
(Number of nodes)
method
NOTEARS
NOTEARS, LOGLL LOGLL
Figure 13: Structural Hamming distance (SHD) for Logistic Model, Row: random graph types, (SF, ER f== (Scale-Free, LOGLL-NOTEARS uses logelikel with quasi-MCP, NOTEARS use log logelikeli with l, Error bars represent standard errors over 10 simulations.
graphs. Columns: kd expected edges. NOTEARS LOGLL (i.c.
36

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
We are confident about this point.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: This is the last section of main paper.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
37

---

Justification:
We include all the important assumptions.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We provide the detailed discussion in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
38

---

Answer:
[Yes]
Justification: is simple adaption of open source code. We will release our the code.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.o cc/
jublic/guides/ htt///uddde//cc/ http///ooeeeaaiiioiinnoao ooiicyyyyyy)o. for more details.
While we encourage the release of code and data, we understand that this might not be possible, so No an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
httpa///ggiieeeeeeeuubiiiiiii / /nipp.cc/// wttpccccppppccccceecc.ccccccccct pttpuulccuuuids..../ /Codeeeooiiisaeeeeeee wwsissssssaaeeeee onnooicy))) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: is included in the experinments section.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification:
We include the error bar to illustrate this point.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
39

---

It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of enrors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative eтror rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [No]
Justification: We run the algorithm on personal computer. It should easy to replicate the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https : http///eerrpppeeeeeeeeeeeeeeeeee cc wwtpppplpccceeeeeeeeeeeeeeeeeeee wttpub/cc//tiisttiiie htt.ii//eeeccccccteeeeeeeeee
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g.. if there is special consid- cration due to laws or regulations their jurisdiction..
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification:
We mainly focus on the theoretical findings.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
40

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This is unrelated to the topic of paper.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: We cite the related paper.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
41

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new assets, Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
42


Reviews:
[{"review_id": "p4iWqpT0TM", "reviewer": "Reviewer_o9h2", "summary": "The authors analyze a framework of sparsity-regularized maximum likelihood learning under the NOTEARS constraint for score-based causal discovery. Drawing from the sparsest permutations principle (a hybrid method), they show that a sparsity regularized likelihood objective is able to recover an element of the MEC even under structural non-identifiability of the SCM. The (non-identifiable) linear Gaussian setting is analyzed in detail, where it is shown that the MCP and SCAD penalties are able to recover the sparsest graph with appropriate hyperparameter settings. The authors then extend this to general likelihoods under the assumption that the (parameter) non-identifiability class (at the ground truth) is finite and other regularity conditions. For the linear-Gaussian case, the authors also prove scale invariance of the structure obtained from the method, addressing concerns of varsortability that similar methods are susceptible to. Finally, an experimental study is conducted on simulated linear and nonlinear ground truths.", "questions": "### Other questions  \n\n- You mention that $\\ell_1$ loss is \"not effective\" compared to MCP, SCAD, is this a computational or theoretical issue? In other words, do you have a conjecture/result on whether the theoretical claims no longer hold when replaced with $\\ell_1$?\n- Although the likelihood is \"structurally\" scale-invariant (Thm. 3), the scale presumably can still affect (uniform over d) thresholding in finite samples by changing the speed at which structural zeros converge. This has always been my intuition on how var-sortability works, and not that the MSE in NOTEARS is actually \"structurally\" affected by scaling. Is there a reference to show that the original MSE loss is not also structurally scale-invariant? If the MSE is also structurally scale-invariant, do you conjecture that the likelihood has desirable properties in terms of thresholding?", "limitations": "The authors do not adequately discuss the limitations of the paper in my opinion (in fact, the paper abruptly ends without even a conclusion), and this is also seen in in the checklist where the justification for this section was left blank. \n\nI think the paper deserves a proper conclusion and discussion of the (potential) limitations of Assumption A1 in Section 5. __I do not think this is a tough ask, so I may lower my score if this remains unaddressed in the rebuttal.__", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "I advocate for the acceptance of this paper based on two simple dimensions:\n\n- (Significance and Motivation): Structural identifiability is a near-universal, yet entirely convenience-driven assumption when using likelihood-based scores for causal discovery. This work theoretically justifies NOTEARS-type approaches under structural non-identifiability, which broadens the class of applicable likelihood-based scores.\n\n- (Clarity and Execution): I found the paper easy to follow, with precise mathematical notation, careful proofs and comprehensive analysis of the proposed framework.", "weaknesses": "However, I do think the paper falls somewhat short in a few (fixable) areas. Actionable questions are __bolded__, and I may raise my score if these points are clarified. \n\n- The section on scale invariance seems out of place and somewhat weaker than the rest of the paper. It's clear that the Gaussian likelihood is scale invariant, and structure invariance is not a terribly surprising conclusion of the regularized objective also. __Unless the authors can clarify the contribution of this section, I feel like this is better off stated as a short paragraph or remark.__\n- The authors state that the finite parameter equivalence class assumption in the general case is \"relatively mild\", and that it is satisfied by a range of models (\"most exponential families\"). I'm not convinced that this is the case. First, typically when statistical models with continuous parameter spaces are non-identifiable, the equivalence class is usually infinite. For example, I believe exponential families in canonical parametrization are either full-rank, where the parameter is identifiable (and hence not relevant to the motivation of the paper), or non-identifiable up to entire subspaces of the parameter space (Theorem 1 of [1]), which are infinite. __Could the authors clarify what they mean on l297-298?__\n- The authors do not provide any examples, or references, of non-Gaussian log-likelihoods scores that would be useful for causal discovery. The non-linear example in the experiments seems to use a MSE loss which is the same as non-linear NOTEARS [2], and the only improvement shown in the experiments is due to changing the optimization scheme to Adam (and adding the quasi-MCD regularizer?). __Could you provide (practical) examples of likelihoods where the theory in Section 5 hold? I.e., other non-identifiable likelihood scores that are useful for causal discovery?__. \n\n[1] Notes by Charles Geyer: https://www.stat.umn.edu/geyer/8053/notes/expfam.pdf\n\n[2] \"Learning Sparse Nonparametric DAGs\", Zheng et al., AISTATS 2020.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to express our gratitude to the reviewer for acknowledging the value of our contributions and clarity of presentation. All the points addressed below will be updated in our paper.\n> The section on scale invariance ….  or remark.\n\nWe will move these details to the appendix, and shorten the discussion in the main paper.\n\n> The authors state that the finite parameter…..mean on l297-298?\n\nFirst, consider simple case of Gaussian models, and then below we discuss how to generalize. Let $X \\sim \\mathcal{N}(0, \\Sigma)$ and $\\Sigma\\succ0$. It belongs to the exponential family and is full-rank. The covariance $\\Sigma$ is identifiable, however, we are interested in the matrix $B$ in the SEM $X = B^\\top X + N$. As in Section 4.1, we know that $\\Sigma = (I - B)^{-\\top} \\Omega (I - B)^{-1}$. As we discuss in Section 4.2, there are at most $p!$ different $B$ that satisfy this condition, making our problem **unidentifiable with respect to $B$, but with a finite number of equivalent parameters**.\n\n**As for unidentifiable exponential families in [4],** we can always reduce it to an identifiable model by projecting canonical statistics to a subspace. For the same example, if $\\Sigma$ is singular, we can decompose $X = (X_a, X_b)$ such that $\\text{Cov}(X_a)\\succ 0$  and has the same rank as $\\Sigma$, and $X_b$ is a linear combination of $X_a$. As discussed before, there is a finite equivalence class for $X_a$, and $X_b$ has a deterministic relationship to $X_a$.\n\nAnother explicit example (GLM with binary data) can be found in the **Common Concern (1)**. This also extends to more general exponential families. When exponential families are full-rank, although their parameters are identifiable, the parameters $(\\psi, \\xi)$ as well as the SEM $B(\\psi)$ may not be identifiable, which are our interest. As before, parameter is $\\Sigma$, $\\psi = B, \\xi = \\Omega$, and $\\Sigma$ is function of $(B,\\Omega)$. Moreover, in many cases, the size of $\\mathcal{E}(\\psi^0, \\xi^0)$ is closely related to $p!$, which is finite since there are $p!$ topological sorts $\\pi$. For each $\\pi$, it is usually associated with a pair $(\\psi(\\pi), \\xi(\\pi))$ in a certain way such that $P(x; \\psi(\\pi), \\xi(\\pi)) = P(x; \\psi^0, \\xi^0)$. This is why we claim Assumption A(1) is reasonable.\n\nRegarding Assumption A and B, more discussion can be found in the **Common Concern (3).** \n\n> The authors do not provide…scores that are useful for causal discovery?.\n\nThe point is that *any* likelihood can be used in conjunction with the quasi-MCP. Two examples of non-Gaussian likelihoods are logistic (binary) models and additive noise models. A third example is any model for the joint distribution with GLMs for CPDs. All told, these comprise a very rich class of models with explicitly derivable score functions. Unfortunately, we have not made this point explicit and will certainly do so in the camera ready.\n\nFor additive noise models, the log-likelihood function is $\\frac{1}{2n}\\sum_{i=1}^d \\log (\\\\|x_i -\\hat{f}\\_i(x)\\\\|^2)$[1], which is used in our experiments (see Line 855 in Appendix). It is different from MSE loss $\\frac{1}{2n}\\sum_{i=1}^d\\\\|x_i - \\hat{f}_i(x)\\\\|_2^2$ used in NOTEARS [2]. For the improvement of the experiments, it is based on the correct use of loss function and penalty, thus it leads to better performance. For other examples, please refer to **Common Concern (1).**\n\n> You mention that $\\ell_1$ loss ... when replaced with $\\ell_1$?\n\nUsing $\\ell_1$ as a penalty can lead to biased parameter estimates because it generally shrinks coefficients to zero. In contrast, MCP provides unbiased estimates. Specifically, with MCP, $\\mathcal{O}_{n,\\lambda,\\delta} = \\mathcal{E}\\_{\\min}(\\Theta)$ when population case is considered. However, this is not true for $\\ell_1$.\n\nConsider $X_1\\sim \\mathcal{N}(0,1)$, $X_2\\sim X_1+\\mathcal{N}(0,1)$, with the known order ($X_1\\rightarrow X_2$), so that the log-likelihood with $\\ell_1$ penalty is\n$$\n\\log((1-a)^2+1)+\\lambda|a|\n$$\n\nIt is easy to see $a = 1$ is not minimizer. This indicates that $\\ell_1$ leads to biased estimates. \n\n> Although the likelihood is \"structurally\"….. in terms of thresholding?\n\nGood point! Although scaling does not affect which entries in the adjacency matrix are zero, it changes the magnitude of the entries. With a finite number of samples, this alters the thresholding we should use, which significantly affect the recovered structure in practice.\n\nWith respect to MSE loss, it depends: It is scale-invariant for a fixed topological sort, but this is not true in general when the sort is unknown. Let us clarify:\n\n1. For fixed topo sort, the MSE in NOTEARS is scale-invariant. Since the matrix $B$ can be recovered by linear regression, and changing the scale of the data $X$ does not change the positions of the nonzero entries in $B$, leading to invariance. This is similar to the intuition of Theorem 3.\n2. However, the MSE loss in NOTEARS is not scale-invariant when the sort is unknown. For instance, consider $X = (X_1,X_2)$ where $X_1 = N_1, X_2 =aX_1+ N_2, N_1, N_2\\overset{i.i.d}{\\sim}\\mathcal{N}(0,1)$. Then the optimal for MSE loss is\n$$\nB^* = \\begin{bmatrix} 0 & a \\\\\\\\ 0 & 0 \\end{bmatrix}\n$$\nHowever, if you standardize $X$ to get $Z =(X_1,X_2/\\sqrt{1+a^2})$ and optimal are\n$$\nB^*_1 = \\begin{bmatrix}0&a/\\sqrt{a^2+1} \\\\\\\\ 0&0\\end{bmatrix}\\qquad B^*_2 = \\begin{bmatrix}0&0 \\\\\\\\ a/\\sqrt{a^2+1}&0\\end{bmatrix}\n$$\nSo MSE is not scale-invariant in the sense that global optimal solution can not have the same structure for different scale of data. However, this is not the case for likelihood (Theorem 3).\n\n> Limitation and conclusion\n\nPlease refer to **Common Concern (4)**. More discussions on Assumption A1 can be found in previous response. We will add this to the paper with the extra page in the camera ready.\n\n[1] Bühlmann, et al \"CAM: Causal additive models, high-dimensional order search and penalized regression\" (2014)\n\n[2] Zheng, et al. \"Learning sparse nonparametric dags\" (2020)", "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": "ht96ErwZZQ", "reviewer": "Reviewer_spWh", "summary": "This paper introduces new identifiability results (MEC) based on maximum likelihood estimation complemented with sparsity regularization (quasi-MCP) for both a Gaussian linear model and a more general, potentially nonlinear models, under the very standard faithfulness assumption. The paper contains a theoretical analysis showing the score function is invariant to rescaling (in the Gaussian linear case only) and provide experiments showing the approach compares favorably to other baselines (in the Gaussian linear case only).\n\n**Review summary:**\nThis paper reads very nicely and some of its proposals are new and interesting (scale-invariant result and the usage of the quasi-MCP to get MEC identifiability), but it missed a very relevant work [a] with very similar contributions (which significantly reduce the novelty factor of this manuscript) and a more complete set of experiments (see below). For these reasons, I believe this work is not ready yet for publication. That being said, I think it could get accepted at another conference by\n- rewriting the paper in a more transparent way, especially by contrasting with the contributions of [b] (and changing the title); \n- highlighting what are the actual contributions (scale-invariant result + MEC identifiability via quasi-MCP regularizer); and \n- providing convincing experiments with nonlinear models.", "questions": "See above.", "limitations": "This paper lacks a conclusion/discussion, so the limitations of the approach are not addressed.", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "- This paper is extremely well written and easy to follow. The notation is always well introduced with kind reminders when less standard notation is used. \n- Overall the paper feels quite pedagogical.\n- I believe the theoretical analysis of scale-invariance, although simple, is novel and very interesting, as it addresses problems raised with existing approaches regarding standardization.\n- I think some aspect of the theoretical results is interesting, like the proof that the quasi-MCP regularizer can be used to get identifiability of the MEC. However, I believe very similar results have been shown previously in the literature (omitted from the related work) which seriously limit the novelty factor of this work.", "weaknesses": "**Important prior work omitted which seriously limits novelty**\n\nThe contribution is motivated by the need to consider more general score functions beyond the mean-squared error used in NOTEARS. The authors missed the work of [a] which proves identifiability for very general likelihood functions, including universal approximators (thanks to normalizing flows). The result of [a] applies to cases where interventions are observed, but covers also the special case without interventions, in which case the result guarantees identifiability of the MEC under the faithfulness assumption. Interestingly, [a] also requires the coefficient of the regularizer to be “small enough”, as is also the case in this manuscript. A key distinction with [a] is that this manuscript is analyzing a differentiable nonconvex sparsity regularizer as opposed to an L0 regularizer, which is not suitable for gradient-descent (in practice, [a] uses gumbel-sigmoid masks and regularizes their probability of being one). Am I missing other key distinctions? Showing that the quasi-MCP regularizer can yield MEC identification is new AFAIK and interesting IMO, but this combined with the simple scale-invariance results makes for a rather limited contribution...\n\nI believe the title reflects the lack of knowledge of the existing literature on continuous structure learning. The title is too general. It sounds like the paper introduces likelihood-based differentiable structure learning, but that’s clearly not the case. [a] is a clear example with very similar theory (although not exactly encompassing this manuscript) as well as [b] and [c] (although without identifiability theory).\n\n**Experiments investigate only the linear case**\n\nThe paper is largely motivated by the need to go beyond linear Gaussian score functions, but the experiments fall completely short of that promise, as only the Gaussian linear case is investigated. In contrast, [a] train neural network based architecture, with similar theoretical guarantees. More experiments are needed to confirm that the quasi-MPC regularizer transfers to nonlinear setting. \n\n**Limitations of the theory:**\n\nAssumption A: Assuming that the cardinality of equivalence of parameters is finite feels like a strong assumption. For instance, if you’d like to parameterize the conditional using neural networks, this assumption wouldn’t be satisfied since most architecture have infinitely many equivalent parameterization (think of ReLU activation where you can rescale the output of a neuron and undo this rescaling in the following layer). The result of [a] does not make this assumption despite their very similar setup. Why is it needed here? Also, can you give a concrete example of function B(\\psi) that is Lipschitz.\n\nAssumption B: The same example with ReLU neural networks is a counterexample to this assumption. Indeed, because of this rescaling, the set of parameters yielding the same network is unbounded. Excluding NN feels limiting, especially since existing results can cover this case.\n\n**Unclear limit statements in Theorem 1 to 4**\n\nThe usage of the limit symbol “->” is a bit strange. Usually, we say that a_n -> b as n -> \\infty, but here the result says a_n = b as n -> \\infty. What does the latter even mean? Does that mean there exists a large enough N such that a_n = b for all n > N? I don’t think it’s true anyway. Why not presenting the result directly for the population likelihood (i.e. with \\ell instead of \\ell_n), as was done in [a]? Skimming the proof very quickly suggests that’s what was proved anyway. Note that this is different from showing consistency of the procedure, i.e. that your estimator converges to the right MEC as the number of samples grows (I suspect the latter will be more difficult to show)\n\nMinor:\n- It would be nice to have a small plot of the quasi-MCP regularizer.\n- Typo in Lemma 3 in appendix? I think an \\mathcal{M} should be replaced by a \\mathcal{G}, no?\n- Text is absurdly small in Figure 1 e.g.. Please fix this.\n- Is the sparsest Markov representation always unique? The phrasing used suggests it is, but it might be worth expliciting.\n\n[a] P. Brouillard, S. Lachapelle, A. Lacoste, S. Lacoste-Julien, and A. Drouin. Differentiable causal discovery from interventional data. In Advances in Neural Information Processing Systems, 2020.\n\n[b] X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing. Learning sparse nonparametric dags. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, 2020.\n\n[c] S. Lachapelle, P. Brouillard, T. Deleu, and S. Lacoste-Julien. Gradient-based neural DAG\nlearning. In Proceedings of the 8th International Conference on Learning Representations,\n2020.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the insightful critiques and comprehensive understanding of our work, and for providing such useful feedback. We will try our best to address the reviewer’s concern.\n\n> **Comparison to previous work**\n> \n\nWe thank the reviewer for highlighting this related work; we will certainly update our paper to cite it and provide a discussion. Although our work and [1] prove similar results on identifiability through optimizing the log-likelihood with sparsity regularization, there are important differences. Most importantly, the use of $\\ell_0$ significantly simplifies the analysis (similar identifiability results as Theorem 1, 2, and 4 are easily obtained, see **Common Concern (3)** for details) and ultimately does *not* lead to a differentiable program. Our focus is on a fully differentiable formulation. In fact, this is precisely what leads to Assumptions A and B, which are *not* needed when $\\ell_0$-regularization is used. This is not surprising since $\\ell_0$ leads to a combinatorial optimization problem which we are trying to avoid. Again, see **Common Concern (3)** for details.\n\nMoreover, our results rely on slightly weaker assumptions to get identifiability results, i.e., Sparsest Markov Representation (SMR) assumption, which has been shown to be weaker than faithfulness (with only observational data, $\\mathcal{I}^*$-faithfulness assumption is regular faithfulness assumption [1]). Our proof is entirely different and more straightforward, relying solely on the score function, without the consideration of different skeletons and immoralities [1]. Our paper also includes new results on scale-invariance to address current concerns with many DAG learning algorithms. \n\n> **Experiments investigate only the linear case**\n> \n\nOriginally, due to space limits, we put the nonlinear experiments results in Appendix D.3.3. (Page 31), which may have been easy to miss. Nonetheless, we have added more nonlinear experiments; please see **Common Concern (2)**.\n\n> **Limitations of the theory**\n> \n\nThanks for these useful comments! Assumptions are needed because the our penalty term is quasi-MCP, introducing fundamental difficulties compared to previous work [1] using $\\ell_0$. If penalty term is $\\ell_0$, all these assumptions can be removed, and similar results hold. Please refer to **Common Concern (3)** for more details.\n\n> **Also, can you give a concrete example of function B(\\psi) that is Lipschitz.**\n> \n\nFirst, in a (general) linear model, $\\psi = B$. In this sense, $B(\\psi) = \\psi = B$, and $B(\\psi)$ is 1-Lipschitz. Another example can be found in **Common Concern (1).** In this case, $\\psi = \\\\{\\beta^S\\\\}_{S\\subseteq [p]}$, $[B(\\psi)]\\_{ji} = \\sum\\_{\\{S:j\\in S, i\\in S, S\\subseteq A\\}}(\\beta^{S})^2$, which is Lipschitz.\n\nMoreover, [2] uses a neural network to model $f_j$ in Equation (2). Let $\\psi_j$ be the weights that connect the input layer and the first hidden layer in the $j$-th neural network. In this case, $\\psi = (\\psi_1, \\ldots, \\psi_p)$. It turns out that $[B(\\psi)]\\_{ij} = \\\\|\\text{i-th-column}(\\psi_j)\\\\|_2$. As a consequence, $\\\\|B(\\psi)\\\\|\\_2 = \\sum\\_{i,j} \\\\|\\text{i-th-column}(\\psi_j)\\\\|_2 \\leq p \\sum_j \\\\|\\psi_j\\\\|_2 \\leq p^2 \\\\|\\psi\\\\|_2$. Therefore, it is $p^2$-Lipschitz. In general, $B(\\psi)$ is a function of $\\psi$, and usually, the nonzero part of $\\psi$ indicates that a certain entry of $B(\\psi)$ is nonzero. $\\ell_1$ or $\\ell_2$ norms are used to ensure nonzeros are mapped to nonzeros and zeros are mapped to zeros, which makes $B(\\psi)$ inherit the Lipschitz property of $\\ell_1$ or $\\ell_2$.\n\n> **Unclear limit statements in Theorem 1 to 4**\n> \n\nThank you for the suggestion, and we apologize for any confusion. Indeed, our results apply to the population case, which we indicate by using $n \\rightarrow \\infty$. For Theorems 1 to 4, we are stating that when considering the population case (infinite number of samples), certain equality relationships hold. We will revise these statement to make this more clear.\n\n> **Minor**\n> \n\n(1) Originally, we include the plot of quasi-MCP, but it is removed due to space limit. But it seems it will be better idea to get it back again.\n\n(2) Yes, this is a typo, it should be $\\mathcal{M}(G^0) = \\mathcal{G}(\\mathcal{E}_{\\min}(\\psi^0,\\xi^0))$\n\n(3) Good suggestion! it will be fixed.\n\n(4) Sparsest Markov representation(SMR) is not unique itself in general, but it is unique up to Markovian class. \n\n> Limitation\n> \n\nThanks for bring it up. We provide conclusions and limitations in **Common Concern (4).**\n\n[1] Brouillard, Philippe, et al. \"Differentiable causal discovery from interventional data.\" *Advances in Neural Information Processing Systems* 33 (2020): 21865-21877.\n\n[2] Zheng, Xun, et al. \"Learning sparse nonparametric dags.\" *International Conference on Artificial Intelligence and Statistics*. Pmlr, 2020.", "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": "llOiOsPC7o", "reviewer": "Reviewer_mGjz", "summary": "This paper proposed a differentiable DAG learning method using the log-likelihood loss and the minimax concave penalty (MCP). The author proved that under such construction, the minimizer of the loss identifies the sparsest graph (i.e. it has the minimal number of edges) which can generate the observational distribution. When faithfulness is further assumed, the identified graph is equivalent to the ground truth graph. The theorems are valid for linear Gaussian DAGs and general nonlinear SEMs where the induced distribution is parametric. The proposed method is validated by extensive experiments involving various dimensions and graph structures.", "questions": "Can you provide further examples/illustrations on whether BIC or $L_1$ penalty can recover minimal models? I think this will strengthen the proposed method when faithfulness fails.\n\nThe log-likelihood for nonlinear models seems improper due to the heteroscedastic noise. Either the setting or the likelihood should be corrected before rebooting the experiments.\n\nIt would be interesting to compare the proposed method with the existing ones on real datasets.", "limitations": "The method assumes a known likelihood function, which can be impractical when the noise distribution or the functional structure is unknown.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and easy to follow. \n\nThe theoretical contribution is solid, demonstrating that a broad class of penalties (SCAD, MCD, quasi-MCD), along with the log-likelihood loss, identify the desired graph structure.\n\nThe simulations are well-designed and illustrate the outperformance of the proposed method.", "weaknesses": "The method can be regarded as a direct combination of several works. For example, the log-likelihood loss was also proposed in GOLEM [43]. The MCD penalty was first proposed in [72].\n\nWhen applied to nonlinear models, the author used $\\frac{1}{2n}\\sum_{i=1}^d\\log\\left(||x_i-\\hat f_i(X)||^2\\right)$ as the log-likelihood loss. This seems only valid for homogeneous Gaussian noise where the variance in the denominator can be canceled out. However, in the experimental setting, the noise variance is heteroscedastic, making the likelihood improper.\n\nReal-world datasets can be added to validate the proposed method.\n\nCode is not available.\n\nSome minor flaws in the presentation:\n+ The conclusion section is missing.\n+ Section D.3.2, time spent on the simulation. The plots are about SHD rather than about the running time.\n+  [72] in the reference should be cited correctly.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to express our gratitude to the reviewer for their time, effort, and valuable suggestions. We would like to take this opportunity to address all the concerns raised in the reviews.\n\n> The method can be regarded as a direct combination of several works. For example, the log-likelihood loss was also proposed in GOLEM [43]. The MCD penalty was first proposed in [72].\n> \n\nThanks for reviewer’s questions on our contribution. We would like to emphasize that our paper is not merely a combination of these previous works.\n\n(1) **New Identifiability Result:** We utilize our log-likelihood loss with a sparsity regularizer (MCP) to prove an identifiability result under common assumptions from the optimization literature. Unlike GOLEM [1], which has specific assumptions (such as the Triangle Assumption [1]) and relies on a linear model, our method does not require these conditions. Moreover, our proof is purely based on the properties of the score function and does not rely on any assumptions about the underlying graph structure, unlike GOLEM. Additionally, as far as we know, this is the first time MCP has been introduced in the continuous DAG learning literature. Previously, $\\ell_1$ was the commonly used sparsity regularizer, which led to biased estimation of the underlying parameters and could negatively affect the recovered structure. By utilizing quasi-MCP, we overcome these drawbacks of $\\ell_1$ and achieve unbiased estimation. Moreover, quasi-MCP also introduces significant technical challenges that our analysis handles (see also **Common Concern (3)**).\n\n(2) **Scale-invariant property:** We address a common concern for many DAG learning algorithms regarding scale invariance, wherein the algorithm may output different DAGs if the data is rescaled. To address this, we demonstrate that the log-likelihood score is scale-invariant under Gaussian noise. By showing this result, we aim to correct the common misunderstanding that the algorithm itself is affected by the scale of the data. Instead, it is the use of an incorrect score function that is sensitive to the scale.\n\nWe believe these points distinguish our contributions from previous works, and are significant.\n\n> When applied to nonlinear models, the author used $\\frac{1}{2n}\\sum_{i=1}^d \\log (\\\\|x_i -\\hat{f}_i(x)\\\\|^2)$ as the log-likelihood loss. This seems only valid for homogeneous Gaussian noise where the variance in the denominator can be canceled out. However, in the experimental setting, the noise variance is heteroscedastic, making the likelihood improper.\n\n> The log-likelihood for nonlinear models seems improper due to the heteroscedastic noise. Either the setting or the likelihood should be corrected before rebooting the experiments.\n\nThank you for this important question! We are happy to clarify the misunderstanding and will include the full derivations in the paper for clarity.\n\nWhen applied to a nonlinear model with heteroscedastic noise, the negative log-likelihood is: \n\n$$\n\\frac{1}{2n}\\sum\\_{i=1}^d \\log (\\\\|x_i -\\hat{f}_i(x)\\\\|^2)\n$$\nHowever, when the noise is homogeneous, the negative log-likelihood becomes: \n$$\n\\frac{d}{2}\\log \\left(\\sum\\_{i=1}^d \\\\|x_i - \\hat{f}_i(x)\\\\|^2\\right)\n$$\nAlthough these formulations appear similar, they are fundamentally different. With heteroscedastic noise, the loss is the log of the sum of least squares. In contrast, with homogeneous noise, the loss is the sum of the logs of least squares. The full derivation can be found in [1] (which covers both heteroscedastic and homogeneous noise in the linear case, with $\\hat{f}_i(x) = B_i^\\top x^{(i)}$) and [2] (which uses heteroscedastic noise in their setting).\n\nThanks for this question! We address this point later in the Question “The log-likelihood….before rebooting the experiments.”\n\n> Real-world datasets can be added to validate the proposed method.\n> \n\n> It would be interesting to compare the proposed method with the existing ones on real datasets.\n> \n\nWe include the real-data application in attached pdf (See **Common Concern (2)**). These results will be updated into our paper.\n\n> Code is not available.\n> \n\nThank you for bringing this up, the code will be released in the camera ready. \n\n> The conclusion section is missing.\n> \n\nThanks for pointing it out. Please refer to **Common Concern (4)** for conclusion and limitation.\n\n> Section D.3.2, time spent on the simulation. The plots are about SHD rather than about the running time.\n> \n\nThanks for bring up this typo. We put the wrong name in the title and y-axis, it should be title: “$d$ vs Time for methods”, y-axis: \"Time (second)”. The value in the plot is correct. This will be update in the paper correspondingly.\n\n> [72] in the reference should be cited correctly.\n> \n\nThanks; this reference will be updated in the final version.\n\n> Can you provide further examples/illustrations on whether BIC or $\\ell_1$ penalty can recover minimal models? I think this will strengthen the proposed method when faithfulness fails.\n> \n\nThanks for the suggestion. Please refer to **Common Concern (2)** for more examples.\n\n[1] Ng, Ignavier, AmirEmad Ghassami, and Kun Zhang. \"On the role of sparsity and dag constraints for learning linear dags.\" *Advances in Neural Information Processing Systems* 33 (2020): 17943-17954.\n\n[2] Bühlmann, Peter, Jonas Peters, and Jan Ernest. \"CAM: Causal additive models, high-dimensional order search and penalized regression.\" (2014): 2526-2556.", "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:
**(1) Another example beyond the linear gaussian model.**
Note that Assumption A holds as long as $P(X_i|X_A)$ has a unique SEM parametrization for any $i$ and $A$, where $A\subseteq[p]\backslash i$.  Because for any fixed topological sort $\pi$, we could let $A = \\{\text{parents node of }i \text{ in }\pi \\}$, and this results in at most $p!$ equivalence class. This is for example true in the Gaussian model.

Another example is a generalized linear model with binary output, i.e., $X = (X_1,\ldots,X_p)$ and $X_i\in\\{0,1\\}$ for $i = 1,\ldots,p$. Let $B = (B_1,\ldots,B_p)$, then  $\mathbb{E}[X_i\mid X_{pa(i)}] = g(B_i^\top X)$ where $g(s) = e^s/(1+e^s)$ which is equivalent to the following SEM

$$
X_i = \text{Bernoulli}(\exp(B_i^\top X)/(1+\exp(B_i^\top X)))\qquad i = 1,\ldots,p
$$

In this case, it can be shown that

$$
X_i\mid X_A \sim \text{Bernoulli}\Big(\frac{1}{1+\exp(- \sum_{r =1}^{|A|+1}(\sum_{j_1 = i,j_2,\ldots,j_r\in A}\beta^{ij_2\ldots j_r}x_{j_2}\ldots x_{j_r}))}\Big)
$$

is the unique parametrization of the conditional distribution (in terms of the coefficients $\beta^{ij_2\ldots j_r}$, which are a function of $B$). Thus, binary models also satisfy Assumption A. We will add this example to the paper.

**(2) More experiment on Nonlinear Neural Network, and General Linear model with binary output**

We conducted more experiments on Nonlinear Neural Networks (see Appendix D for details), and generalized linear models with binary output. The results are attached in the pdf.

**(3) Assumption A, B for Theorem 4**

First of all, we would like to emphasize that when using $\ell_0$ as the penalty term, Assumptions A and B are not needed at all. The proof for this case would be significantly simplified (see below for a proof sketch). Replacing $\ell_0$ with the differentiable quasi-MCP introduces significant complications, requiring some additional assumptions that distinguish our work from [1]. In fact, our assumptions are exactly what is needed to make the problem amenable to gradient-based optimization.

Considering Assumption A, the finiteness can be relaxed. What is really needed is that the minimal nonzero edge has enough “signal,” i.e.,

$$
\min\_{(\psi,\xi)\in \mathcal{E}(\psi^0,\xi^0)}\min\_{\\{(i,j):B(\psi)_{ij}\ne 0\\}}|B(\psi)|\_{ij}>0
$$

This is trivially true when $|\mathcal{E}(\psi^0,\xi^0)|$ is finite. When it is infinite, each $|B(\psi)|\_{ij}$ could be positive, but it is possible $\lim\inf_{(\psi,\xi)\in \mathcal{E}(\psi^0,\xi^0)}\min\_{\\{(i,j):B(\psi)_{ij}\ne 0\\}}|B(\psi)|\_{ij}=0$, because $|B(\psi)|\_{ij}$ can be arbitrarily small. The $\ell_0$ penalty deals with this with its discontinuity at zero, whereas the continuity of quasi-MPC makes this more challenging. This is the cost of differentiability, which we argue is worthwhile.

As for Assumption B, this is a common assumption in the optimization literature and quite weak in general. Moreover, this is nearly necessary because quasi-MCP cannot count the number of edges in $B(\psi)$ exactly: The magnitude of the quasi-MCP penalty does not reveal the number of edges. This is the price to pay when we replace $\ell_0$ with a fully differentiable sparsity penalty. Finally, we can point out that this can also be relaxed: What is needed is that for any $\epsilon>0,$ there exists  $\delta>0$

$$
\ell(\psi,\xi)-\ell(\psi^0,\xi^0)>\delta\quad \text{for }\\{(\psi,\xi)\mid \text{dist}((\psi,\xi),\mathcal{E}(\psi^0,\xi^0))>\epsilon\\}
$$

In other words, it requires a loss gap when $(\psi,\xi)$ is not in $\mathcal{E}(\psi^0,\xi^0)$.

Finally, we include a brief proof sketch of a similar result to Theorem 4 when $\ell_0$ is used. This illustrates that the quasi-MCP indeed introduces fundamental difficulties. Moreover, the reason the proof simplifies compared to [1] is a) The proof below does not consider interventions (which is a major innovation of [1]), and b) The use of the weaker SMR assumption (vs faithfulness) simplifies the analysis.

**Proof sketch:** We can also assume $|\mathcal{E}\_{\min}(\psi^0,\xi^0)|=1$, following the same reason as for Theorem 4.

When $s_{B(\psi^0)} = 0,$ the result is obvious. When $s_{B(\psi^0)}>0$ , divide the parameter space into $A_1 = \\{(\psi,\xi)\mid s_{B(\psi)} > s_{B(\psi^0)}\\}, A_2 = \\{(\psi,\xi)\mid s_{B(\psi)} = s_{B(\psi^0)}\\}, A_3 = \\{(\psi,\xi)\mid s_{B(\psi)} < s_{B(\psi^0)}\\}$. It is straightforward to verify each case. For example, for $A_2$, since $|\mathcal{E}\_{\min}(\psi^0,\xi^0)|=1$, that implies $\forall(\psi,\xi)\in\mathcal{E}(\psi^0,\xi^0)$ and  $(\psi,\xi)\ne (\psi^0,\xi^0)$, it holds $s_{B(\psi)}>s_{B(\psi^0)}$. Therefore, $\forall (\psi,\xi)\in A_2$, it holds $\ell(\psi^0,\xi^0)<\ell(\psi,\xi)$. As consequence, for any $\lambda>0$, result holds.

The other cases are similar with slight modifications.

(4) **Missing conclusion and limitation:** We proposed a fully differentiable score function for causal discovery, composed of log-likelihood and quasi-MCP. We demonstrated that the global solution corresponds to the sparsest DAG structure that is Markov to the data distribution. Under mild assumptions, we conclude that all optimal solutions are the sparsest within the same Markov Equivalence Class. Additionally, the proposed score is scale-invariant, producing the same structure regardless of the data scale under the linear Gaussian model. Experimental results validate our theory, showing that our score provides better and more robust structure recovery compared to other scores.

However, there are limitations to our work. We focus on parametric models and rely on assumptions such as the finiteness of the equivalence class and the boundedness of the level set of the log-likelihood. These assumptions limit the applicability of our theorem. Future work should explore ways to relax these assumptions to extend our method’s applicability to broader scenarios.

[1] Brouillard, et al. "Differentiable causal discovery from interventional data." (2020)

Meta-review: This paper presents a fully-differentiable approach for DAG structure learning based on the log-likelihood loss and quasi-MCP (minimax concave penalty) regularization. The main contributions are in terms of identifiability results (the method identifies the sparest graph that is Markov to the data distribution) and scale invariance properties (the proposed method produces the same DAG regardless of the data scale). 

Several concerns were raised during the discussion period regarding (1) experiments only investigating the linear Gaussian case; (2) the omission (and lack of discussion) of closely related work by Brouillard et al (NeurIPS 2020); (3) the limitations of the theory (3) and the lack of conclusions in the manuscript and discussion of limitations of the proposed approach. I believe the authors have provided a comprehensive rebuttal addressing the above and the reviewers have acknowledged that in their comments and scores updates.

Therefore, I recommend acceptance.
Decision: Accept (poster)