Markov Equivalence and Consistency in Differentiable Structure Learning

Authors: chang  deng
kevin  bello
pradeep  ravikumar
bryon  aragam
Venue: NeurIPS_cc_2024_Conference
Published: 2024-04-23 23:19:34

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, thus paving the way for differentiable structure learning under general models and losses.

Main paper: 
 introduction: 
 Directed acyclic graphs (DAGs) are 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 scorebased approaches to DAG learning [11, 23] , where the structure learning problem is formulated as optimizing a given score or loss function s(B; X) that measures how well the graph, represented as an adjacency matrix B ∈ {0, 1} p×p , fits the observed data X, constrained to the graphical structure B being acyclic. This combinatorial optimization problem is generally known to be NP-complete [10, 12] .Recent advances in score-based methods have introduced a 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, 1} p×p is first relaxed to the space of real matrices, i.e., B ∈ R p×p , and then a differentiable function h : R p×p → [0, ∞) is devised so that h(B) = 0 if and only if B is a DAG [73, 4] . This results in the following optimization problem: min B∈R p×p s(B; X) subject to h(B) = 0.(1) Considering a differentiable score function s, the differentiable program (1) facilitates the use of gradient-based optimization techniques along with the use of richer models, such as neural networks, 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 functions under this framework, we note: [41] studied the properties of the ℓ 1 -regularized profile log-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 an ℓ 0 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-Softmax) 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 ℓ 0 function, it provably leads to useful solutions that are minimal models and Markov equivalent to the underlying structure, without assuming faithfulness. 
preliminaries: 
 We let G = (V, E) denote a directed graph on p nodes, with vertex set V = [p] := {1, . . . , p} and edge set E ⊂ V × V , where (i, j) ∈ E indicates the presence of a directed edge from node i to node j. We associate each node i ∈ V to a random variable X i , and let X = (X 1 , . . . , X p ).). An SEM (X, f, P (N )) over the random vector X = (X 1 , . . . , X p ) is a collection of p structural equations of the form:EQUATIONwhere f = (f j ) p j=1 is a collection of functions f j : R p+1 → R, here N = (N 1 , . . . , N p ) is a vector of independent noises with distribution P (N ), and PA j denotes the set of parents of node j. Here, ∂ k f j denotes the partial derivative of f j w.r.t. X k , which is identically zero when f j is independent of X k , i.e. f j (X, N j ) = f j (PA j , N j ). The graphical structure induced by the SEM, assumed to be a DAG, will be represented by the following p × p weighted adjacency matrix B: The SEM (2) is general enough to include many well-known models, such as linear SEMs [e.g., 33, 49] , generalized linear models [47, 45, 18] , and additive noise models [24, 51] , post-nonlinear models [70, 71] and general nonlinear SEM [38, 21, 26, 74] . To illustrate some of these models: In linear SEMs we have X j = f j (PA j ) + N j , where f j is a linear map; in causal additive models (CAM) we have X j = k∈PAj f j,k (X k ) + N j , where f j,k is a univariate function; in post-nonlinear models we have X j = f j,1 (f j,2 (PA j ) + N j ). In fact, essentially any distribution can be represented as an SCM of the form (2) ; see Proposition 7.1 in Peters et al. [50] .EQUATIONFaithfulness and sparsest representations. It is well-known that the DAG G is not always identifiable 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 implioed by the distribution P , and let I(G) denote the set of d-separations implied by the graph G. Then P is Markov to G if I(G) ⊂ I(P ), and faithful to G if I(P ) ⊂ I(G). When both conditions hold, i.e. I(P ) = I(G), then G is called a perfect map of P . Following common convention, we will simply call P faithful when I(G) = I(P ). When P is faithful to G, the Markov equivalence class (MEC) of G is identifiable and can be represented by a CPDAG. Since faithfulness may not always hold, there has been progress in understanding what can be identified 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 a 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 4 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, s, we will use ψ ∈ Ψ ⊆ R m and ξ ∈ Ξ ⊆ R s to denote the model parameters for f = (f 1 , . . . , f p ) and N , respectively. 2Then we denote the distribution of X by P (X; ψ, ξ). Let x ∈ R p denote one observation of X. Given n i.i.d. samples X = (x 1 , . . . , x n ) ⊤ where x i ∼ P (X; ψ, ξ), the negative log-likelihood and expected negative log-likelihood can be written as:EQUATIONwhere the subscript n in ℓ n is used to indicate the sample version of the log-likelihood.Identifiability. Let ψ 0 (resp. ξ 0 ) denote the model parameters for the ground truth f 0 (resp. N 0 ), let B 0 = B 0 (ψ 0 ) ∈ R p×p denote the induced weighted adjacency matrix, and let G(B 0 ) ∈ {0, 1} p×p denote the induced binary adjacency matrix. For example, in the general linear Gaussian model ( 6), ψ = B represents the adjacency matrix, and ξ = Ω denotes the variance of the Gaussian noise. In another case, if f j is approximated by a multilayer perceptron (MLP), with N j as Gaussian noise, then ψ includes all the parameters of the MLP, while ξ represents the variance of the Gaussian noise. Additionally,EQUATIONj ∥, where Aj is the first hidden layer in f j [74] . Thus, by our definitions, P (X; ψ 0 , ξ 0 ) is the true distribution. Here, there are two types of identifiability questions:1. Parameter identifiability: Is it possible to uniquely determine the parameters (ψ 0 , ξ 0 ) based on observations from P (X; ψ 0 , ξ 0 )? Formally, is there any ( ψ, ξ) ̸ = (ψ 0 , ξ 0 ), such that P (X, ψ 0 , ξ 0 ) = P (X, ψ, ξ) almost surely? 2. Structural identifiability: Is it possible to uniquely determine the DAG G(B 0 ) based on observations from P (X; ψ 0 , ξ 0 )? In other words, is there any( ψ, ξ) ̸ = (ψ 0 , ξ 0 ) such that P (X, ψ 0 , ξ 0 ) = P (X, ψ, ξ) but G(B 0 ) ̸ = G(B( ψ)).In general, parameter identifiability implies 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.e. structural identifiability does not always imply parameter identifiability, as different parameter values can lead to the same structure. Classical results on identifiability of SEMs include: linear SEM with equal variance [33] , linear SEM with 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 ℓ(ψ, ξ) serves as a natural choice for the score function because it attains a unique minimum at the true parameters (ψ 0 , ξ 0 ). However, this approach is not straightforward for nonidentifiable models, where multiple parameter sets can induce the same data distribution P (X; ψ 0 , ξ 0 ), leading to ambiguities in parameter or structure estimation. In such cases, regularizing the log-likelihood 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. 
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 to theoretically show how the NLL with nonconvex differentiable regularizers can lead to meaningful solutions such as minimal-edge models and elements of Markov equivalent classes. We also discuss and prove the scale invariance of NLL, 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.A linear SEM (B, Ω) over X with independent Gaussian noises N , a special case of (2), is wellknown to be nonidentifiable in terms of parameters and structure [see 2, for discussion]. We write the model as follows:EQUATIONwhere B ∈ R p×p is a matrix of coefficients with G(B) being a DAG, and N ∈ R p is the vector of independent noises with covariance matrix Ω = diag(ω 2 1 , . . . , ω 2 p ).3 Given the model (6) it is easy to see that the distribution P (X) is Gaussian and is fully characterized by the pair (B, Ω). That is:EQUATION)where Σ is the covariance matrix of X. In the sequel, we use the subscript f to refer to a function. In this case, Σ f denotes a function with arguments (B, Ω) and returns the covariance matrix. Moreover, we use Θ to denote the corresponding precision matrix (inverse of the covariance matrix):EQUATION)Let X = (x 1 , . . . , x n ) ⊤ be n i.i.d. samples of X. Then, let the sample covariance matrix beΣ = 1 n n i=1 x i x ⊤ i .The sample NLL function is given by:ℓ n (B, Ω) = - 1 n log n i=1 P (x i ; B, Ω) = 1 2 log det Ω -log det(I -B) + 1 2 Tr( ΣΘ(B, Ω)) + const.The corresponding population NLL function isℓ(B, Ω) = -E X log P (X; B, Ω) = 1 2 log det Ω -log det(I -B) + 1 2Tr(ΣΘ(B, Ω)) + const.The full derivation can be found in Appendix C.1. Here, it is important to note that the distribution of X is fully determined by either the precision matrix Θ or the covariance matrix Σ.Our goal is to identify (B, Ω): Unfortunately, the model is inherently nonidentifiable in terms of both parameter and structure. This means that multiple pairs (B, Ω) 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 E(Θ) as the set of all pairs (B, Ω)EQUATIONIt is worth noting that the size of E(Θ) is finite and at most p!, which corresponds to the number of permutations for p variables [2] . For more comprehensive details on this class, see Appendix C.2. This ambiguity naturally leads to the question: which pair (B, Ω) 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, a DAG that induces the precision matrix Θ 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 B = |{(i, j) :B ij ̸ = 0}|. Definition 2 (Minimality). (B, Ω) is called a minimal-edge I-map 4 in the equivalence class E(Θ) if s B ≤ s B , ∀( B, Ω) ∈ E(Θ).The set of all minimal-edge I-maps in the equivalence class E(Θ) is referred to as the minimal equivalence class E min (Θ):E min (Θ) = {(B, Ω) : (B, Ω) is minimal-edge I-map, (B, Ω) ∈ E(Θ)}.In the sequel, for brevity, we will often refer to such models as "minimal models".Unlike faithfulness, which may not always hold, the minimal equivalence class E min (Θ) 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, we have M(G) = E min (Θ), i.e., E min (Θ) is the Markov equivalence class of G. However, there could be multiple pairs (B, Ω) within E min (Θ). Nevertheless, our goal is to recover one element from E min (Θ). The elements in E min (Θ) 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 0 , Ω 0 ) andΘ 0 = Θ f (B 0 , Ω 0 ). Assume that P (X) is faithful to G 0 := G(B 0 ). Then M(G 0 ) = E min (Θ 0 ).Recall our convention that this means that M(G 0 ) = G(E min (Θ 0 )), i.e. the DAG structures contained in E min (Θ 0 ) coincide with M(G 0 ). Thus, under the faithfulness assumption, recovering E min (Θ 0 ) is the same as recovering the MEC, which is the usual goal in causal discovery. Moreover, we 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-likelihood when the noise variances are different; see Appendix C.3.In order to distinguish elements in E(Θ) from the minimal elements in E min (Θ), 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 equivalently the ℓ 0 penalty. Although both approaches effectively penalize the number of nonzero entries in B, their non-differentiability makes them unsuitable for differentiable structure learning. The ℓ 1 penalty, while amenable to differentiable approaches, 5 is not effective in precisely counting the number of edges, and also biased in parameter estimation 6 . 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 to use a reparametrized version of MCP, termed quasi-MCP, defined as follows:quasi-MCP:EQUATIONHere, 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, δ] and remains constant for values greater than δ. The function is smooth, and for values |t| > δ, it approximates the behavior of the ℓ 0 penalty, thus serving to penalize the number of non-zero coefficients in B.The score function in (1) can be naturally written asEQUATIONwhere p λ,δ (B) = i̸ =j p λ,δ (B ij ). Then, the optimization problem can be written asEQUATIONIt is worth noting that for any B, the corresponding optimal Ω that minimizes s(B, Ω; λ, δ, X) can be easily be expressed in terms of B as Ω f (B) (see Appendix C.1). Therefore, we can always plug Ω f (B) into (12) to profile out Ω.Solving problem (12) requires minimizing ℓ n (B, Ω) and p λ,δ simultaneously. To study the behavior of these minimizers, let us define the set of global minimizers,EQUATION)Ideally, we would like O n,λ,δ = E min (Θ 0 ), however, it is unclear whether there exist values of λ and δ such that any optimal solution (B * , Ω * ) lies within E min (Θ). The following theorem provides an affirmative answer to this question. In the sequel, we say that a property S(x) holds for all sufficiently small x > 0 if there is some fixed ϵ > 0 such that for every x ≤ ϵ, the property S(x) holds. Theorem 1. Let X follow model (6) with (B 0 , Ω 0 ) and Θ 0 = Θ f (B 0 , Ω 0 ). Let X be n i.i.d. samples from P (X), and O n,λ,δ be defined as in (13) . Then, for all sufficiently small λ, δ > 0 (independent of n), it holds that P (O n,λ,δ = E min (Θ 0 )) → 1 as n → ∞.In other words, we can always guarantee that O n,λ,δ = E min (Θ 0 ) by taking λ, δ sufficiently small, which is easily accomplished in practice. In the following, we use the superscript 0 to denote ground truth parameters. Additionally, we can assume that B 0 always belongs to E min (Θ 0 ), ensuring that our reference to the ground truth aligns with the simplest or minimal representation within the equivalence class. Moreover, by Lemma 1, under the faithfulness assumption, Theorem 1 can be interpreted as recovering the Markov equivalence class M(G 0 ): Theorem 2. Consider the setup in Theorem 1 and assume additionally that P (X) is faithful to G 0 := G(B 0 ). Then, for all sufficiently small λ, δ > 0 (independent of n), it holds thatP (O n,λ,δ = M(G 0 )) → 1 as n → ∞.Theorem 2 indicates with properly chosen hyperparameters, the optimal solution from optimization (12) will produce a graph that adheres to the same independence statements as G 0 . 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 1 in Appendix A for details.It is known that the LS loss is not scale-invariant, i.e. re-scaling the data (and in particular, standardizing it) can drastically change the structure [33] , a fact which Reisach et al. [55] use to argue that differentiable DAG learning with the LS Loss is also not scale-invariant. Here we show that by using a different score-in this case the log-likelihood-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 are invariant to re-scaling (i.e. re-scaling does not change the support for any (B, Ω) ∈ E(Θ)) using our notation: Lemma 2. Let X ∼ N (0, Σ), suppose Σ is a positive definite covariance matrix and let Θ := Σ -1 , suppose D is a diagonal matrix with positive diagonal entries. Then G(E(Θ)) = G(E(DΘD)).Lemma 2 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 whether X or Z is used as input, and Lemma 2 suggests that re-scaling X will not alter the structure of the DAG that is recovered from optimizing (12) . The following theorem formalizes this: Theorem 3. Under the same setting as Theorem 1, the solutions to (12) are scale-invariant. That is, for any n ≥ 0, letO n,λ,δ (X) ={(B * , Ω * ) : (B * , Ω * ) is a minimizer of (12) with data X}, O n,λ,δ (Z) ={(B * , Ω * ) : (B * , Ω * ) is a minimizer of (12) with data Z},where Z is the standardized version of X. Then, for all sufficiently small λ, δ ≥ 0 and all n, we haveG(O n,λ,δ (X)) = G(O n,λ,δ (Z)).Moreover, for all sufficiently small λ, δ > 0 we haveP G(O n,λ,δ (X)) = G(O n,λ,δ (Z)) = G(E min (Θ f (B 0 , Ω 0 ))) → 1 as n → ∞.Thus, even on finite samples, the set of DAG structures G(O n,λ,δ (X)) derived from the raw (unstandardized) data X will always be the same as G(O n,λ,δ (Z)), which is derived from standardized data Z. As a result, standardizing Gaussian data does not affect the recovered DAG structure if the optimization problem ( 12) can be solved exactly. Remark 2. Theorem 3 applies to global optimization of the objective (12) . Of course, in practice, algorithms can get stuck in local optima, but the global solutions (even for finite samples n) will always be scale invariant. 
nonconvex regularized log-likelihood for general models: 
 The results in the previous section are not specific to 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 X follows model ( 2) and the induced distribution is denoted by P (X; ψ 0 , ξ 0 ). Let us define the equivalence class E(ψ 0 , ξ 0 ),E(ψ 0 , ξ 0 ) = {(ψ, ξ) : P (x; ψ, ξ) = P (x; ψ 0 , ξ 0 ), ∀x ∈ R p }.That is, E(ψ 0 , ξ 0 ) is a set of pairs (ψ, ξ) that induce the same distribution P (X; ψ 0 , ξ 0 ). As a result, any pair (ψ, ξ) within this equivalence class will be a minimizer of the NLL ℓ(ψ, ξ). Analogously to Definition 2, we can also define the collection of minimal elements in the equivalence class E(ψ 0 , ξ 0 ).Definition 3. (ψ, ξ) is called a minimal-edge I-map in the equivalence class E(ψ 0 , ξ 0 ) if s B(ψ) ≤ s B( ψ) , ∀( ψ, ξ) ∈ E(ψ 0 , ξ 0 ). We further define E min (ψ 0 , ξ 0 ) = {(ψ, ξ) : (ψ, ξ) is minimal-edge I-map, (ψ, ξ) ∈ E(ψ 0 , ξ 0 )}.Here, it is crucial that our concept of minimality concerns s B(ψ) , which is the number of nonzero entries in the weighted adjacency matrix B(ψ), rather than the number of nonzero entries in the parameter ψ itself. Therefore, s B(ψ) essentially counts the number of edges in the adjacency matrix.Assumption A. (1) |E(ψ 0 , ξ 0 )| is finite. (2) B(ψ) is L-Lipschitz w.r.t. ψ, i.e. ∥B(ψ1)-B(ψ2)∥2 ∥ψ1-ψ2∥2 ≤ L.Assumption B. For any α such that ℓ(ψ 0 , ξ 0 ) < α, the level set {(ψ, ξ) : ℓ(ψ, ξ) ≤ α} is bounded, where ℓ(ψ, ξ) 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 2) is a mild continuity requirement on B(ψ). Assumption B 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.|E(ψ 0 , ξ 0 )| = 1. Assumption A(Similar in spirit to Theorem 1, we can show that by combining the NLL with quasi-MCP for appropriate λ, δ, solving the following problem, we recover elements of E min (ψ 0 , ξ 0 ):EQUATIONwhere p λ,δ (•) is quasi-MCP defined in (10) . Next, define its set of global minimizers.O n,λ,δ = {(ψ * , ξ * ) : (ψ * , ξ * ) is minimizer of (14)}. Theorem 4. Let X follow model (2)with parameters (ψ 0 , ξ 0 ) and let X be n i.i.d. samples from P (X; ψ 0 , ξ 0 ). Under Assumptions A-B, for all sufficiently small λ, δ > 0 (independent of n), it holds thatP (O n,λ,δ = E min (ψ 0 , ξ 0 )) → 1 as n → ∞.Theorem 5. Under the setting in Theorem 4 and assuming that P (X; ξ 0 , ψ 0 ) is faithful with respect to G 0 := G(B(ψ 0 )). Then, for all sufficiently small λ, δ > 0 (independent of n), it holds thatP (O n,λ,δ = M(G 0 )) → 1 as n → ∞. 
experiments: 
 To solve ( 12) and ( 14), we employ the augmented Lagrangian algorithm [5] from NOTEARS [73, 74] , modifying their least squares score with ℓ 1 penalty into the log-likelihood with MCP (10) . We compare our approach to relevant baselines, e.g. NOTEARS [73] , GOLEM [41] , DAGMA [4] , VarSort [55] , FGES [52] and PC [59] . For our variation of NOTEARS that employs a score function based on the NLL with MCP, we name it as LOGLL-NOTEARS. The suffixes 'POPULATION' and 'SAMPLE' denote the use of the population and sample covariance matrix, respectively. Full details of the experiments are given in Appendix D.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) is nonidentifiable, we compare the CPDAGs of the estimated graph and the ground truth graph.In Figure 1 (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 {ER1, SF1}. 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-SAMPLE and LOGLL-NOTEARS-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 1 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(-NOTEARS)-SAMPLE and LOGLL(-NOTEARS)-POPULATION exhibit the scale-invariant property alongside PC and FGES, in alignment 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 LOGLL(-NOTEARS)-SAMPLE and LOGLL(-NOTEARS)-POPULATION. This can be attributed to the increased complexity of optimization as the size of the graph grows.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 δ 0 that can be applied according to Theorem 1. We generate 10 5 initializations B random with weight for each edge uniformly sampled within [-5, 5], and perform optimization using LOGLL-NOTEARS starting from these points. The "maximal δ" is the theoretical maximum δ 0 that ensures the validity of Theorem 1. We computed the SHD and the distances between the estimated M(B est ) and M(B 0 ). The red line in Figure 3 represents the average SHD and distances. The distribution of these 10 5 estimated SHD and distances is visualized using dots of varying sizes, where larger dots indicate a higher frequency of points. In some cases where SHD takes a value of -1, this value is used to indicate that the estimated B est does not form a 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 j ) and neural networks. Further details on the experimental settings and additional experiments can be found in Appendix D and E. 
conclusion: 
 Continuous score-based 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 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.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 log-likelihood become more interesting in this regime. We have mentioned already that extensions to richer data types including interventions is an 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.  
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This is unrelated to the topic of paper.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 to 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 a best faith effort.12. Licenses for existing assets Question: 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?Answer: [NA] Justification: We cite the related paper.Guidelines:• The answer NA means that the paper does not use 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 a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.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, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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. 
 
 Appendix: Differentiable Structure LearningIn 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) as a 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 1 (MCP/SCAD). Under the same setting as Theorem 1. Let optimal solutions collection beEQUATIONThen, for all sufficiently small λ, a > 0 (independent of n), it holds that O n,λ,a = E min (Θ 0 ) as n → ∞, where MCP p M CP λ,a (•) and SCAD p SCAD λ,a (•) are defined in Appendix C.5.The following lemma is a generalization of Lemma 1. Even for the general model, under the faithfulness assumption, all elements in the minimal equivalence class E min (ψ 0 , ξ 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 (ψ 0 , ξ 0 ). Assume that P (X;ξ 0 , ψ 0 ) is faithful to G 0 := G(B(ψ 0 )). Then M(G 0 ) = G(E min (ψ 0 , ξ 0 )) = {G(B(ψ)) : (ψ, ξ) ∈ E min (ψ 0 , ξ 0 )}where B(ψ) is the adjacency matrix implied by the parameterization (ψ, ξ), see (3) . M(G 0 ) is the Markov equivalence class of G 0 , see Definition 1.Under the Sparsest Markov representation assumption, all elements in the minimal equivalence class are 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 ) satisfies the Sparsest Markov representation assumption. The following lemma provides the formulation for the standardization of X, along with its covariance and precision matrices. Lemma 5 (standardization). Let X ∼ N (0, Σ), σ 2 i := Var(X i ) and D := diag(σ 1 , . . . , σ p ). Then the standardization of X, corresponding covariance matrix and precision matrix can be expressed asX std := D -1 (X -EX), Cov(X std ) = D -1 ΣD -1 [Cov(X std )] -1 = DΘDThe following lemma establishes an useful identity that holds for any adjacency matrix of a DAG, which is used in the derivation of the log-likelihood function for the model in Equation (6) . Lemma 6. If B is adjacency matrix of a DAG, then log det(I -B) = 0.The following lemma provides a condition under which the optimization problem ( 12) is well-defined, ensuring that ℓ(B, Ω) > -∞ for any (B, Ω). Lemma 7. For any (B, Ω), if Ω > 0, then Σ := Σ f (B, Ω) is positive definite. Moreover, if X is generated by Equation ( 6) with (B 0 , Ω 0 ), then ℓ(B, Ω) > -∞ for any (B, Ω).The following lemma is used in the proof of Theorem 1. It justifies that the loss of every element in A 3 is strictly greater than the loss of the ground truth, i.e., (B 0 , Ω 0 ). Lemma 8. Under the same setting and notation as in the proof of Theorem 1, see Section B.1. If for any ( B, Ω) ∈ E(Θ 0 ), it holds that dist( B, A 3 ) > 0, there exists α > 0 such that ℓ(B, Ω)ℓ(B 0 , Ω 0 ) > α for all (B, Ω) ∈ A 3 . By convention, we assume that (B 0 , Ω 0 ) ∈ E min (Θ 0 ). Now, consider the case where p λ,δ (B 0 ) = 0, which is equivalent to B 0 = 0, since ℓ(B, Ω) ≥ ℓ(B 0 , Ω 0 ) and p λ,δ (B) ≥ p λ,δ (B 0 ) for any B. Therefore, for all λ > 0 and δ > 0, B 0 is the unique optimal solution to optimization problem (12) , proving the conclusion.In the subsequent proof, we assume that|E min (Θ 0 )| = 1, that is, E min (Θ 0 ) = {(B 0 , Ω 0 )}.This assumption simplifies the proof because any element of E min (Θ 0 ) is indistinguishable based on the value of ℓ(B, Ω) and the penalty for the chosen (δ, λ), as shown below. Our goal is to identify one element via optimization problem (12) , which significantly simplifies the argument.First, let us defineδ 0 = τ 1 + ∆where τ := min (B,Ω)∈E(Θ 0 ) min {(i,j)|Bij ̸ =0} |B ij | (a) = min π min {(i,j)|[ B 0 (π)]ij ̸ =0} [ B 0 (π)] ijp λ,δ (B) = λδ 2 s B > λδ 2 s B 0 = p λ,δ (B 0 ) ∀(B, Ω) ∈ E(Θ 0 ) \ {(B 0 , Ω 0 )}.As a consequence, this implies that ℓ(B 0 , Ω 0 ) < ℓ(B, Ω), ∀(B, Ω) ∈ A 1 . Therefore,ℓ(B 0 , Ω 0 ) + p λ,δ (B 0 ) < ℓ(B, Ω) + p λ,δ (B) ∀(B, Ω) ∈ A 1 .Next, we defineA 2 = {(B, Ω) | p λ,δ (B) > p λ,δ (B 0 )}. Since ℓ(B 0 , Ω 0 ) ≤ ℓ(B, Ω), it follows that for all (B, Ω) ∈ A 2 , the following inequality holds:ℓ(B 0 , Ω 0 ) + p λ,δ (B 0 ) < ℓ(B, Ω) + p λ,δ (B) ∀(B, Ω) ∈ A 2 .Therefore, we need to examine the setA 3 = {(B, Ω) | p λ,δ (B) < p λ,δ (B 0 )}.For (B 0 , Ω 0 ) to achieve the minimum value of the score function, it is crucial that the following condition is satisfied:ℓ(B 0 , Ω 0 ) + p λ,δ (B 0 ) < ℓ(B, Ω) + p λ,δ (B) ∀(B, Ω) ∈ A 3 .This condition guarantees that the ground truth parameters (B 0 , Ω 0 ) correspond to the optimal solution by comparing their score with any other parameters in the subset A 3 .It is important to note that p λ,δ (t) = λp 1,δ (t), ∀t. Thus, a necessary and sufficient condition for this to hold is:λ < min (B,Ω)∈A3 ℓ(B, Ω) -ℓ(B 0 , Ω 0 ) p 1,δ (B 0 ) -p 1,δ (B) .Note that for all (B, Ω) ∈ A 3 , we have p 1,δ (B 0 ) -p 1,δ (B) ≤ δ 2 s 0 , with equality achieved when B = 0. Therefore, the denominator on the RHS cannot be arbitrarily large. Moreover, since (B, Ω) ∈ A 3 , it follows that (B, Ω) / ∈ E(Θ 0 ), as A 3 ∩ E(Θ 0 ) = ∅.We define the distance from B to the set A 3 as:dist( B, A 3 ) = inf (B,Ω)∈A3 ∥B -B∥ 2 .For all ( B, Ω) ∈ E(Θ 0 ), it turns out that dist( B, A 3 ) must be positive due to the design of δ 0 , giving:dist( B, A 3 ) > min (B,Ω)∈E(Θ 0 ) min {(i,j)|Bij ̸ =0} |B ij | -δ 0 =τ - τ 1 + ∆ = ∆ 1 + ∆ τ > 0.By Lemma 8, there exists some α > 0 such that ℓ(B, Ω) -ℓ(B 0 , Ω 0 ) > α for all (B, Ω) ∈ A 3 . Consequently, we have:inf (B,Ω)∈A3 ℓ(B, Ω) -ℓ(B 0 , Ω 0 ) p 1,δ (B 0 ) -p 1,δ (B) > 0.Thus, we can defineλ 0 = inf (B,Ω)∈A3 ℓ(B, Ω) -ℓ(B 0 , Ω 0 ) p 1,δ (B 0 ) -p 1,δ (B) > 0.In summary, for all 0 < λ < λ 0 and 0 < δ < δ 0 , for any ( B, Ω) ∈ E min (Θ 0 ), and for all (B, Ω) / ∈ E min (Θ 0 ), the following holds:ℓ( B, Ω) + p λ,δ ( B) < ℓ(B, Ω) + p λ,δ (B).This concludes the proof.Proof. Here, Θ f (B 0 , Ω 0 ) = Θ 0 . From Theorem 1, we know that when n → ∞ O n,λ,δ = E min (Θ 0 ).Given the additional assumption that p(X) is faithful with respect to G 0 := G(B 0 ), by Lemma 1, we haveM(G 0 ) = E min (Θ 0 ) = {G(B) : (B, Ω) ∈ E min (Θ 0 )}.Note that{G(B) : (B, Ω) ∈ E min (Θ 0 )} = {G(B) : (B, Ω) ∈ O n,λ,δ }.Thus, we conclude thatM(G 0 ) = O n,λ,δ .as n → ∞ This completes the proof.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 Σ based on X as follows:Σ = 1 n [X -1 n • ( µ 1 , . . . , µ p )] ⊤ [X -1 n • ( µ 1 , . . . , µ p )]and plug it into the negative sample log-likelihood function. The same procedure applies to Z:Σ std = 1 n D -1 [X -1 n ( µ 1 , . . . , µ p )] ⊤ [X -1 n ( µ 1 , . . . , µ p )] D -1 = D -1 ΣD -1 .Denote Θ = ( Σ) -1 and Θ std = ( Σ std ) -1 = D ΘD. For Θ, by applying Theorem 1, there exist λ raw,n 0 > 0 and δ raw,n 0 > 0 such that for any 0 < λ < λ raw,n 0 and 0 < δ < δ raw,n 0 , we have O n,λ,δ (X) = E min ( Θ). For D ΘD, we apply Theorem 1 again, and there exist λ std,n 0 > 0 and δ std,n 0 > 0 such that for any 0 < λ < λ std,n 0 and 0 < δ < δ std,n 0 , we have O n,λ,δ (Z) = E min (D ΘD).We can select δ 0 = min{δ raw,n 0 , δ std,n 0 } and λ 0 = min{λ raw,n 0 , λ std,n 0 } in optimization (12) to ensure that O n,λ,δ (X) = E min ( Θ) and O n,λ,δ (Z) = E min (D ΘD) hold simultaneously.By applying Lemma 2, we conclude that:G(O n,λ,δ (X)) = G(E min ( Θ)) = G(E min (D ΘD)) = G(O n,λ,δ (Z)).Furthermore, as n → ∞, we have Σ → Σ and Θ → Θ. Therefore, E min ( Θ) → E min (Θ), andE min ( Θ std ) → E min (Θ std ) as n → ∞. Thus, G(O n,λ,δ (X)) = G(E min (Θ)) = G(O n,λ,δ (Z)).where n → ∞ Note that we use n → ∞ to indicate we consider the result in population level.Proof. This proof shares many similarities with the proof of Theorem 1. First, when n → ∞, we consider the result at the population level. Thus, ℓ n (ψ, ξ) → ℓ(ψ, ξ), and we will focus on ℓ(ψ, ξ) in the following. As a result, we only work with ℓ(ψ, ξ) instead of ℓ n (ψ, ξ).By convention, we can always assume that (ψ 0 , ξ 0 ) ∈ E min (ψ 0 , ξ 0 ). Now, consider the case where p λ,δ (B(ψ 0 )) = 0, which implies that B(ψ 0 ) = 0. Since X ∼ P (X, ψ 0 , ξ 0 ), we have ℓ(ψ 0 , ξ 0 ) ≤ ℓ(ψ, ξ), and thusℓ(ψ 0 , ξ 0 ) + p λ,δ (B(ψ 0 )) ≤ ℓ(ψ, ξ) + p λ,δ (B(ψ)) ∀(ψ, ξ) ∈ Ψ × Ξ, ∀λ > 0, ∀δ > 0.Therefore, (ψ 0 , ξ 0 ) is the optimal solution to optimization (12) .As we iterate, we can assume that |E min (ψ 0 , ξ 0 )| = 1, meaning E min (ψ 0 , ξ 0 ) = {(ψ 0 , ξ 0 )}. This assumption simplifies the proof because any element in E min (ψ 0 , ξ 0 ) is indistinguishable based on the value of ℓ(ψ, ξ) and the penalty for the chosen parameters (δ, λ). Our goal is to find one element by solving the optimization problem (14) , and this assumption simplifies the argument significantly.First, we defineδ 0 = τ 1 + ∆ τ := min (ψ,ξ)∈E(ψ 0 ,ξ 0 ) min {(i,j):B(ψ)ij ̸ =0} |B(ψ)| ij .It is important to note that, under Assumption A (1), since |E(ψ 0 , ξ 0 )| is finite, we have τ > 0.Then, for any λ > 0 and 0 < δ < δ 0 , consider the setA 1 = {(ψ, ξ) | p λ,δ (B(ψ)) = p λ,δ (B(ψ 0 ))}.It is clear that for any (ψ, ξ) ∈ A 1 , we must have (ψ, ξ) / ∈ E(ψ 0 , ξ 0 ), since for any (ψ, ξ) ∈ E(ψ 0 , ξ 0 ), we have p λ,δ (B(ψ)) > p λ,δ (B(ψ 0 )). As a result,ℓ(ψ 0 , ξ 0 ) < ℓ(ψ, ξ) ∀(ψ, ξ) ∈ A 1 .Therefore,ℓ(ψ 0 , ξ 0 ) + p λ,δ (B(ψ 0 )) < ℓ(ψ, ξ) + p λ,δ (B(ψ)) ∀(ψ, ξ) ∈ A 1 .Next, we consider the set A 2 = {(ψ, ξ) | p λ,δ (B(ψ)) > p λ,δ (B(ψ 0 ))}, and we know that ℓ(ψ 0 , ξ 0 ) ≤ ℓ(ψ, ξ). Therefore,ℓ(ψ 0 , ξ 0 ) + p λ,δ (B(ψ 0 )) < ℓ(ψ, ξ) + p λ,δ (B(ψ)) ∀(ψ, ξ) ∈ A 2 .Consequently, we need to check A 3 = {(ψ, ξ) | p λ,δ (B(ψ)) < p λ,δ (B(ψ 0 ))}. For (ψ 0 , ξ 0 ) to be the minimizer of ( 14), we require that the following condition holds:ℓ(ψ 0 , ξ 0 ) + p λ,δ (B(ψ 0 )) < ℓ(ψ, ξ) + p λ,δ (B(ψ)) ∀(ψ, ξ) ∈ A 3 .This condition ensures that the ground truth parameters (ψ 0 , ξ 0 ) correspond to the optimal solution by comparing their score with that of any other parameters in the subset A 3 .It is also worth noting that p λ,δ (t) = λp 1,δ (t) for all t. Therefore, a necessary and sufficient condition for this to hold is:λ < inf (ψ,ξ)∈A3 ℓ(ψ, ξ) -ℓ(ψ 0 , ξ 0 ) p 1,δ (B(ψ 0 )) -p 1,δ (B(ψ)).Note that for (ψ, ξ) ∈ A 3 , we have p 1,δ (B(ψ 0 )) -p 1,δ (B(ψ)) ≤ δ 2 s B(ψ 0 ) . Therefore, the denominator on the RHS cannot be arbitrarily large. Moreover, for any 0 < δ < δ 0 , the following holds:∆ 1 + ∆ τ ≤ ∥B(ψ 0 ) -B(ψ)∥ 2 ≤ L∥ψ 0 -ψ∥ 2 ∀(ψ, ξ) ∈ A 3 .The second inequality follows from Assumption A (b). As a consequence,∥(ψ, ξ) -(ψ 0 , ξ 0 )∥ 2 ≥ ∥ψ 0 -ψ∥ 2 ≥ τ ∆ L(1 + ∆).Thus, we obtainA 3 ⊆ {(ψ, ξ) | ∥(ψ, ξ) -(ψ 0 , ξ 0 )∥ 2 ≥ ∆τ L(1 + ∆) } = A 4 .First, note that A 3 is a nonempty set. Otherwise, the conclusion would hold immediately. Let us select any ( ψ, ξ) ∈ A 3 ̸ = ∅, and defineA 5 = {(ψ, ξ) | ℓ(ψ, ξ) ≤ ℓ( ψ, ξ)}. Then, inf (ψ,ξ)∈A3 ℓ(ψ, ξ) = inf (ψ,ξ)∈A3∩A5 ℓ(ψ, ξ) ≥ inf (ψ,ξ)∈A4∩A5 ℓ(ψ, ξ).It is important to note that A 4 ∩ A 5 is nonempty, since ( ψ, ξ) ∈ A 4 ∩ A 5 . By Assumption B and the properties of ℓ(ψ, ξ), we know that A 5 is a bounded and closed set, and A 4 is a closed set. Consequently, A 4 ∩A 5 is compact. Furthermore, for all (ψ, ξ) ∈ E(ψ 0 , ξ 0 ), we have (ψ, ξ) / ∈ A 4 ∩A 5 . All of this leads to the following conclusion:inf (ψ,ξ)∈A3 ℓ(ψ, ξ) ≥ min (ψ,ξ)∈A4∩A5 ℓ(ψ, ξ) = inf (ψ,ξ)∈A4∩A5 ℓ(ψ, ξ) > ℓ(ψ 0 , ξ 0 ).As a result, we define λ 0 as follows:λ 0 = inf (ψ,ξ)∈A3 ℓ(ψ, ξ) -ℓ(ψ 0 , ξ 0 ) p 1,δ (B(ψ 0 )) -p 1,δ (B(ψ)) > 0.Proof. The proof is combination of Lemma 3 and Theorem 4, similar to Proof of Theorem 2.Before proving the result, we introduce the definitions of the Sparsest Markov representation assumption and restricted faithfulness, along with a few useful theorems. Definition 4 (Sparsest Markov representation [54] ). A pair (G 0 , P ) satisfies the Sparsest Markov Representation (SMR) assumption if (G 0 , P ) satisfies the Markov property and |G| > |G 0 | for every DAG G such that (G, P ) satisfies the Markov property and G ̸ ∈ M(G 0 ).In other words, the SMR assumption asserts that the true DAG G 0 is the (unique up to Markov equivalence) sparsest DAG satisfying the Markov property.Definition 5 (Restricted-faithfulness [53, 54] ). A distribution P satisfies the restricted-faithfulness assumption with respect to a DAG G if it is Markov to G and following two conditions hold:• Adjacency-faithfulness: for all (j, k) ∈ E and all subsets S ⊂ [p]\{j, k} it holds thatX j ̸ ⊥ ⊥ X k | X S• Orientation-faithfulness: for all triples (j, k, l) with skeleton j -l -k and all subsets S ⊂ [p]\{j, k} such that j is d-connected to k given S it holds that X j ̸ ⊥ ⊥ X k | X S Theorem 6 ( [53] ). If a distribution P is faithful to G, then such distribution P also satisfies the restricted-faithfulness assumption with respect to G. Theorem 7 (Theorem 2.4 in [54] ). Let (G, P ) satisfy the Markov property. Then the restrictedfaithfulness assumption implies the SMR assumption.Proof. First, by Theorem 6, the faithfulness assumption implies the restricted faithfulness assumption. Second, by Theorem 7, the restricted faithfulness assumption implies the Sparsest Markov representation assumption. Furthermore, note that for any (B, Ω) ∈ G(E min (Θ 0 )), the distribution P (X) is Markov to G(B), since (B, Ω) ∈ G(E(Θ 0 )). According to the definition of the Sparsest Markov representation assumption, all sparsest DAGs that satisfy the Markov property must belong to the same Markov equivalence class. In our case, this means G(E min (Θ 0 )) = M(G 0 ).Proof. For X ∼ N (0, Σ), let Θ = DΘD, and denote the inverse of Θ as Σ = ( Θ) -1 = D -1 ΣD -1 . It follows that X := D -1 X ∼ N (0, Σ). Now, consider the following least squares regression for j ∈ {1, . . . , p} and S ⊆ {1, . . . , p} \ {j}. Let β ∈ R |S| . Then the following relationships hold:β Sj = arg min β E∥X j -β ⊤ X S ∥ 2 2 ⇒ β Sj = Σ -1 SS Σ Sj βSj = arg min β E∥ Xj -β⊤ XS ∥ 2 2 ⇒ βSj = Σ-1 SS ΣSj Σ-1 SS =([D -1 ΣD -1 ] SS ) -1 = D SS Σ -1 SS D SS ΣSj =[D -1 ΣD -1 ] Sj = D -1 SS Σ Sj D -1 jj βSj = Σ-1 SS ΣSj = D SS Σ -1 SS D SS D -1 SS Σ Sj D -1 jj = D SS β Sj D -1 jjAs a consequence, supp(β Sj ) = supp( βSj ). Note that for all (B, Ω) ∈ E(Θ), we know from Section C.2 that there exists a π ∈ P such that B = B(π). Moreover, B can be recovered by least squares regression using X with its topological sort [2, 13] that is consistent with π. For such a π, we can find a pair ( B, Ω) ∈ E(DΘD), where B has the same topological sort as π, and it can be recovered by least squares regression on X. We have shown that, for the same S, j, supp(β Sj ) = supp( βSj ). Therefore, supp(B) = supp( B).Proof. The proof is the same as Lemma 1.Proof. This follows directly from the definition of the Sparsest Markov Representation (SMR) assumption. Since for all (B, Ω) ∈ E min (Θ), G(B) is Markovian to P and G(B) is the Sparsest, by Definition 2, all G(B) must belong to the same Markov equivalence class by the definition of SMR.Proof. X std = D -1 (X -EX) is based on definition of standardization.Cov(X std ) = Cov(D -1 (X -EX)) = D -1 Cov((X -EX))D -1 = D -1 ΣD -1 [Cov(X std )] -1 = [D -1 ΣD -1 ] -1 = DΣ -1 D = DΘD.Proof. Detailed proof can be found in [41] , Appendix Section D.Proof. From the definition of Σ f (B, Ω):Σ f (B, Ω) := (I -B) -⊤ Ω(I -B) -1 = (I -B) -⊤ Ω 1/2 Ω 1/2 (I -B) -1 ,where Ω 1/2 = diag(ω 1 , . . . , ω p ). It is clear that Σ(B, Ω) is positive semidefinite, asx ⊤ Σ(B, Ω)x = ∥Ω 1/2 (I -B) -1 x∥ 2 2 ≥ 0, ∀x ∈ R p .Next, we just need to show that Ω 1/2 (I -B) -1 x ̸ = 0 for all x ̸ = 0.Ω 1/2 (I -B) -1 x ̸ = 0 ⇔ (I -B) -1 x ̸ = 0 ⇔ x ̸ = 0.Here, ω 2 j > 0 for all j, so Ω 1/2 is invertible. As (I -B) is a full rank matrix, then (I -B) -1 is also a full rank matrix, it indicates that Σ(B, Ω) is positive definite matrix.Since Ω 0 > 0, it follows that Σ 0 is positive definite. By Lemma 6, we have:ℓ(B, Ω) = 1 2 log det Ω -log det(I -B) + 1 2 Tr(Σ 0 Θ(B, Ω)) + const. = 1 2 log det Ω + 1 2 Tr(Σ 0 Θ(B, Ω)) + const. ≥ℓ(B, Ω f (B)) = ℓ(B)where Ω f (B) and ℓ(B) are defined in Equations ( 15) and ( 17), respectively. The last inequality follows from Section C.1. Next, we need to prove that ℓ(B) > -∞.ℓ(B) = 1 2 log det diag((I -B) ⊤ Σ 0 (I -B)) + const. = 1 2 p j=1 log E∥X j -B ⊤ j X∥ 2 2 + const. = 1 2 p j=1 log E∥(e j -B j ) ⊤ X∥ 2 2 + const. = 1 2 p j=1log(e j -B j ) ⊤ Σ 0 (e j -B j ) + const.≥ 1 2 p j=1 log ∥(e j -B j )∥ 2 2 Λ min (Σ 0 ) + const. ≥ 1 2 p j=1 log Λ min (Σ 0 ) > -∞Here, e j ∈ R p is a unit vector with the j-th position equal to 1 and all other positions being zero, and Λ min (Σ 0 ) is the minimum eigenvalue of Σ 0 . Since Σ 0 is positive definite, we have Λ min (Σ 0 ) > 0.Because B is the adjacency matrix of a DAG, it follows that B jj = 0, which implies ∥e j -B j ∥ ≥ ∥e j ∥ = 1. As a result,ℓ(B, Ω) ≥ ℓ(B) > -∞.Proof. Note that for a fixed B, the corresponding optimalΩ f (B) = diag((I -B) ⊤ Σ 0 (I -B))is the solution with respect to ℓ(B, Ω). Therefore, without causing confusion, we take Ω f (B) = diag((I -B) ⊤ Σ 0 (I B)) and consider the log-likelihood as a function of B only, i.e., ℓ(B), for simpler representation. See Equation (17) in Section C.1 for details. It is clear that 0 ∈ A 3 , so we defineA 4 = {B | ℓ(B) ≤ ℓ(0)}. Note that ℓ(0) is finite. ℓ(0) ≥ ℓ(B) = p j=1log(e j -B j ) ⊤ Σ 0 (e j -B j )≥ p j=1 log ∥e j -B j ∥ 2 Λ min (Σ 0 ) = log   (Λ min (Σ 0 )) p p j=1 ∥e j -B j ∥ 2   .This indicates thatp j=1 ∥e j -B j ∥ 2 ≤ exp(ℓ(0)) (Λ min (Σ 0 )) p Moreover, ∥e k -B k ∥ 2 ≤ p j=1 ∥e j -B j ∥ 2 ≤ exp(ℓ(0)) Λ p min ∀k ∈ {1, . . . , p}This implies that B k must be bounded, and therefore everyB in A 4 is bounded. It is clear that arg min B∈A3 ℓ(B) ∈ A 4 . Thus, we need to show that min B∈A3 ℓ(B) = min B∈A3∩A4 ℓ(B) > ℓ(B 0 ). Define A 5 = { B | dist( B, E(Θ 0 )) ≥ 1 2 min B∈E(Θ 0 ) dist(B, A 3 )}It is easy to see that A 3 ⊆ A 5 . Then,min B∈A3∩A4 ℓ(B) ≥ min B∈A4∩A5 ℓ(B).Note that A 4 ∩ A 5 is closed, bounded, and nonempty (0 ∈ A 4 ∩ A 5 ), and ℓ(B) is a continuous function of B. Consequently, there exists at least one minimizer of ℓ(B) in A 4 ∩ A 5 . Combining this with the fact that B / ∈ A 4 ∩ A 5 for all B ∈ E(Θ 0 ), we conclude that:min B∈A3 ℓ(B) ≥ min B∈A4∩A5 ℓ(B) > ℓ(B 0 ).Proof. By Theorem 1, we know there exists λ 0 > 0 and δ 0 > 0. For MCP, it can be transformed into quasi-MCP, by reparameterization from Section C.5. Then, combining these results together.0 < λ = λ mcp < λ 0 , 0 < δ = aλ mcp < δ 0 ⇒ 0 < λ mcp < λ 0 , 0 < a < δ 0 λ mcpWe could simple set a 0 : δ0 λ0 and (λ mcp ) 0 = λ 0 . For SCAD, we just requires the following is satisfied to satisfies the pattern in the proof of Theorem 1.0 < aλ scad < δ 0 , 0 < λ 2 scad (a + 1) 2 < λ 0One simple choice is to leta 0 = (λ scad ) 0 < min{ δ 0 , λ 0 , 1}This completes the proof of Corollary 1.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-likelihood function for the model in Equation ( 6) (Appendix C.1).• A brief introduction to the characterization of the equivalence class E(Θ) (Appendix C.2).• Examples demonstrating that the optimal solution for the least squares loss differs from the optimal solution of the log-likelihood (Appendix C.3). • An example illustrating the estimation bias when the ℓ 1 penalty is applied (Appendix C.4).• The formulations for quasi-MCP, MCP, and SCAD (Appendix C.5).• The standardization of the random variable X and the dataset X (Appendix C.6).• A detailed discussion of Assumptions A and B (Appendix C.7).C.1 Log-likelihood of Model (6) In this subsection, we detail the negative log-likelihood of the model in Equation (6) . ℓ n (B, Ω) = - 1 n log n i=1 f (x i ; B, Ω) = - 1 n log n i=1 1 (2π) p/2 (det Σ(B, Ω)) 1/2 exp - x ⊤ i Θ(B, Ω)x i 2 = p 2 log 2π + 1 2 log det Σ(B, Ω) + 1 2n n i=1 x ⊤ i Θ(B, Ω)x i = 1 2 log det(I -B) -⊤ Ω(I -B) -1 + 1 2 Tr(Θ(B, Ω)( n i=1 x ⊤ i x i n )) + const.((I -B) ⊤ Σ(I -B)) jj ω 2 j -log det(I -B) + const.It is easy to know that for any fix B, the optimal solution of (ω * j ) 2 can be written as:(ω * j ) 2 = [(I -B) ⊤ Σ(I -B)] jj Therefore, optimal solution Ω f (B) for any fixed B can be written as:EQUATIONLet us define profile sample log-likelihood ℓ n (B) as function of B with such optimal Ω(B) plugged inEQUATIONWhere X = (X 1 , . . . , X p ) = (x 1 , . . . , x n ) ⊤ and corresponding profile population log-likelihoodEQUATIONWe provide a brief introduction to the equivalence class E(Θ), 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 is an ordering on the nodes, often denoted by ≺, such that the existence of a directed edge X k → X j implies that X k ≺ X j in the ordering.Let P denote the collection of all permutations of the indices {1, . . . , p}. For an arbitrary matrix A and any π ∈ P, let P π A represent the matrix obtained by permuting the rows and columns of A according to π, such that (P π A) ij = a π(i)π(j) .A DAG B is said to be compatible with permutation π if P π B is a lower-triangular matrix, which is equivalent to saying thatX k → X j in B implies that π -1 (k) > π -1 (j).Similarly, π is also called compatible with B.For any positive definite matrix Θ and π ∈ P, the matrix P π Θ represents the same covariance structure as Θ, up to a reordering of the variables. The Cholesky decomposition of P π (Θ) can be uniquely written as:P π Θ = (I -L)D -1 (I -L) ⊤ = Θ f (L, D),where L is strictly lower triangular and D is diagonal. By Lemma 8 in [2] , the following holds:P π Θ(L, D) = Θ(P π L, P π D) ∀π ∈ P.Therefore,Θ = Θ f (P π -1 L, P π -1 D).For each π, we define:B(π) := P π -1 L, Ω(π) := P π -1 D.This suggests that for any π ∈ P, there exists a pair ( B(π), Ω(π)) ∈ E(Θ), where B(π) can be uniquely determined based on the permutation π and Θ [2] . It is important to emphasize that different permutations, π 1 ̸ = π 2 , can still result in the same pairs, i.e., ( B(π 1 ), Ω(π 1 )) = ( B(π 2 ), Ω(π 2 )). Furthermore, this indicates that for any (B, Ω) ∈ E(Θ), there exists at least one permutation π such that (B, Ω) = ( B(π), Ω(π)). Moreover, it turns out that the collection of pair of ( B(π), Ω(π)) forms the entire equivalence class E(Θ). Lemma 9 (Lemma 1, 2). Suppose Σ is a positive definite covariance matrix and Θ = Σ -1 . Then,E(Θ) ={(P π -1 L, P π -1 D) : P π Θ = Θ f (L, D), π ∈ P} ={( B(π), Ω(π)) : π ∈ P}This result indicates that the size of E(Θ) is at most p!, which is large but finite.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 is just Example 1 in [33] . Suppose (X 1 , X 2 ) is distributed according to the following linear SEM with unequal variances:X 1 = ϵ 1 , X 2 = - X 1 2 + ϵ 2 , ϵ 1 ∼ N (0, 1), ϵ 2 ∼ N (0, 1/4).ThusB 0 = 0 -1/2 0 0 , Ω 0 = 1 0 0 1/4 , Σ 0 = Σ f (B 0 , Ω 0 ) = 1 -1/2 -1/2 1/2and alsoE(Θ 0 ) = E min (Θ 0 ) = {B 0 , B 1 } where B 1 = 0 0 -1 0 , Ω 1 = 1 0 0 1 .Moreover, ℓ(B 0 , Ω 0 ) = ℓ(B 1 , Ω 1 ), since both SEM represent the same covariance.But it turns out that Similar calculations can be carried out for d > 2, but are tedious owing to the size of E(Θ). For example, here is an example of an 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:E[∥X -B ⊤ 1 X∥ 2 ] < E[∥X -(B 0 ) ⊤ X∥ 2 ]: More precisely, it is easy to check that E[∥X -B ⊤ 1 X∥ 2 ] = Tr((I -B 1 ) ⊤ Σ 0 (I -B 1 )) = 1, E[∥X -(B 0 ) ⊤ X∥ 2 ] = Tr((I -B 0 ) ⊤ Σ 0 (I -B 0 )) = 5/4,B 0 = 0 0 -3/10 0 0 -2 0 0 0 , Ω 0 = 7 0 0 0 3 0 0 0 2 , Σ 0 = Σ f (B 0 , Ω 0 ) = 7 0 -2 0 3 -5 -2 -5 10.For this model, LS loss selects the following SEM with 3 edges:B 1 = 0 0 0 -1.197 0 -1.589 -0.7532 0 0 . We have B 0 ∈ E min (Θ 0 ), but B 1 / ∈ E min (Θ 0 ).We provide an example showing that when the ℓ 1 penalty is applied, the estimation becomes biased. Therefore, ℓ 1 should not be used. Consider the following linear Structural Equation Model (SEM):X 1 = N (0, 1) X 2 = X 1 + N (0, σ 2 )If the topological sort is known, i.e., X 1 → X 2 , and an ℓ 1 penalty is used for minimizing the negative log-likelihood:min a log E[∥X 2 -aX 1 ∥ 2 2 ] + log E[∥X 1 ∥ 2 2 ] + λ|a|Ideally, we would expect a = 1 to be the minimal solution to the loss function. However, the problem is equivalent to:log((1 -a) 2 + σ 2 ) + λ|a|It is clear that a = 1 is not the minimal solution to the loss function, as the derivative at a = 1 is nonzero for any λ > 0, indicating that the ℓ 1 penalty leads to a biased estimator. This bias does not occur when using MCP or SCAD with appropriate hyperparameters.We present the formulas for quasi-MCP, MCP, and SCAD, and demonstrate that quasi-MCP and MCP are equivalent.p λ,δ (t) = λ |t| - t 2 2δ 1(|t| < δ) + δ 2 1(|t| > δ) MCP p M CP λ,a (t) = 1(|t| < aλ) λ|t| - t 2 2a + 1(|t| ≥ λa) λ 2 a 2 SCAD p SCAD λ,a (t) = λ|t|1(|t| < λ) + 1(λ < |t| < aλ) 2aλ|t| -t 2 -λ 2 2(a -1) + 1(|t| ≥ λa)λ 2 (a + 1) 2 We present the formulas for the standardization of X and the standardization of the corresponding dataset X.i := Var(X i ) and D := diag(σ 1 , . . . , σ p ). Denote the standardized version of X as X std , which can be expressed as:X std = D -1 (X -EX).For X ∈ R n×p , we can write X = (X ij ) = (x 1 , . . . , x n ) ⊤ = (X 1 , . . . , X p ), and define sample average for node j as µ jµ j = 1 n n i=1 X ij ∀j ∈ [p].Next, we define the sample variance for node j asσ 2 j = 1 n -1 n i=1 (X ij -µ j ) 2 .The diagonal matrix of sample standard deviations is thenD = diag( σ 1 , . . . , σ p ).Finally, we standardize X by subtracting the sample means and scaling by the inverse of D:Z = [X -1 n • ( µ 1 , . . . , µ p )] D -1 ,where 1 n ∈ R n is an n-dimensional vector with all entries equal to 1.In this subsection, we provide a more detailed discussion of Assumptions A and B.First, it is important to emphasize that if p λ,δ is replaced by the ℓ 0 penalty in (12) or (14) , then Assumptions A and B 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 ℓ 0 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 is truly needed is that the minimal nonzero edge has sufficient "signal," i.e., min(ψ,ξ)∈E(ψ 0 ,ξ 0 ) min {(i,j):B(ψ)ij ̸ =0} |B(ψ)| ij > 0 This is trivially true when |E(ψ 0 , ξ 0 )| is finite. When |E(ψ 0 , ξ 0 )| is infinite, each |B(ψ)| ij could be positive, but it is possible lim inf (ψ,ξ)∈E(ψ 0 ,ξ 0 ) min {(i,j):B(ψ)ij ̸ =0} |B(ψ)| ij = 0, because |B(ψ)| ijcan be arbitrarily small. The ℓ 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 Assumption B is a 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 B(ψ). The magnitude of the quasi-MCP penalty does not directly reveal the number of edges. This is the trade-off for replacing the ℓ 0 penalty with a 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 ϵ > 0, there exists δ > 0 such that ℓ(ψ, ξ) -ℓ(ψ 0 , ξ 0 ) > δ for all {(ψ, ξ) | dist((ψ, ξ), E(ψ 0 , ξ 0 )) > ϵ}.In other words, we require a loss gap when (ψ, ξ) is not in E(ψ 0 , ξ 0 ). This can be inferred from Assumption B.The following is the proof when p λ,δ is replaced with the ℓ 0 penalty.Proof when p λ,δ is replaced with ℓ 0 :Proof. We can also assume that |E min (ψ 0 , ξ 0 )| = 1, meaning E min (ψ 0 , ξ 0 ) = {(ψ 0 , ξ 0 )}. In other words, there is a unique element in the minimal equivalence class. This is because any element in E min (ψ 0 , ξ 0 ) is indistinguishable based on the score function, i.e., the value of ℓ(ψ, ξ) and the penalty for the number of edges in B(ψ). Our objective is to find this unique element by solving equation (12) in the paper, which simplifies the proof.When s B(ψ 0 ) = 0, the result is straightforward:ℓ(ψ 0 , ξ 0 ) + s B(ψ 0 ) ≤ ℓ(ψ, ξ) + s B(ψ) .Now, let us consider the more general case where s B(ψ 0 ) > 0 and divide the parameter space into three regions:A 1 = {(ψ, ξ) | s B(ψ) > s B(ψ 0 ) }, A 2 = {(ψ, ξ) | s B(ψ) = s B(ψ 0 ) }, A 3 = {(ψ, ξ) | s B(ψ) < s B(ψ 0 ) }.Case 1: Consider A 1 . Since ℓ(ψ 0 , ξ 0 ) ≤ ℓ(ψ, ξ), the following holds for any λ > 0:ℓ(ψ 0 , ξ 0 ) + λs B(ψ 0 ) < ℓ(ψ, ξ) + λs B(ψ) ∀(ψ, ξ) ∈ A 1 .Case 2: Consider A 2 . Since |E min (ψ 0 , ξ 0 )| = 1, it follows that for all (ψ, ξ) ∈ E(ψ 0 , ξ 0 ) and (ψ, ξ) ̸ = (ψ 0 , ξ 0 ), we have s B(ψ) > s B(ψ 0 ) . Therefore, for all (ψ, ξ) ∈ A 2 , it holds that ℓ(ψ 0 , ξ 0 ) < ℓ(ψ, ξ). Consequently, for any λ > 0:ℓ(ψ 0 , ξ 0 ) + λs B(ψ 0 ) < ℓ(ψ, ξ) + λs B(ψ) ∀(ψ, ξ) ∈ A 2 .Case 3: Consider A 3 . We need to prove that:ℓ(ψ 0 , ξ 0 ) + λs B(ψ 0 ) < ℓ(ψ, ξ) + λs B(ψ) ∀(ψ, ξ) ∈ A 3 .• LOGLL(-NOTEARS/DAGMA)-SAMPLE/POPULATION 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 ℓ 1 penalty, it employs a log-likelihood function that includes a quasi-MCP penalty, as defined in (10) . For LOGLL(-NOTEARS/DAGMA)-SAMPLE, we use the sample covariance matrix Σ in the score function. In contrast, LOGLL(-NOTEARS/DAGMA)-POPULATION uses the true covariance matrix Σ as a baseline approach.In this paper, LOGLL-NOTEARS refers to solving (12) using NOTEARS with NLL and quasi-MCP as the score function, while LOGLL-DAGMA 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 B(π) as 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 an average basis. This involves identifying the set M min (Θ) or M min ( Θ), calculating the SHD for each DAG within this set, and then computing the average SHD.Linear model There are two main challenges in 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/DAGMA often gets stuck at a local optimal solution. The second challenge arises from Theorem 1, where we are advised to select λ and δ such that 0 < λ < λ 0 and 0 < δ < δ 0 . Theoretically, smaller values for λ and δ should be used to adhere to the theorem'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 ℓ 1 penalty) or DAGMA (with least squares loss and ℓ 1 penalty) to obtain a "good" initialization point. Then, we apply LOGLL(-NOTEARS/DAGMA)-POPULATION/SAMPLE to obtain the final output.To address second challenge, we use warm starts. We begin with larger values for λ and δ and solve (12) using LOGLL-NOTEARS/DAGMA to obtain an initial B est . We then reduce λ and δ by a 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 ℓ n (B est , Ω est ) 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 1.In Algorithm 1, we detail the complete implementation of LOGLL(-NOTEARS/DAGMA)-SAMPLE. By replacing Σ with Σ and substituting ℓ n with ℓ, this algorithm is adapted to the full implementation of LOGLL(-NOTEARS/DAGMA)-POPULATION.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 hyperparameter unchanged. The score function (negative log-likelihood) we uses N LL (f, X) = 1 2n d i=1 log ∥x i -f i (X)∥ 2Here fi 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 is scale-invariant for the linear model with Gaussian noise (see Theorem 3), it was observed that using standardized • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a 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 is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are 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 to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).8: Graph: structure X 0 → X 1 , X 0 → X 2 , X 1 → X 3 , X 2 → X 3 . For 0 < δ < δ 0 , the estimated (B est , Ω est ) ∈ E min (Θ 0 ) because SHDGiven that ψ describes all the parameters for the functions f , we will also use B(ψ) to denote B(f ) in (3).In terms of notation given in Section 3, we have parameters ψ = B and ξ = (ω 2 1 , . . . , ω 2 p ), and parameter spaces Ψ = R p×p , Ξ = R p >0 .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 error variance Ω.Although ℓ1 is nonsmooth, standard smoothing techniques can be applied to ℓ1 regularizers as in[73,4].See Appendix C.4 for examples. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

- 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.__
- 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?__
- 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?__. 

[1] Notes by Charles Geyer: https://www.stat.umn.edu/geyer/8053/notes/expfam.pdf

[2] "Learning Sparse Nonparametric DAGs", Zheng et al., AISTATS 2020.
questions: questions: ### Other questions  

- 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$?
- 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: 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. 

I 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.__
paper_summary: paper_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).

**Review summary:**
This 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
- rewriting the paper in a more transparent way, especially by contrasting with the contributions of [b] (and changing the title); 
- highlighting what are the actual contributions (scale-invariant result + MEC identifiability via quasi-MCP regularizer); and 
- providing convincing experiments with nonlinear models.
strength_weakness: strength_weakness: **Important prior work omitted which seriously limits novelty**

The 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...

I 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).

**Experiments investigate only the linear case**

The 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. 

**Limitations of the theory:**

Assumption 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.

Assumption 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.

**Unclear limit statements in Theorem 1 to 4**

The 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)

Minor:
- It would be nice to have a small plot of the quasi-MCP regularizer.
- Typo in Lemma 3 in appendix? I think an \mathcal{M} should be replaced by a \mathcal{G}, no?
- Text is absurdly small in Figure 1 e.g.. Please fix this.
- Is the sparsest Markov representation always unique? The phrasing used suggests it is, but it might be worth expliciting.

[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.

[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.

[c] S. Lachapelle, P. Brouillard, T. Deleu, and S. Lacoste-Julien. Gradient-based neural DAG
learning. In Proceedings of the 8th International Conference on Learning Representations,
2020.
questions: questions: See above.
limitations: limitations: This paper lacks a conclusion/discussion, so the limitations of the approach are not addressed.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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].

When 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.

Real-world datasets can be added to validate the proposed method.

Code is not available.

Some minor flaws in the presentation:
+ The conclusion section is missing.
+ Section D.3.2, time spent on the simulation. The plots are about SHD rather than about the running time.
+  [72] in the reference should be cited correctly.
questions: 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.

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.

It would be interesting to compare the proposed method with the existing ones on real datasets.
limitations: limitations: The method assumes a known likelihood function, which can be impractical when the noise distribution or the functional structure is unknown.

Comments: title: Thanks!
Thank you so much!!!
Thank you!
Thank you!
Thank you for valuable comment!
comment: Dear Reviewer,

We sincerely appreciate your positive feedback on our rebuttal. We will incorporate [1] and related works into the paper and expand the related work section to discuss the similarities and differences between our work and [1].

We acknowledge that the current title is too general, and we will refine it to better reflect the paper’s content. Additionally, all statements involving $n\rightarrow \infty$ will be revised to refer to "population" results.

Thank you once again for your invaluable suggestions and the time you’ve dedicated to reviewing our work.
I thank the authors for engaging with my review and providing a convincing rebuttal. The additional discussion contrasting their theoretical result with that of [1] was enlightening and very clear. I urge the authors to integrate these discussions in the next version of their work to better contextualize their contribution within existing literature. The authors should also consider changing the title to something less general, one way to do that would be to mention the quasi-MCP regularizer somehow (this paper is not the first differentiable likelihood based structure learning approach proposed!). Also, please remove the limit statement as these are inaccurate and should be replaced by "population results" (\ell instead of \ell_n). 

Assuming these changes will be properly addressed in the next revision, I raise my score from 3 to 5. I do not raise further as I still believe the novelty factor fairly low.
Dear authors,

Thank you very much for your efforts in the discussion. Explicit examples and mathematical justification have shown that the model assumptions can incorporate general cases. Therefore, I will update my score to 6.
We sincerely appreciate the time and effort you put into reviewing our paper and providing valuable feedback. We will incorporate your suggestions to enhance the clarity and strength of our work for our readers!
Dear Authors, thank you for the clarification. This addresses my main concern, so I will raise the score to 7 assuming it is clarified in the paper.
Thank you for the reviewer’s comment. We hope you can have a more comprehensive assessment after discussing with the other reviewers. We appreciate your time and effort.
Thank you for the reviewer’s further comments, and we apologize for the unclear part. We are happy to address any concerns you have and ensure a clear understanding of our work. Please rest assured, we will incorporate all the suggestions into the paper, including moving Section C.2 to the main paper to provide readers with more background knowledge.



> Is it the case that for a given sparsity pattern (i.e., DAG), only one $B$ setting in the Gaussian example will yield the same distribution?
> 

You are correct! From my understanding, if a sparsity pattern here refers to a topological sort $\pi$. For each $\pi$, there will be **only one** $B$ that is consistent with $\pi$ and can generate the distribution.

Consider a simple three-node example. Let $X \sim \mathcal{N}(0, \Sigma)$, where $\Sigma \succ 0$:
$$
\Sigma = \begin{bmatrix}
1/2&1&1\\\\1&3&3\\\\1&3&4
\end{bmatrix}
$$
We would like to recover $B$ (DAG) such that $X = B^\top X + N$ and $X \sim \mathcal{N}(0, \Sigma)$. In this case, there are $p! = 3! = 6$ different topological sorts, and at most six different $B$s can generate $X$. Based on Section C2, **it is worth noting that aside from these different $B$s, no other $B$ can generate the distribution $X$ [1]**. These $B$s and $\pi$s are:
$$
B(\pi_1) = \begin{bmatrix}
0&0&0\\\\1/3&0&0\\\\0&3/4&0
\end{bmatrix}\qquad \pi_1 = [2,1,0]
$$
$$
B(\pi_2) = \begin{bmatrix}
0&0&0\\\\1/3&0&1\\\\0&0&0
\end{bmatrix}\qquad \pi_2 = [1,2,0]
$$
$$
B(\pi_3) = \begin{bmatrix}
0&1&0\\\\0&0&0\\\\1/4&1/2&0
\end{bmatrix}\qquad \pi_3 = [2,0,1]
$$
$$
B(\pi_4) = \begin{bmatrix}
0&1&2\\\\0&0&0\\\\0&1/2&0
\end{bmatrix}\qquad \pi_4 = [0,2,1]
$$
$$
B(\pi_5) = \begin{bmatrix}
0&0&0\\\\1/3&0&1\\\\0&0&0
\end{bmatrix}\qquad \pi_5 = [1,0,2]
$$
$$
B(\pi_6) = \begin{bmatrix}
0&2&0\\\\0&0&1\\\\0&0&0
\end{bmatrix}\qquad \pi_6 = [0,1,2]
$$
Here $B(\pi_5) = B(\pi_2)$, although $\pi_2\ne \pi_5$. This is why we say it at most $p!$ different $B$. 

For unidentifiable exponential families, consider the case where $X \sim \mathcal{N}(0,\Sigma)$ and $\Sigma$ is singular. For example,

$$
\Sigma = \begin{bmatrix}
1&1&1\\\\1&1&1\\\\1&1&2
\end{bmatrix}
$$

We can reduce this to a lower dimension. Note that here $X_1 = X_2$, so we only need to consider $(X_2, X_3)$:

$$
\text{Cov}(X_2,X_3) = \begin{bmatrix}
1&1\\\\1&2
\end{bmatrix}\succ0
$$

Then it reduces to previous case.

As for more general model, for each $\pi$, it is usually associated with **only one** pair $(\psi(\pi), \xi(\pi))$ in a certain way such that $P(x; \psi(\pi), \xi(\pi)) = P(x; \psi^0, \xi^0)$, similar to Linear Guassian case. Another example in **Common concern (1)** follows such property.

> I'm wondering if you mean C.2? I can't actually find this statement anywhere.
> 

Apologies for the misunderstanding. You are correct; it should be Section C.2. We will include more details in Section 4.2 for clarity!

As stated in Section C.2, for any topological sort $\pi$, it corresponds to a pair $(\tilde{B}(\pi), \tilde{\Omega}(\pi)) \in \mathcal{E}(\Theta)$. For all DAGs with $p$ nodes, there are at most $p!$ different topological sorts. This explains why there are at most $p!$ different $B$ matrices that satisfy this condition.

[1] Aragam, Bryon, and Qing Zhou. "Concave penalized estimation of sparse Gaussian Bayesian networks." *The Journal of Machine Learning Research* 16.1 (2015): 2273-2328.
Thank you for the response. I will retain my score for now, to discuss with other reviewers, who also found the equivalence classes confusing. Based on your $p!$ characterization, it seems equivalence classes contain either DAGs, or parameter settings in 1-to-1 correspondence with DAGs. Upon reading C.2., is it the case that for a given sparsity pattern (i.e., DAG), only one B setting in the Gaussian example will yield the same distribution? I was originally thinking that rotations gave infinitely many equivalent parameters (and this generalizes in EFs) but this would break the sparsity. Could the authors comment on whether this understanding is correct?

Note this was not at all clear to me on the first few passes so I would appreciate if you could put a version of C.2 in the main text. 

> As we discuss in Section 4.2, there are at most $p!$ different $B$ that satisfy this condition

I'm wondering if you mean C.2? I can't actually find this statement anywhere.
Thank you for your response. Most of my concerns are addressed. At this moment, I would like to retain my score before further discussions with other reviewers.