Score matching through the roof: linear, nonlinear, and latent variables causal discovery

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 12:46:38

Abstract: Causal discovery from observational data holds great promise, but existing methods rely on strong assumptions about the underlying causal structure, often requiring full observability of all relevant variables. We tackle these challenges by leveraging the score function $\nabla \log p(X)$ of observed variables for causal discovery and propose the following contributions. First, we generalize the existing results of identifiability with the score to additive noise models with minimal requirements on the causal mechanisms. Second, we establish conditions for inferring causal relations from the score even in the presence of hidden variables; this result is two-faced: we demonstrate the score's potential as an alternative to conditional independence tests to infer the equivalence class of causal graphs with hidden variables, and we provide the necessary conditions for identifying direct causes in latent variable models. Building on these insights, we propose a flexible algorithm for causal discovery across linear, nonlinear, and latent variable models, which we empirically validate.

Main paper: 
 introduction: 
 The inference of causal effects from observations holds the potential for great impact arguably in any domain of science, where it is crucial to be able to answer interventional and counterfactual queries from observational data [1, 2, 3] . Existing causal discovery methods can be categorized based on the information they can extract from the data [4] , and the assumptions they rely on. Traditional causal discovery methods (e.g. PC, GES [5, 6] ) are general in their applicability but limited to the inference of an equivalence class. Additional assumptions on the structural equations generating effects from the cause are, in fact, imposed to ensure the identifiability of a causal order [7, 8, 9, 10] .As a consequence, existing methods for causal discovery require specialized and often untestable assumptions, preventing their application to real-world scenarios.Further, the majority of existing approaches are hindered by the assumption that all relevant causes of the measured data are observed, which is necessary to interpret associations in the data as causal relationships. Despite the convenience of this hypothesis, it is often not met in practice, and the solutions relaxing this requirement face substantial limitations. The FCI algorithm [11] can only return an equivalence class from the data. Appealing to additional restrictions ensures the identifiability of some direct causal effects in the presence of latent variables: RCD [12] relies on the linear non-Gaussian additive noise model, whereas CAM-UV [13] requires nonlinear additive mechanisms. Nevertheless, the strict conditions on the structural equations hold back their applicability to more general settings.Our paper tackles these challenges and can be put in the context of a recent line of work that derives a connection between the score function ∇ log p(X) and the causal graph underlying the data-generating process [14, 15, 16, 17, 18, 19] . The use of the score for causal discovery is practically appealing, as it yields advantages in terms of scalability to high dimensional graphs [16] Submitted to 38th Conference on Neural Information Processing Systems (NeurIPS 2024). Do not distribute. and guarantees of finite sample complexity bounds [20] . Instead of imposing assumptions that ensure strong, though often impractical, theoretical guarantees, we organically demonstrate different levels of identifiability based on the strength of the modeling hypotheses, always relying on the score function to encode all the causal information in the data. Starting from results of Spantini et al. [21] and Lin [22] , we show how constraints on the Jacobian of the score ∇ 2 log p(X) can be used as an alternative to conditional independence testing to identify the Markov equivalence class of causal models with hidden variables. Further, we prove that the score function identifies the causal direction of additive noise models, with minimal assumptions on the causal mechanisms. This extends the previous findings of Montagna et al. [17] , limited by the assumption of nonlinearity of the causal effects, and Ghoshal and Honorio [14] , limited to linear mechanisms. On these results, we build the main contributions of our work, enabling the identification of direct causal effects in hidden variables models.Our main contributions are as follows: (i) We present the necessary conditions for the identifiability of direct causal effects and the presence of hidden variables with the score in the case of latent variables models. (ii) We propose AdaScore (Adaptive Score-based causal discovery), a flexible algorithm for causal discovery based on score matching estimation of ∇ log p(X) [23] . Based on the user's belief about the plausibility of several modeling assumptions on the data, AdaScore can output a Markov equivalence class, a directed acyclic graph, or a mixed graph, accounting for the presence of unobserved variables. To the best of our knowledge, the broad class of causal models handled by our method is unmatched by other approaches in the literature. 
model definition and related works: 
 In this section, we introduce the formalism of structural causal models (SCMs), separately for the the cases with and without hidden variables.Let X be a set of random variables in R defined according to the set of structural equationsEQUATIONN i ∈ R are mutually independent random variables with strictly positive density, known as noise or error terms. The function f i is the causal mechanism mapping the set of direct causes X PA G i of X i and the noise term N i , to X i 's value. A structural causal model (SCM) is defined as the tuple (X, N, F, P N ), where F = (f i ) k i=1 is the set of causal mechanisms, and P N is the joint distribution relative to the density p N over the noise terms N ∈ R k . We define the causal graph G as a directed acyclic graph (DAG) with nodes X = {X 1 , . . . , X k }, and the set of edges defined as Under this model, the probability density of X satisfies the Markov factorization (e.g. Peters et al.{X j → X i : X j ∈ X PA G i },[1] Proposition 6.31):EQUATIONwhere we adopt the convention of lowercase letters referring to realized random variables, and use p to denote the density of different random objects, when the distinction is clear from the argument. This factorization is equivalent to the global Markov condition (e.g. Peters et al. [1] Proposition 6.22) that demands that for all {X i , X j } ∈ X, X Z ⊆ X \ {X i , X j }, thenX i |= d G X j |X Z =⇒ X i |= X j |X Z ,where (• |= • |•) denotes probabilistic conditional independence of X i , X j given X Z , and(• |= d G • |•)is the notation for d-separation, a criterion of conditional independence defined on the graph G (Definition 5 of the appendix). As it is commonly done, we assume that the reverse directionX i |= X j |X Z =⇒ X i |= d G X j |X Z hold, and we say that the density p is faithful to the graph G [2, 24] (hence the faithfulness assumption). Together with the global Markov condition, faithfulness implies an equivalence between the probabilistic and graphical notions of conditional independence:EQUATIONIn general, several DAGs may entail the same set of d-separations: graphs sharing such common structure form a Markov equivalence class (see Definition 6 in the appendix).The above model assumes that there aren't any unobserved causes of variables in X, other than the noise terms in N . As we are interested in distributions with potential hidden variables, we will now generalize our model to represent data-generating processes that may involve latent causes.As graphs play a central role in our work, Appendix A.1 provides a detailed overview of the fundamental notation and definitions that we rely on in the remainder of the paper. For the next section, we advise the reader to be comfortable with the notions of ancestors (Definition 2) and inducing paths (Definition 3) in DAGs.Closely related works. Several methods for the causal discovery of fully observable models using the score have been recently proposed. Ghoshal and Honorio [14] demonstrates the identifiability of the linear non-Gaussian model from the score, and it is complemented by Rolland et al. [15] , whichshows the connection between score matching estimation of ∇ log p(X) and the inference of causal graphs underlying nonlinear additive noise models with Gaussian noise terms, also allowing for sample complexity bounds [20] . Montagna et al. [17] provides identifiability results in the nonlinear setting, without posing any restriction on the distribution of the noise terms. Montagna et al. [16] is the first to show that the Jacobian of the score provides information equivalent to conditional independence testing in the context of causal discovery, limited to the case of additive noise models.All of these studies make specialized assumptions to find theoretical guarantees of identifiability, whereas our paper provides a unifying view of causal discovery with the score function, which generalizes and expands the existing results.Under the model (1), we consider the case where the set of variables X is partitioned into the disjoint subsets of observed random variables V = {V 1 , . . . , V d } and unobserved (or latent) random variablesU = {U 1 , . . . , U p }.We assume that the following set of structural equations is satisfied:EQUATIONwhere U i stands for the set of unobserved parents of V i , andV PA G i = {V k |k ∈ PA G i , V k ∈ V } arethe observed direct causes of V i . Some of the causal relations and the conditional independencies implied by the set of equations ( 4) can be summarized in a graph obtained as a marginalization of the DAG G onto the observable nodes V .Definition 1 (Marginal graph, Zhang [25] ). Let X = V ∪U and G be a DAG over X. The following construction gives the marginal graph M G V , with nodes V and edges found as follows:• pair of nodes V i , V j are adjacent in the graph M G V if and only if there is an inducing path between them relative to U in G;• for each pair of adjacent nodesV i , V j in M G V , orient the edge as V i → V j if V i is an ancestor of V j in G, else orient it as V i ↔ V j .We define the map G → M G V as the marginalization of the DAG G onto V , the observable nodes.The graph resulting from the above construction is a maximal ancestral graph (MAG, Definition 4), hence we will often refer to it as the marginal MAG of G. Intuitively, a directed edge denotes the presence of an ancestorship relation, whereas bidirected edges represent dependencies that can not be removed by conditioning on any of the variables in the graph.In the case of DAGs, d-separation encodes the probabilistic conditional independence relations between the variables of X in the graph G, as explicit by Equation (3). Such notion of graphical separation has a natural generalization to maximal ancestral graphs, known as m-separation (Definition 5 of the appendix). Zhang [25] shows that m-separation and d-separation are in fact equivalent (see Lemma 1 of the appendix), such that given V Z ⊂ V and {V i , V j } ⊂ V , the following holds:EQUATIONwhere(• |= m M G V • | •) denotes m-separation relative to the graph M G V .Just like with DAGs, MAGs that imply the same set of conditional independencies define an equivalence class. Usually, the common structure of these graphs is represented by partial ancestral graphs (PAGs, Definition 7 of the appendix). We useP M G V to denote the PAG relative to M G V .Problem definition. In this work, our goal is to provide theoretical guarantees for the identifiability of the Markov equivalence class of the marginal graph M G V and its direct causal effects with the score, where variables V i are defined according to Equation (4).Without further assumptions on the data-generating process, we can identify the graph M G V only up to its partial ancestral graph, as discussed in the next section.Closely related works. Causal discovery with latent variables have been first studied in the context of constraint-based approaches with the FCI algorithm [11] , which shows the identifiability of the equivalence class of a marginalized graph via conditional independence testing. The RCD and CAM-UV [12, 13] approaches instead demonstrate the inferrability of directed causal edges via regression and residuals independence testing. Both methods rely on strong assumptions on the causal mechanisms: their theoretical guarantees apply to models where the effects are generated by a linear (RCD) or nonlinear (CAM-UV) additive contribution of each cause. Our work demonstrates that using the score function for causal discovery unifies and generalizes these results, presenting an alternative to conditional independence testing for constraint-based methods, and being agnostic about the class of causal mechanisms of the observed variables, under the weaker requirement of additivity of the noise terms.3 Theory for a score-based test of separationIn this section, we show that for V ⊆ X generated according to Equation (4) the Hessian matrix of log p(V ) identifies the equivalence class of the marginal MAG M G V . It has already been proven that cross-partial derivatives of the log-likelihood are informative about a set of conditional independence relationships between random variables: Spantini et al. [21] (Lemma 4.1) shows that, givenEQUATIONEquation ( 3) resulting from faithfulness and the directed global Markov property immediately implies that this expression can be used as a test of conditional independence to identify the Markov equivalence class of the graph M G V , as commonly done in constraint-based causal discovery (for reference, see e.g. Section 3 in Glymour et al. [4] ). This result generalizes Lemma 1 of Montagna et al. [16] , where it is used to define constraints to infer edges in the causal structure without latent variables.Proposition 1 (Adapted1 from [21] ). Let V be a set of random variables with strictly positive density generated according to model (4) .For each set V Z ⊆ V of nodes in M G V such that {V i , V j } ⊆ V Z , the following holds for each supported value v Z : ∂ 2 ∂V i ∂V j log p(v Z ) = 0 ⇐⇒ V i |= m M G V V j |V Z \ {V i , V j }.The result of Proposition 1 presents an alternative to conditional independence testing in constraintbased approaches to causal discovery, showing that the equivalence class of the graph M G V can be identified using the cross partial derivatives of the log-likelihood as a test of conditional independence between variables, much in the spirit of the Fast Causal Inference algorithm [11] . Identifying the Markov equivalence class is the most we can hope to achieve without further hypotheses. As we will see in the next section, the score function can also help leverage additional restrictive assumptions on the causal mechanisms of Equation ( 4) to identify direct causal effects. 
theory assumptions and proofs: 
 Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: All our theoretical results make explicit the assumptions for which they are valid. Plus, we in section Appendix B we provide the proofs of our theoretical results.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 crossreferenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to 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. 
a theory of identifiability from the score: 
 In this section, we show that, under additional assumptions on the data-generating process, we can identify the direct causal relations that are not influenced by unobserved variables, as well as the presence of unobserved active paths (Definition 5) between nodes in the marginalized graph M G V .As a preliminary step before diving into causal discovery with latent variables, we show how the properties of the score function identify edges in directed acyclic graphs, that is in the absence of latent variables (when U = ∅ and G = M G V ). The goal of the next section is two-sided: first, it introduces the fundamental ideas connecting the score function to causal discovery that also apply to hidden variable models, second, it extends the existing theory of causal discovery with score matching to additive noise models with both linear and nonlinear mechanisms.In this section, we summarise and extend the theoretical findings presented in Montagna et al. [17] ,where the authors show how to derive constraints on the score function that identify the causal order of the DAG G where all the variables in the set X are observed. Define the structural relations of (1) as:EQUATIONwith three times continuously differentiable mechanisms h i , noise terms centered at zero, and strictly positive density p X . Given the Markov factorization of Equation (2), the components of the score function ∇ log p(x) are:EQUATIONwhere CH G i denotes the set of children of node X i . We observe that if a node X s is a sink, i.e. a node satisfying CH G s = ∅, then the summation over the children vanishes, implying that:EQUATIONThe key point is that the score component of a sink node is a function of its structural equation noise term, such that one could learn a consistent estimator of ∂ Xs log p X from a set of observations of the noise term N s . Given that, in general, one has access to X samples rather than observations of the noise random variables, authors in Montagna et al. [17] show that N s of a sink node can be consistently estimated from i.i.d. realizations of X. For each node X 1 , . . . , X k , we define the quantity:EQUATIONwhere X \Xi are the random variables in the setX \ {X i }. E[X i |X \Xi ]is the optimal least squares predictor of X i from all the remaining nodes in the graph, and R i is the regression residual. For a sink node X s , the residual satisfies:EQUATIONwhich can be seen by rewritingE[X s |X \Xs ] = h s (X PA G s ) + E[N s |X DE G s , X ND G s ] = h s (X PA G s ) + E[N s ],where X DE G s and X ND G s denotes the descendants and non-descendants of X s , respectively. Equations ( 9) and (11) together imply that the score ∂ Ns log p(N s ) is a function of R s , such that it is possible to find a consistent approximator of the score of a sink from observations of R s .Proposition 2 (Generalization of Lemma 1 in Montagna et al. [17] ). Let X be a set of random variables, generated by a restricted additive noise model (Definition 9) with structural equations (7) , and let X j ∈ X. Consider r j in the support of R j . Then:EQUATIONOur result generalizes Lemma 1 in Montagna et al. [17] , as they assume X generated by an identifiable additive noise model with nonlinear mechanisms. Instead, we remove the nonlinearity assumption and make the weaker hypothesis of a restricted additive noise model, which is provably identifiable [9] , in the formal sense defined in the appendix (Definition 8). This result doesn't come as a surprise, given the previous findings of Ghoshal and Honorio [14] showing that the score infers linear non-Gaussian additive noise models: Proposition 2 provides a unifying and general theory for the identifiability of models with potentially mixed linear and nonlinear mechanisms.Based on these insights, Montagna et al. [17] propose the NoGAM algorithm to exploit the condition in (12) for identifying the causal order of the graph: being E [∂ Xi log p(X) | R i ] the optimal least squares estimator of the score of node X i from R i , a sink node is characterized as theargmin i E [E [∂ Xi log p(X) | R i ] -∂ Xi log p(X)] 2, where in practice the residuals R i , the score components and the least squares estimators are replaced by their empirical counterparts. After a sink node is identified, it is removed from the graph and assigned a position in the order, and the procedure is iteratively repeated up to the source nodes. Being the score estimated by score matching techniques [23] , we usually make reference to score matching-based causal discovery.In the next section, we show how we can generalize these results to identify direct causal effects between a pair of variables in the marginal MAG M G V when U ̸ = ∅We now introduce the last of our main theoretical results, that is: given a pair of nodes V i , V j that are adjacent in the graph M G V with U ̸ = ∅, we can use the score function to identify the presence of a direct causal effect between V i and V j , or that of an active path that is influenced by unobserved variables. Given that the causal model of Equation ( 4) ensures identifiability only up to the equivalence class, we need additional restrictive assumptions. In particular, we enforce an additive noise model with respect to both the observed and unobserved noise variables. This corresponds to an additive noise model on the observed variables with the noise terms recentered by the latent causal effects.Assumption 1 (SCM assumptions). The set of structural equations of the observable variables specified in (4) is now defined as:EQUATIONassuming the mechanisms f i to be of class C 3 (R|V PA G i |), and mutually independent noise terms with strictly positive density function. The N i 's are assumed to be non-Gaussian when f i is linear in some of its arguments.Crucially, our hypothesis is weaker than those required by two state-of-the-art approaches, CAM-UV [13] and RCD [12] : CAM-UV assumes a Causal Additive Model (CAM) with structural equations with nonlinear mechanisms in the formV i := k∈PA G i f ik (V k ) + U i k g ik (U i k ) + N i ,requires an additive noise model with linear effects of both the latent and observed causes. Thus, our model encompasses and extends the nonlinear and linear settings of CAM-UV and RCD, such that the theory developed in the remainder of the section is valid for a broader class of causal models.Our first step is rewriting the structural relations in (13) as:EQUATIONwhich provides an additive noise model in the form of (7) . Next, we define the following regression residuals for any nodeEQUATIONwhere V Z\{k} denotes the set of random variablesV Z \ {V k }.Given these definitions, we are ready to show how directed edges, and the presence of unobserved variables can be identified from the score of linear and nonlinear additive noise models.Consider V i , V j adjacent nodes in the PAG P M G V : we want to investigate when a direct causal effect V i ∈ V PA G j can be identified from the score. We make the following observations: forEQUATIONwhere we useV PA G j |= G d U j to write E[ Ñj |V Z\{j} ] = E[ Ñj ]. Moreover, we note that V j is a sink node relative to M G V Z , the marginalization of G onto V Z .In analogy to the case without latent variables, we can show that ∂ Vj log p(V Z ) is a function of Ñj , the error term in the additive noise model of Equation (14) , such that the score of V j can be consistently predicted from observations of the residual R j (V Z ).Proposition 3. Let X be generated by a restricted additive noise model with structural equations (7) ,and causal graph G. Consider V i , V j adjacent in M G V , marginalization of G. Further, assume that the score component ∂ Vj log p(V Z ) is not constant for uncountable values of V Z . (i) Let V Z = V PA G j ∪ {V i , V j }, and r j ∈ R in the support of R j (V Z ). Then: V PA G j |= d G U j ∧ V i ∈ V PA G j ⇐⇒ E[∂ Vj log p(V Z ) -E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 = 0. (ii) Let V Z ⊆ V , such that {V i , V j } ⊆ V Z . Then: V PA G j ̸ |= d G U j ∨ V i ̸ ∈ V PA G j ⇐⇒ E[∂ Vj log p(V Z ) -E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 ̸ = 0.Intuitively, the proposition has two essential implications. Part (i) provides the condition for the identifiability of the potential direct causal effect between a pair V i , V j , that is, when the association between V j and its observed parents is not influenced by active paths that involve latent variables.This condition is necessary: given an active path such thatV PA G j ̸ |= d G U j, the score could not identify a direct causal effect V i → V j , which is the content of the second part of the proposition.We have established theoretical guarantees of identifiability for linear and nonlinear additive noise models, even in the presence of hidden variables: we find that the score function is a means for the identifiability of all direct parental relations that are not influenced by unobserved variables; all the remaining arrowheads of the edges in the graph M G V are identified no better than in the equivalence class. Based on these insights, we propose AdaScore, a score matching-based algorithm for the inference of Markov equivalence classes, direct causal effects, and the presence of latent variables.Building on our theory, we propose AdaScore, a generalization of NoGAM to linear and nonlinear additive noise models with latent variables. The main strength of our approach is its adaptivity with respect to structural assumptions: based on the user's belief about the plausibility of several modeling assumptions on the data, AdaScore can output an equivalence class (using the condition of Proposition 1 instead of conditional independence testing in an FCI-like algorithm), a directed acyclic graph (as in NoGAM), or a mixed graph, accounting for the presence of unobserved variables.We now describe the version of our algorithm whose output is a mixed graph, where we rely on score matching estimation of the score and its Jacobian (Appendix C.2). At an intuitive level, we find unoriented edges using Proposition 1, i.e. checking for dependencies in the form of non-zero entries in the Jacobian of the score via hypothesis testing on the mean, and find the edges' directions via the condition of Proposition 3, i.e. by estimating residuals of each node X i and checking whether they can correctly predict the i-th score entry (the vanishing mean squared errors are verified by hypothesis test of zero mean). It would be tempting to simply find the skeleton (i.e. the graphical representation of the constraints of an equivalence class) first via the well-known adjacency search of the FCI algorithm and then iterate through all neighborhoods of all nodes to orient edges using Proposition 3. This would be prohibitively expensive, as finding the skeleton is well-known to have super-exponential computational complexity [11] . Instead, we propose an alternative solution: exploiting the fact that some nodes may not be influenced by latent variables, we first use Proposition 2 to find sink nodes that are not affected by latents (using hypothesis testing to find vanishing mean squared error in the score predictions from the residuals), in the spirit of the NoGAM algorithm. If there is such a sink, we search all its adjacent nodes via Proposition 1 (plus an optional pruning step for better accuracy, Appendix C.2), and orient the inferred edges towards the sink. Else, if no sink can be found, we pick a node in the graph and find its neighbors by Proposition 1, orienting its edges using the condition in Proposition 3 (score estimation by residuals under latent effects). This way, we get an algorithm that is polynomial in the best case (Appendix C.3). Details on AdaScore are provided in Appendix C, while a pseudo-code summary is provided in the Algorithm 1 box.Algorithm 1 Simplified pseudo-code of AdaScore while nodes remain do if Proposition 3 finds a sink with all parents observed then add edges from adjacent nodes to sink else pick some remaining node V i ∈ V prune neighbourhood of V i using Proposition 1 orient edges adjacent to V i using Proposition 3 if V i has outgoing directed edge to some V j ∈ V then continue with V j else remove V i form remaining nodes prune remaining bidirected edges using Proposition 1Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a 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 as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We use the causally2 Python library [26] to generate synthetic data with known ground truths, created as Erdös-Rényi sparse and dense graphs, respectively with probability of edge between pair of nodes equals 0.3 and 0.5. We sample the data according to linear and nonlinear mechanisms with additive noise, where the nonlinear functions are parametrized by a neural network with random weights, a common approach in the literature [18, 26, 27, 28, 29] . Noise terms are sampled from a uniform distribution in the [-2, 2] range. Hidden causal effects are obtained by randomly picking two nodes and dropping the corresponding column from the data matrix. See Appendix D.1 for further details on the data generation. As metric, we consider the structural Hamming distance (SHD) [30, 31] , a simple count of the number of incorrect edges, where missing and wrongly directed edges count as one error. We fix the level of the hypothesis tests of AdaScore to 0.05, which is a common choice in the absence of prior knowledge. We compare AdaScore to NoGAM, CAM-UV, RCD, and DirectLiNGAM, whose assumptions are detailed in Table 1 . In the main manuscript, we comment on the results on datasets of 1000 observations from dense graphs, with and without latent variables. Additional experiments including those on sparse networks are presented in Appendix E.Our synthetic data are standardized by their empirical variance to remove shortcuts in the data [18, 32] .Discussion. Our experimental results on models without latent variables of Figure 1a show that when causal relations are linear, AdaScore can recover the causal graph with accuracy that is comparable with all the other benchmarks, with the exception of DirectLiNGAM. On nonlinear data AdaScore presents better performance than CAM-UV, RCD, and DirectLiNGAM while being comparable to NoGAM in accuracy. This is in line with our expectations: in the absence of finite sample errors and in the fully observable setting, NoGAM and AdaScore are indeed the same algorithms.When inferring under latent causal effects, Figure 1b , our method performs comparably to CAM-UV and RCD on graphs up to seven nodes while slightly degrading on nine nodes. Additionally, AdaScore outperforms NoGAM in this setting, as we would expect according to our theory. Overall, we observe that our method is robust to a variety of structural assumptions, with accuracy that is often comparable and sometimes better than competitors (as in nonlinear observable settings). We remark that although AdaScore does not clearly outperform the other baselines, its broad theoretical guarantees of identifiability are not matched by any available method in the literature; this makes it an appealing option for inference in realistic scenarios that are hard to investigate with synthetic data, where the structural assumptions of the causal model underlying the observations are unknown. Figure 1 : Empirical results on dense graphs with different numbers of nodes, on fully observable (no hidden variables) and latent variable models. We report the SHD accuracy (the lower, the better). We note that DirectLiNGAM is surprisingly robust to different structural assumptions, and AdaScore is generally comparable or better (as in nonlinear observable data) than the other benchmarks.Table 1 : Experiments causal discovery algorithms. The content of the cells denotes whether the method supports (✓) or not (✗) the condition specified in the corresponding row.Linear additive noise model✗ ✓ ✗ ✓ ✓Nonlinear additive noise model✗ ✗ ✓ ✗ ✓ Nonlinear CAM ✓ ✗ ✓ ✗ ✓ Latent variables effects ✓ ✓ ✗ ✗ ✓ Output Mixed Mixed DAG DAG Mixed 
conclusion: 
 The existing literature on causal discovery shows a connection between score matching and structure learning in the context of nonlinear ANMs: in this paper, (i) we formalize and extend these results to linear SCMs, and (ii) we show that the score retains information on the causal structure even in the presence of unobserved variables. Additionally, while previous works posit the accent on finding the causal order through the score, we study its potential to identify the Markov equivalence class with a constraint-based strategy that does not explicitly require tests of conditional independence, as well as to identify direct causal effects. Our theoretical insights result in AdaScore: unlike existing approaches for the estimation of causal directions, our algorithm provides theoretical guarantees for a broad class of identifiable models, namely linear and nonlinear, with additive noise, in the presence of latent variables. Even though AdaScore does not clearly outperform the existing baselines on our synthetic benchmark, its adaptivity to different structural hypotheses is a step towards causal discovery that is less reliant on prior assumptions, which are often untestable and thus hindering reliable inference in real-world problems. While we do not touch on the task of causal representation learning [33] , where causal variables are learned from data, we believe this is a promising research direction in relation to our work due to the specific interplay between score-matching estimation and generative models.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: In our experiments section 5, we present all the necessary details on the data generation procedure and a description of the empirical results that are necessary for understanding our findings. Additionally, a comprehensive overview of our experimental design is presented in Appendix D.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 a 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, or as supplemental material. 
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 report all our experimental results in the form of boxplots.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, confidence intervals, or statistical significance tests, at least for the experiments that support the 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, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• 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 a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error 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. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?• Researchers should communicate the details of the dataset/code/model 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. 
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 as details about compensation (if any)?Answer: [NA] Justification: We do not work with human subjects or crowdsourcing.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.• According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor 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 NeurIPS 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. 
 
 Appendix: In this section, we provide a collection of results and definitions relevant to the theory of this paper.Let X = X 1 , . . . , X d a set of random variables. A graph G = (X, E) consists of finitely many nodes or vertices X and edges E. We now provide additional definitions, separately for directed acyclic and mixed graphs.Directed acyclic graph. In a directed graph, nodes can be connected by a directed edge (→), and between each pair of nodes there is at most one directed edge. We say that X 1 is a parent of X j ifX i → X j ∈ E,in which case we also say that X j is a child of X i . Two nodes are adjacent if they are connected by an edge. Three nodes are called a v-structure if one node is a child of the other two, e.g. asX i → X k ← X j is a collider. A path in G is a sequence of at least two distinct vertices X i1 , .. . , X im such that there is an edge between X i k and X i k+1 . If X i k → X i k+1 for every node in the path, we speak of a directed path, and call X i k an ancestor of X i k+1 , X i k+1 a descendant ofX i k .Given the set DE G i of descendants of a node X i , we define the set of non-descendants of X i asND G i = X \ (DE G i ∪{X i }).A node without parents is called a source node. A node without children is called a sink node. A directed acyclic graph is a directed graph with no cycles.Mixed graph. In a mixed graph nodes can be connected by a directed edge (→) or a bidirected edge (↔), and between each pair of nodes there is at most one directed edge. Two vertices are said to be adjacent in a graph if there is an edge (of any kind) between them. The definitions of parent, child, ancestor, descendant, path provided for directed acyclic graph also apply in the case of mixed graphs. Additionally, X i is a spouse of X j (and vice-versa) ifX i ↔ X j ∈ E. An almost directed cycle occurs when X i ↔ X j ∈ E and X i is an ancestor of X j in G.For ease of reference from the main text, we separately provide the definition of inducing paths and ancestors in directed acyclic graphs.Definition 2 (Ancestor). Consider a DAG G with set of nodes X, and X i , X j elements of X. We say that X i is an ancestor of X j if there is a directed path from X i to X j in the graph, as inX i → . . . → X j .Definition 3 (Inducing path). Consider a DAG G with set of nodes X, and Y, Z disjoint subsets such that X = Y ∪Z. We say that there is an inducing path relative to Z between the nodes Y i , Y j if every node on the path that is not in Z ∪ {Y i , Y j } is a collider on the path (i.e. for each Y k ∈ Y on the path the sequence Y i . . . → Y k ← . . . Y j appears) and every collider on the path is an ancestor of Y i or Y j .One natural way to encode inducing paths and ancestral relationships between variables is represented by maximal ancestral graphs.. A maximal ancestral graph (MAG) is a mixed graph such that:1. there are no directed cycles and no almost directed cycles;2. there are no inducing paths between two non-adjacent nodes.Next, we define conditional independence in the context of graphs.Definition 5 (m-separation). Let M be a mixed graph with nodesX. A path π in M between X i , X j elements of X is active w.r.t. Z ⊆ X \ {X i , X j } if:1. every non-collider on π is not in Z 2. every collider on π is an ancestors of a node in Z.X i and X j are said to be m-separated by Z if there is no active path between X i and X j relative to Z.Two disjoint sets of variables W and Y are m-separated by Z if every variable in W is m-separated from every variable in Y by Z.If m-separation is applied to DAGs, it is called d-separation.The set of directed acyclic graphs that satisfy the same set of conditional independencies form an equivalence class, known as the Markov equivalence class.Definition 6 (Markov equivalence class of a DAG). Let G be a DAG with nodes X. We denote with[G] the Markov equivalence class of G. A DAG G with nodes X is in [G]if the following conditions are satisfied for each pair X i , X j of distinct nodes in X:• there is an edge between X i , X j in G if and only if there is an edge between X i , X j in G;• let Z ⊆ X \ {X i , X j }. Then X i |= d G X j |Z ⇐⇒ X i |= d G X j |Z;• let π be a path between X i and X j . X k is a collider in the path π in G if and only if it is a collider in the path π in G.In summary, graphs in the same equivalence class share the edges up to direction, the set of dseparations, and the set of colliders.Just as for DAGs, there may be several MAGs that imply the same conditional independence Intuitively, a PAG represents an equivalence class of MAGs by displaying all common edge marks shared by all members of the class and displaying circles for those marks that are not in common.In this section, we provide a proof for equation (5) , stating the equivalence between m-separation and d-separation in a formal sense.Lemma 1 (Adapted from Zhang [25] ). Let G be a DAG with nodes X = V ∪ U , with V and U disjoint sets, and M G V the marginalization of G onto V . For any {V i , V j } ∈ V and V Z ⊆ V \{V i , V j }, the following equivalence holds:V i |= d G V j |V Z ⇐⇒ V i |= m M G V V j |V Z . Proof. The implication V i |= d G V j |V Z =⇒ V i |= m M G V V j |V Z is a direct consequence of Lemma 18from Spirtes and Richardson [34] , where we set S = ∅, since we do not consider selection bias. The implication V i |= d G V j |V Z ⇐= V i |= m M G V V j |V Z followsWe study the identifiability of the additive noise model, reporting results from Peters et al. [9] . We start with a formal definition of identifiability in the context of causal discovery.Definition 8 (Identifiable causal model). Let (X, N, F, p N ) be an SCM with underlying graph G and p X joint density function of the variables of X. We say that the model is identifiable from observational data if the distribution p X can not be generated by a structural causal model with graph G ̸ = G.First, we consider the case of models of two random variablesEQUATION)Condition 1 (Condition 19 of Peters et al. [9] ). Consider an additive noise model with structural equations (17) . The triple (f, p X1 , p N ) does not solve the following differential equation for all pairsEQUATIONHere, ξ := log p X1 , ν := log p N , the logarithms of the strictly positive densities. The argumentsx 2 -f (x 1 ), x 1 , and x 1 of ν, ξ and f respectively, have been removed to improve readability.Next, we show that a structural causal model satisfying Condition 1 is identifiable, as in Definition 8Theorem 1 (Theorem 20 of Peters et al. [9] ). Let p X1,X2 the joint distribution of a pair of random variables generated according to the model of equation ( 17) that satisfies Condition 1, with graph G.Then, G is identifiable from the joint distribution.Finally, we show an important fact, holding for identifiable bivariate models, which is that the score∂ ∂X1 log p(x 1 , x 2 ) is nonlinear in x 1 .Lemma 2 (Sufficient variability of the score). Let p X1,X2 the joint distribution of a pair of random variables generated according to a structural causal model that satisfies Condition 1, with graph G.Then:∂ ∂X 1 (ξ ′ (x 1 ) -f ′ (x 1 )ν ′ (x 2 -f (x 1 ))) ̸ = 0,for all pairs (x 1 , x 2 ).Proof. By contradiction, assume that there exists(x 1 , x 2 ) such that ∂ ∂X1 (ξ ′ (x 1 ) -f ′ (x 1 )ν ′ (x 2 - f (x 1 ))) = 0. Then: ∂ ∂X 1   ∂ 2 ∂X 2 1 π(x 1 , x 2 ) ∂ 2 ∂X1∂X2 π(x 1 , x 2 )   = 0,where π(x 1 , x 2 ) = log p(x 1 , x 2 ). By explicitly computing all the partial derivatives of the above equation, we obtain that equation 18 is satisfied, which violates Condition 1.These results guaranteeing the identifiability of the bivariate additive noise model can be generalized to the multivariable case, with a set of random variables X = {X 1 , . . . , X k } that satisfy:EQUATIONwhere G is the resulting causal graph directed and acyclic. The intuition is that, rather than studying the multivariate model as a whole, we need to ensure that Condition 1 is satisfied for each pair of nodes, adding restrictions on their marginal conditional distribution.Definition 9 (Definition 27 of Peters et al. [9] ). Consider an additive noise model with structural equations (19) . We call this SCM a restricted additive noise model if for all X j ∈ X, X i ∈ X PA G j , and all sets X S ⊆ X, S ⊂ N, withX PA G j \ {X i } ⊆ X S ⊆ X G NDj \ {X i , X j }, there is a value x S with p(x S ) > 0, such that the triplet(f j (x PA G j \{i} , •), p Xi|X S =x S , p Nj ) satisfies Condition 1. Here, f j (x PA G j \{i} , •) denotes the mechanism function x i → f j (x PA G j ).Additionally, we require the noise variables to have positive densities and the functions f j to be continuous and three times continuously differentiable.Then, for a restricted additive noise model, we can identify the graph from the distribution.Theorem 2 (Theorem 28 of Peters et al. [9] ). Let X be generated by a restricted additive noise model with graph G, and assume that the causal mechanisms f j are not constant in any of the input arguments, i.e. for X i ∈ X PA G j , there existx i ̸ = x ′ i such that f j (x PA G j \{i} , x i ) ̸ = f j (x PA G j \{i} , x ′ i ).Then, G is identifiable.We state several results that hold for a pair of random variables that are not connected by an active path that includes unobserved variables (active paths are introduced in Definition 5). For the remainder of the section, let V, U be a pair of disjoint sets of random variables, X = V ∪ U generated according to the structural causal model defined by the set of equations ( 1), G the associated causal graph, andM G V the marginalization onto V .The first statement provides under which condition the unobserved parents of two variables in the marginal MAG are mutually independent random vectors.Lemma 3. Let V j ∈ V , and Z ⊂ N such that V Z = V PA G j ∪ {V j }. Assume V PA G j |= d G U j . Then U j |= d G U Z k for each index Z k ̸ = j.Proof.The assumptionV PA G j |= d G U j implies that there is no active path in G between nodes in V PA G jand nodes in U j . Given that for eachZ k ∈ Z, Z k ̸ = Z, nodes in U Z k are direct causes of at least one node in V PA G j, any active path between nodes in U Z k and nodes in U j would also be an active path between V PA G j and U j , which is a contradiction. HenceU j |= d G U Z k .The previous lemmas allow proving the following result, which will be fundamental to demonstrate the theory of Proposition 3.Lemma 4. Let V j ∈ V , and Z ⊂ N such that V Z = V PA G j ∪ {V j }. Assume V PA G j |= d G U j . W.l.o.g., let the j-th element of V Z be V Zj = V j . Denote as U Z the set of unobserved parents of nodes in V Z , and U Z\{j} the unobserved parents of nodes in V Z\{j} := V Z \ V j . Then, the following holds for each v Z , u Z values:log p(v Z ) = log p(v j |v PA G j ) + log Q(v Z ),whereQ(v Z ) = u Z\{j} p(u Z\{j} ) |Z| k̸ =j p(v Z k |v Z1 , . . . , v Z k-1 , u Z k ).Proof. By the law of total probability and the chain rule, we can write p(v Z ) as:EQUATIONBy Lemma 3, U Zj |= U Z k , k ̸ = j, where U Z k denotes unobserved parents of the node V Z k . Then, we can factorize p(u) = p u Zj p u Z\{j} . Plugging the factorization in equation ( 20) we findp(v Z ) = u p u Zj p u Z\{j} p(v Zj |u, v Z\{j} )p(v Z\{j} |u) = u p u Zj p u Z\{j} p(v Zj |u Zj , v PA G Z j )p(v Z\{j} |u),where the latter equation comes from the global Markov property on the graph G. Further, by assump-tion of V PA G j |= d G U j , we know that U Zj |= V Z k , k ̸ = j, such that p(v Z\{j} |u) = p(v Z\{j} |u Z\{j} ).Then:p(v Z ) = u p u Zj p u Z\{j} p(v Zj |u Zj , v PA G Z j )p(v Z\{j} |u Z\{j} ) = u Z j p u Zj p(v Zj |u Zj , v PA G Z j ) u Z\{j} p u Z\{j} p(v Z\{j} |u Z\{j} ) = p(v Zj |v PA G Z j ) u Z\{j} p u Z\{j} p(v Z\{j} |u Z\{j} ),which proves the claim.Intuitively, Lemma 4 shows that given a node V j without children and bidirected edges in a marginalized graph M G V Z , the kernel of node V j in the Markov factorization of p(v Z ) is equal to the kernel of the same node in the Markov factorization of p(x) of equation ( 2), relative to the graph without latent confounders G.Proof of Proposition 1. Observe that∂ 2 ∂V i ∂V j log p(v Z ) = 0 ⇐⇒ V i |= d G V j |V Z \ {V i , V j } ⇐⇒ V i |= m M G V V j |V Z \ {V i , V j },where the first equivalence holds by a combination of the faithfulness assumption with the global Markov property, as explicit in equation ( 3), and the second due to Lemma 1. Then, the claim is proven.Proof. The forward direction is immediate from equation ( 9) and R j = N j , when X j is a sink (equation ( 11)). Thus, we focus on the backward direction. GivenE E ∂ Xj log p(X) | R j = r j -∂ Xj log p(X) 2 = 0,we want to show that X j has no children, which we prove by contradiction.Let us introduce a function q : R → R such that:E ∂ Xj log p(X) | R j = r j = q(r j ),and s j : R |X| → R, s j (x) = ∂ Xj log p(x).The mean squared error equal to zero implies that s j (X) is a constant, once R j is observed. Formally, under the assumption of p(x) > 0 for each x ∈ R k , this implies thatp(s j (x) ̸ = q(R j )|R j = r j ) = 0, ∀x ∈ R k .By contradiction, we assume that X j is not a leaf, and want to show that s j (X) is not constant in X, given R j fixed. Let X i such that X j ∈ X PA G i . Being the structural causal model identifiable, there is no model with distribution p X whose graph has a backward edge X i → X j : thus, the Markov factorization of equation ( 2) is unique and implies:∂ Xj log p(X) = ∂ Nj log p(N j ) - k∈CH G j ∂ Xj h k (X PA k )∂N k log p(N k ).We note that, by definition of residual in equation ( 10), R j = r j fixes the following distance:R j = N j -E[N j |X \Xj ].Hence, conditioning on R j doesn't restrict the support of X: given R j = r j , for any x \Xj (value of the vector of elements in X \ {X j }), ∃n j with p(n j > 0) (by the hypothesis of strictly positive densities of the noise terms) that satisfiesr j = n j -E[N j |x \Xj ].Next, we condition on all the parents of X i , except for X j , to reduce our problem to the simpler bivariate case. Let S ⊂ N and X S ⊆ X such that X PA G i \ {X j } ⊆ X S ⊆ X ND G i \ {X i , X j }, and consider x S such that p(x S > 0). Let X PA G i = x PA G i hold under X S = x S . We define X j |x s := X j |(X S = x S ),g i (x j ) = h i (x PA G i \{Xj } , x j ). Consider X i = x i , and the pair of values (x j , x * j ) such that x j ̸ = x * j andν ′′ Ni (x i -g i (x j ))g ′ i (x j ) ̸ = 0, ν ′′ Ni (x i -g i (x * j ))g ′ i (x * j ) ̸ = 0, where we resort to the usual notation ν Ni := log p Ni . By Lemma 2, (x i , x j ) and (x i , x * j ) satisfy:∂ Xj (ξ ′ (x j ) -ν ′ Ni (x i -g i (x j ))g ′ i (x j )) ̸ = 0, ∂ Xj (ξ ′ (x * j ) -ν ′ Ni (x i -g i (x * j ))g ′ i (x * j )) ̸ = 0, where ξ := log p Xj |x s. Thus, we can fix x j and x * j (which are arbitrarily chosen) such thatEQUATIONFixing X |x S ,xj = x and X | x S ,x * j = x * , where the two values differ only in their j-th component, we find the following difference:s j (x) -s j (x * ) = ∂ Xj (ξ ′ (x j ) -ν ′ Ni (x i -g i (x j ))g ′ i (x j )) -∂ Xj (ξ ′ (x * j ) -ν ′ Ni (x i -g i (x * j ))g ′ i (x * j )), which is different from 0 by equation ( 21). This contradicts the fact that the score s j is constant once R j is fixed, which proves our claim.In this proof, we use several ideas from the demonstration of Proposition 2. We demonstrate the forward and the backward parts of the two statements separately.Proof of part (i), forward direction. Given V Z = V PA G j ∪ {V i , V j } and r j ∈ R in the image of R j ,we want to show:V PA G j |= d G U j ∧ V i ∈ V PA G j =⇒ E[∂ Vj log p(V Z ) -E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 = 0.By Lemma 4, the score of V j is∂ Vj log p(V Z ) = ∂ Vj log p(V j |V PA G j ) + ∂ Vj log Q(V Z ) = log p( Ñj ),for some Q map acting on V Z . The latter equality holds because all variables in V Z are nondescendants of V j , such that ∂ Vj Q(V Z ) = 0. Further, by equation ( 16) we know thatR j (V Z ) = Ñj + c,where c = -E[ Ñj ] is a constant. It follows that the least square estimator of the score of V j from R j (V Z ) satisfies the following equation:E[∂ Vj log p(V Z )|R j (V Z )] = E[∂ Vj log p( Ñj )| Ñj ] = ∂ Vj log p( Ñj ),where the first equality holds becauseEQUATIONwe want to show thatV PA G j |= d G U j ∧ V i ∈ V PA G j, meaning that there is a direct causal effect that is not biased by unobserved variables. We provide the proof by contradiction, in analogy to the demonstration of the backward direction of Proposition 2.j : R |V Z | → R, s j (v Z ) = ∂ Vj log p(V Z ).The mean squared error equal to zero implies that s j (V Z ) is constant in V Z , once R j is observed. By contradiction, we assume thatV PA G j ̸ |= d G U j ∨ V i ̸ ∈ V PA G j, and want to show that s j (V Z ) is not constant in V Z , given R j fixed. In this regard, we make the following observation: by definition of residual in equation ( 15), R i (V Z ) = r i fixes the following distance:R j (V Z ) = Ñj -E[ Ñj |V Z\{j} ].Hence, conditioning on R j (V Z ) doesn't restrict the support of V Z : given R j (V Z ) = r j , ∃ñ j with p(ñ j ) > 0 (by assumption of strictly positive densities p Nj and p X ), that satisfiesr j = ñj -E[ Ñj |v Z\{j} ],for all v Z\{j} . Hence, the random variable V Z |R j (V Z ) = r j has strictly positive density on all points v Z where p V Z (v Z ) > 0. Now, consider v Z and v * Z , taken from the set of uncountable values such that the score s j function is not a constant, meaning thats j (v Z ) ̸ = s j (v * Z ), where V Z is sampled givenR j (V Z ) = r j .Given that different v Z and v * Z are selected from an uncountable subset of the support, we conclude that the scores j |(R j (V Z ) = r j ) = ∂ Vj log p(V Z |R j (V Z ) = r j )is not a constant for at least an uncountable set of points, which contradicts equation (22) .Proof of part (ii), forward direction. Given that V i is connected to V j in the marginal MAG and thatV PA G j ̸ |= d G U j ∨ V i ̸ ∈ V PA G j, we want to show that for each V Z ⊆ V with {V i , V j } ⊆ V Z , the following holds:EQUATIONLet us introduce h : R → R such that:E[∂ Vj log p(V Z )|R j (V Z ) = r j ] = h(r j ),and further define:s j (V Z ) = ∂ Vj log p(V Z ).Having the mean squared error in equation ( 23) equals zero implies that s j (V Z ) is a constant, once R j (V Z ) is observed. Thus, the goal of the proof is to show that there are values of V Z such that the score is not a constant once R j is fixed. By definition of residual in equation ( 15), R j (V Z ) = r j fixes the following distance:R j (V Z ) = Ñj -E[ Ñj |V Z\{j} ].Hence, conditioning on R j (V Z ) doesn't restrict the support of V Z : given R j (V Z ) = r j , ∃ñ j with p(ñ j ) > 0 (by assumption of positive density of the noise N j on the support R), that satisfiesr j = ñj -E[ Ñj |v Z\{j} ],for all v Z\{j} . Hence, the random variable V Z |R j (V Z ) = r j has strictly positive density on all points v Z where p V Z (v Z ) > 0. Now, consider v Z and v * Z , taken from the set of uncountable values such that the score s j function is not a constant, meaning thats j (v Z ) ̸ = s j (v * Z ), where V Z is sampled given R j (V Z ) = r j .Given that different v Z and v * Z are selected from an uncountable subset of the support, we conclude that the scores j |(R j (V Z ) = r j ) = ∂ Vj log p(V Z |R j (V Z ) = r j )is not a constant for at least an uncountable set of points, such that the claim follows.Proof of part (ii), backward direction.Given that E[∂ Vj log p(V Z ) -E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 ̸ = 0 for all V Z ⊆ V such that {V i , V j } ∈ V Z ,and given V i and V j adjacent in the marginal MAG, we want to show thatV PA G j ̸ |= d G U j ∨ V i ̸ ∈ V PA G j .The prove comes easily by contradiction: say thatV PA G j |= d G U j ∧ V i ∈ V PA G j. Then, by the forward direction of part (i) of Proposition 3, we know thatV Z = V PA G j ∪ {V j } satisfies E[∂ Vj log p(V Z ) - E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 = 0, leading to a contradiction.In Proposition 1 we have seen that score matching can detect m-separations and therefore the skeleton of the PAG describing the data. If one is willing to make the assumptions required for Proposition 3 it could be desirable to use this to orient edges, since the interpretation of PAG edges might be cumbersome for people not familiar with ancestral models. Therefore, one could simply find the skeleton of the PAG using the fast adjacency search [5] and then orient the edges by applying Proposition 3 on every subset of the neighbourhood of every node. This would yield a very costly algorithm. But if we make the assumptions required to orient edges with Proposition 3 we can do a bit better. In Algorithm 2 we present an algorithm that still has the same worst case runtime but runs polynomially in the best case. The main intuition is that we iteratively remove irrelevant nodes in the spirit of the original SCORE algorithm [15] . To this end, we first check if the is any unconfounded sink if we consider the set of all remaining variables. If there is one, we can orient its parents and ignore it afterwards. If there is no such set, we need to fall back to the procedure proposed above, i.e.we need to check the condition of Proposition 3 on all subsets of the neighbourhood of a node, until we find no node with a direct outgoing edge. In Proposition 4 we show that this way we do not fail orient edge or fail to remove any adjacency. In the following discussion, we will use the notationδ i (X Z ) := E[∂ Vj log p(V Z ) -E[∂ Vj log p(V Z )|R j (V Z ) = r j ]] 2 ,for the second residual from Proposition 3 and alsoδ i,j (X Z ) := ∂ 2 ∂V i ∂V j log p(v Z )for the cross-partial derivative, where X i , X j ∈ V and Z ⊆ V .Proposition 4 (Correctness of algorithm). Let X = V ∪U be generated by the SCM in Equation ( 4)with non-constant scores for uncountably many values. Let G X be the causal DAG of X and G V be the marginal MAG of G X . Then Algorithm 2 outputs a directed edge from X i ∈ V to X j ∈ V iff there is a direct edge in G X between them and no unobserved backdoor path w.r.t. U . Further, the output of Algorithm 2 has the same skeleton as G V .Proof. We proof the statement by induction over the steps of the algorithm. Let S be the set of remaining nodes in an arbitrary step of the algorithm. Our induction hypothesis is that for X i , X j ∈ S and X k ∈ B i we have 1. X i is an unconfounded sink w.r.t. to some set S ′ ⊆ S iff X i is an unconfounded sink w.r.t. someS ′′ ⊆ V 2. if there is no S ′ ⊆ V \ {X i , X j } such that X i |= X j | S ′ then X j ∈ B iClearly, this holds in the initial step as S = V .Suppose we find δ i (X S ) = 0 for X i ∈ S. If X i has at least one adjacent node in M G V , by Proposition 3, we know that X i does not have any children and is also not connected to any other node in S via a hidden mediator or unobserved confounder. This means, all nodes that are not separable from X i must be direct parents of X i , which are by our induction hypothesis 2) the nodes in B i . Since X i does not have children, it also suffices to checkX i |= X j |S \ {X i , X j } for X j ∈ B i(instead of conditioning on all subsets of B i ). So we can already add these direct edges to the output.If, on the other hand, X i has no adjacent nodes in M G V , we haveX i |= X j |S \ {X i , X j } for X j ∈ B i ,so in both cases we add the correct set of parents. Since X i is not an ancestor of any of the nodes in S \ {X i }, X i cannot be a hidden mediator or hidden confounder between nodes in S \ {X i } and conditioning on X i cannot block an open path. Thus, the induction hypothesis still holds in the next step.Suppose now there is no unconfounded sink and we explore X i . By our induction hypothesis 2), B i contains the parents of X i and by Proposition 3 it suffices to only look at subsets of B i to orient direct edges. And also due to the induction hypothesis 2) B i contains all nodes that are not separable from X i . So by adding bidirected edges to all nodes in B i can only add too many edges but not miss some.Algorithm 2 AdaScore Algorithm procedure ADASCORE(p, X 1 , . . . , X d ) S ← {X 1 , . . . , X d } ▷ Remaining nodes E ← { } ▷ Edges for X i ∈ S do B i ← {X 1 , . . . , X d } ▷ Neighbourhoods while S ̸ = ∅ do ▷ While nodes remain if ∃X i ∈ S : δ i (X S ) = 0 then ▷ If there is an unconfounded sink S ← S \ {X i } E ← E ∪ {X j → X i : δ i,j (X S ) ̸ = 0} ▷ Add edges like DAS else for X i ∈ S do for X j ∈ B i do ▷ Prune neighbourhoods if δ i,j (X S ) = 0 then B i ← B i \ {X j } B j ← B j \ {X i } for X j ∈ B i do ▷ Orient edges in B i m i = min S ′ ⊆Bi δ i (X S ′ ∪{Xi} ) m j = min S ′ ⊆Bj δ j (X S ′ ∪{Xj } )) if m i = 0 ∧ m j ̸ = 0 then E ← E ∪ {X j → X i } else if m i ̸ = 0 ∧ m j = 0 then E ← E ∪ {X i → X j } else E ← E ∪ {X i ↔ X j } if ∃X j ∈ B i : (X i → X j ) ∈ E then continue with X j else ▷ X i has no unconfounded outgoing edge S ← S \ {X i } ▷ Remove X i break for X i ↔ X j ∈ E do ▷ Prune bidirected edges if min S ′ ⊆Adj(Xi) δ i,j (X S ′ ∪{Xi} ) = 0 ∨ min S ′ ⊆Adj(Xj ) δ i,j (X S ′ ∪{Xi} ) = 0 then E ← E \ {X i ↔ X j } return ENow it remains to show that the induction hypothesis holds if we set S to S \ {X i }. For 1) we need to show that X i cannot be a hidden mediator or hidden confounder w.r.t. S \ {X i } (since ignoring X i won't change whether there is a direct edge or not). Suppose X i is on a unobserved causal pathX k → • • • → U m → X l with X k , X l ∈ S \ {X i } and U m ∈ X \ (S \ {X i }). This path must have been a unobserved causal path before, unless X i = U m . But then there is a direct edge X i → X l . We would not remove X i from S if this edge was unconfounded, so there must a hidden confounder between X i and X l . But in this case, Proposition 3 wouldn't allow us to direct the edge anyway, sinceV PA l ̸ |= d G U l . Suppose there is confounding path X k ← • • • → U m → X l with X k , X l ∈ S \ {X i } and U m ∈ X \ (S \ {X i }). If X i ̸ = U mthe path was already been a confounding path without X i being unobserved. So again, there must be a confounder between X i and X l , as otherwise we would not remove X i . And analogously to before, we could not have oriented the edge even withX i ∈ S since V PA l ̸ |= d G U l .For 2) we only have to see that we just remove nodes from B i if we found an independence.For |S| < 2, the algorithm enters the final pruning stage. From the discussion above it is clear, that we already have the correct result, up to potentially too many bidirected edges. In the final step we certainly remove all these edges X i ↔ X j , as we check m-separation for all subsets of the neighbourhoods Adj(X i ) and Adj(X j ), which are supersets of the true neighbourhoods.All theoretical results in the paper have assumed that we know the density of our data. Obviously, in practise we have to deal with a finite sample instead. Especially, in Proposition 1 and Proposition 3 we derived criteria that compare random variables with zero. Clearly, this condition is never met in practise. Therefore, we need find ways to reasonably set thresholds for these random quantities.First note, that we use the Stein gradient estimator [35] to estimate the score function. This means especially that for a node V i we get a vectorEQUATIONi.e. an estimate of the score for every one of the m samples. Analogously, we get a m × d × d tensor for the estimates of ∂ 2 ∂Vi∂Vj log p(v).In Proposition 1 we showed that∂ 2 ∂V i ∂V j log p(v Z ) = 0 ⇐⇒ X i |= m M G V V j |V Z \ {V i , V j }.In the finite sample version, we use a one sample t-test on the vector of estimated cross-partial derivatives with the null-hypothesis that the means is zero. Due to the central limit theorem, the sample mean follows approximately a Gaussian distribution, regardless of the true distribution of the observations.For Proposition 3 we need to do some additional steps. Recall, that the relevant quantity in Proposition 3 is the mean squared error of a regression, which is always positive. Therefore, a test for mean zero is highly likely to reject in any case. We decided to employ a two-sample test in a similar (but different) manner as Montagna et al. [17] . As test, we used the Mann-Whitney U-test. Note, that Algorithm 2 employs Proposition 3 in two different ways: first, to decide whether there is an unconfounded sink and second, to orient edges in case there is no unconfounded sink. We pick a different sample as second sample of the Mann-Whitney U-test.Analogously to before, this is a vector with m entries, one for every sample.Note, that in the case where we want to check if there is an unconfounded sink, we do not make any mistake by rejecting too few hypotheses, i.e. if we miss some unconfounded sinks (instead, we only lose efficiency, as we do the costly iteration over all possible sets of parents). Therefore, for this test we chose a a second sample that yields a "conservative" test result.As candidate sink for set S ⊆ V , we pick the node X i = min i mean(δ i (X S )). In fact, we want to know whether the mean of δ i is significantly lower than all other means. But we empirically observed that choosing the concatenated δs of all nodes as second sample makes the test reject with very high probability, which would lead our algorithm to falsely assume the existence of an unconfoudned sink.Instead, we then pick as second "reference node" X j = min j̸ =i mean(δ j (X Z )). We then do the two sample test between δ i (X Z ) and δ j (X Z ). The intuition is that the test will reject the hypothesis of identical means, if X i is an unconfounded sink but X j is not.In the case where we use Proposition 3 to orient edges, we only need to decide whether an not previsouly directed edge X i -X j needs to be oriented one way, the other way, or not at all. Instead, here the issue lies in the fact that we need to iterate over possible sets of parents of the nodes. Let B i be the set of nodes that have not been m-separated from X i by any test so far. We pick the subset Z i = min Z ′ ⊆Bi mean(δ Z ′ i ), i.e. the set with the lowest mean error. We then conduct the test with δ i (X Zi ) and δ j (X Zj ). If there is a directed edge between them, one of the residuals will be significantly lower than the other.Just like Montagna et al. [17] we use a cross-validation scheme to generate the residuals, in order to prevent overfitting. We split the dataset into several equally sized, disjoint subsamples. For every residual we fit the regression on all subsamples that don't contain the respective target.Also, just like in the NoGAM algorithm Montagna et al. [17] we add a pruning step for the directed edges to the end. The idea is to use a feature selection method to remove insignificant edges. Just like In the best case, we always find an unconfounded sink. Then our algorithm reduces to NoGAM.In this section, we present the details of our experiments in terms of synthetic data generation and algorithms hyperparameters.In this work, we rely on synthetic data to benchmark AdaScore's finite samples performance. For each dataset, we first sample the ground truth graph and then generate the observations according to the causal graph.Erdös-Renyi graphs. The ground truth graphs are generated according to the Erdös-Renyi model.It allows specifying the number of nodes and the probability of connecting each pair of nodes). In ER graphs, a pair of nodes has the same probability of being connected.Nonlinear causal mechanisms. Nonlinear causal mechanisms are parametrized by a neural network with random weights. We create a fully connected neural network with one hidden layer with 10 units, Parametric ReLU activation function, followed by one normalizing layer before the final fully connected layer. The weights of the neural network are sampled from a standard Gaussian distribution.This strategy for synthetic data generation is commonly adopted in the literature [26, 18, 28, 29, 27] .Linear causal mechanisms. For the linear mechanisms, we define a simple linear regression model predicting the effects from their causes and noise terms, weighted by randomly sampled coefficients.Coefficients are generated as samples from a Uniform distribution supported in the range [-3, -0.5]∪ [0.5, 3]. We avoid too small coefficients to avoid close to unfaithful datasets Uhler et al. [24] .Noise terms distribution. The noise terms are sampled from a Uniform distribution supported between -2 and 2.Finally, we remark that we standardize the data by their empirical data. This is known to remove shortcuts that allow finding a correct causal order sorting variables by their marginal variance, as in varsortability, described in Reisach et al. [32] , or sorting variables by the magnitude of their score |∂ Xi log p(X)|, a phenomenon known as scoresortability analyzed by Montagna et al. [18] .For AdaScore, we set the α level for the required hypothesis testing at 0.05. For the CAM-pruning step, the level is instead set at 0.001, the default value of the dodidscover Python implementation of the method, and commonly found in all papers using CAM-pruning for edge selection [15, 16, 17, 36] .For the remaining parameters. The regression hyperparameters for the estimation of the residuals are found via cross-validation during inference: tuning is done minimizing the generalization error on the estimated residuals, without using the performance on the causal graph ground truth. Finally, for the score matching estimation, the regularization coefficients are set to 0.001.All experiments have been run on an AWS EC2 instance of type p3.2xlarge. These machines contain Intel Xeon E5-2686-v4 processors with 2.3 GHz and 8 virtual cores as well as 61 GB RAM.All experiments can be run within a day.In this section, we provide additional experimental results. All synthetic data has been generated as described in Appendix D.1.In Figure 1 we have demonstrated the performance of our proposed method on data generated by linear SCMs and non-linear SCMs with additive noise. But Proposition 1 also holds for any faithful distribution generated by an acyclic model. Thus, we employed as mechanism a neural network-based approach similar to the non-linear mechanism described in Appendix D. Instead of adding the noise term, we feed it as additional input into the neural network. Results in this setting are reported in Figure 2 . As neither AdaScore nor any of the baseline algorithms has theoretical guarantees for the orientation of edges in this scenario, we report the F 1 -score (popular in classification problems) w.r.t.to the existence of an edge, regardless of orientation. Our experiments show that AdaScore can, in general, correctly recover the graph's skeleton in all the scenarios, with an F 1 score median between 1 and ∼ 0.75, respectively for small and large numbers of nodes.In this section, we present the experiments on sparse Erdös-Renyi graphs where each pair of nodes is connected by an edge with probability 0.3. The results are illustrated in Figure 3 . For sparse graphs, recovery results are similar to the dense case, with AdaScore generally providing comparable performance to the other methods.In the following series of plots we demonstrate the scaling behaviour of our method w.r.t. to the number of samples. Figure 5 shows results with edge probability 0. Figure 5 : Empirical results on dense graphs with different numbers of samples and seven nodes, on fully observable (no hidden variables) and latent variable models. We report the SHD accuracy (the lower, the better).In this section, we remark the limitations of our empirical study. It is well known that causal discovery lacks meaningful, multivariate benchmark datasets with known ground truth. For this reason, it is common to rely on synthetically generated datasets. We believe that results on synthetic graphs should be taken with care, as there is no strong reason to believe that they should mirror the benchmarked algorithms' behaviors in real-world settings, where often there is no prior knowledge about the structural causal model underlying available observations. method does not provide performance that clearly improves on the existing literature. These limitations of our work are discussed in the discussion of the experiments in Section 5, as well as in the "Limitations" appendix section E.4. Concerning the assumptions required by our method, we thoroughly discuss them in the theoretical sections of the paper, where we define the results that are later used for the definition of the AdaScore. Finally, computational complexity is discussed in Appendix C.3.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a 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 on 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 a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not 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 to address problems of privacy and fairness.• 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 acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.In their Lemma 4.1 Spantini et al.[21] provides the connection between vanishing cross-partial derivatives of the log-likelihood and conditional independence of random variables. Note that this result does not depend on the assumption of a generative model, thus holding beyond the set of structural equations (4). Our result adapts their finding to the case when observations are generated according to a fully observable causal model.https://causally.readthedocs.io/en/latest/ 


Decision: This manuscript was actively discussed at length. First, with the authors during the rebuttal period, and afterwards, again among the reviewers and myself. It is fair to say that all reviewers agree that the manuscript puts forward an interesting idea, and everyone also appreciated that the authors did their due diligence during the rebuttal by providing additional results. Considering the manuscript, reviews, and discussion, there remain a number of concerns that withhold me from recommending acceptance at this time. First and foremost, the authors claim the adaptivity and flexibility of their method is unmatched, although the experiments do not support that the proposed method convincingly outperforms the baselines esp. in scenarios where it should. Second, originally raised by reviewer mmww, scalability of the method to larger numbers of variables is a major concern; the authors show equivalance of sparsity of the Jacobian to that of conditional independence, and as we all know, methods based on the latter often suffer in terms of scalability. Third, the manuscript at times leans towards overselling. It would certainly benefit from a clear discussion of the actual contributions (e.g. in light of existing work such as NoGAM) and in particular the limitations of the method (including the two previous points). All in all, I recommend the authors to improve their manuscript in these regards, and I am looking forward to seeing the improved version published at a later date.

Reviews: review: paper_summary: paper_summary: The paper extends theoretical results about causal discovery through score matching to encompass both linear and non-linear SCMs and lift the sufficiency assumption. The theoretical results relax the non-linearity assumption of Montagna et al 2023 by swapping it with the less restrictive one of restricted ANM from Peters et al 2009. As for the latent confounder detection a parallel with m-separation is drawn using results from Spantini et al 2018, to establish that the score will be non-zero in the presence of an active path. Following the theoretical results, an algorithm to estimate causal graphs from data is proposed and evaluated, generalizing the NoGAM algorithm of Montagna et al 2023, which only covers the non-linear case.
strength_weakness: strength_weakness: The paper motivation is basically the weakening of current assumptions for causal discovery methods. However, assumptions and benefits of the proposed methodology are not clearly specified. In line 74, the authors state that faithfulness is assumed (I believe, it is kind of hidden in the background notions). If that is the case, the method adopts the same assumption as FCI, plus ANM. So the proposed method is relaxing assumptions compared to CAM-UV, RCD and NoGAM, but adding onto FCI. Regarding the benefits, the alleged flexibility of the method to output DAGs, MECs, MAGs, PAGs, which should make it preferable to FCI, is merely touched upon in the contributions and the experiment section.  

Proposition 1 is a rather trivial application of the more general lemma in Spantini et al. and it does not specify the required faithfulness assumption to obtain the result from Eq. 6 in the paper. 

The experimental results show limited added value according to the one metric chosen (SHD) in a synthetic setting. They are not comprehensive enough, with no application to common (pseudo-)real benchmarks (e.g. from bnlearn). More experiments and more metrics are needed as, as it stands, the proposed method seems to add no real value compared to the baselines. Additionally, it is not clear from the experiments if it is really able to identify confounders. Breaking down precision and recall by mark would show this. FCI and a random baseline should be also added for reference.  

Experiments are conducted on data with at most 9 variables, and the scalability of the method is not shown nor discussed. 

The model used to estimate residuals is not discussed, nor the assumption that the chosen model fits the data adequately to correctly estimate residuals, and what is needed to assess this.
questions: questions: - Line 28: FCI weakness is highlighted to be that it outputs PAGs, does the proposed method improve on this at parity of assumptions? 

- Line 44: ANM is an assumption on the causal mechanism, is it not?

- Line 217: "the noise terms recentered by the latent causal effects" what does this mean? 

- Line 221: "N_i are assumed to be non-Gaussian when f_i is linear" This reduces to RCD, right? 

- Line 257: "the remaining arrowheads … are identified no better than in the equivalence class" what does this mean? 

- Line 279: "using hypothesis testing to find vanishing MSE..." how is this set up? 

- Line 280: how much added value does pruning have? E.g. CAM pruning can have huge impact https://arxiv.org/pdf/1310.1533 

- Alg. 1: do you start from a disconnected graph? Where do the "remaining" bidirected edges come from, do you reduce them prior to this? 

- Line 292: exogenous variables are selected at random, hence could be sink nodes with no consequence on the discovery process. Alternative and more meaningful testing strategies could be to select confounders (e.g. https://proceedings.neurips.cc/paper_files/paper/2021/file/144a3f71a03ab7c4f46f9656608efdb2-Paper.pdf) 

- Line 305: "presents better performance" I am not sure about this, they seem pretty aligned considering quartiles. Statistical tests to ascertain this are missing. 

- Line 315: some realistic scenarios are available through e.g. the bnlearn repository (https://www.bnlearn.com/bnrepository/). Any reason not to test on some of the data available there? 

Minor points and typos 

- Line 73: reference for d-separation is missing 

- Line 79: reference for MEC is missing 

- Line 130: have -> has 

- Line 147: "directed global Markov property": not introduced. I guess it is just the global Markov condition introduced before? 

- Footnote 1: provides -> provide 

- Line 162: "unobserved active paths" are not defined in Definition 5, m-separation.  

- Line 209: full stop missing 

- Line 236: the PAG -> a PAG 

- Line 267: "Appendix C.2" what should the reader expect to find in there? 

- Line 296: "We fix the […] level" significance
limitations: limitations: yes
paper_summary: paper_summary: The authors propose AdaScore, a method for causal discovery that generalizes previous work based on score matching for SCMs with possibly latent nodes. They combine connections of the score to conditional independence as well as to additive noise SCMs and show that a NoGAM-type procedure works to recover the direction of non-confounded edges of the corresponding partial ancestral graph.
strength_weakness: strength_weakness: The novelty of the paper lies primarily in the application of NoGAM to orient very specific edges in a partial ancestral graph (PAG) (the ancestral graph that represents the Markov equivalence class, in analogue to a CPDAG). This falls significantly short of the main contributions as described by the authors on l33--55. Specifically:

- The authors state that they show how constraints on the Jacobian of the score can be used as conditional independence testing. However, the extent to which this is done is only by noticing the equivalence between conditional independence and the corresponding zero in the Jacobian term (previously noted in [1,2]), without any formal analysis of the proposed t-test (Appendix C) as a statistical test of conditional independence (which happens to be a notoriously difficult test).

- The authors state that their identification results for additive noise models generalize the previous results obtained by previous works. In l193, the authors state "we remove the nonlinearity assumption (of [3]) and make the weaker hypothesis of a restricted additive noise model", but 1) this is a stronger assumption than additive noise, not a weaker one, and 2) the authors in [3] also consider the same restricted additive noise model. 

- The authors claim that AdaScore is able to handle a broad class of causal models (l54), but three out of four possible situations are direct applications of existing work. 1) Under no structural assumptions with or without latent confounders, AdaScore simply performs constraint-based causal discovery (FCI) using the conditional independence properties of the Jacobian of the score, a straightforward application of [1] also previously noticed in [2]. 2) Under an additive noise assumption, AdaScore is exactly equivalent to NoGAM. 3) Only under an additive noise assumption with hidden confounders, does AdaScore generalize NoGAM to orient unconfounded edges of the PAG returned by FCI, which may be very few of the discovered adjacencies.

### Other comments

- The experiments do not seem to suggest that AdaScore out performs other methods in any meaningful way---in fact, in Figure 1 a) AdaScore is completely equivalent to NoGAM, and is thus redundant. In Figure 1 b), where AdaScore should distinguish itself, it does not appear to be consistently better than CAM-UV. 

- Much of the paper (> 5 pages) is spent on directly describing previous works, NoGAM[3] and/or provide basic background on DAGs and MAGs. 

[1] Spantini et al., "Inference via low-dimensional couplings." JMLR 2018.
[2] Montagna et al., "Scalable causal discovery with score matching." CLeaR 2023.
[3] Montagna et al., "Causal discovery with score matching on additive models with arbitrary noise." CLeaR 2023.
questions: questions: - Is it correct to say that the unconfounded edges that are oriented by AdaScore are the purely _undirected_ edges of the PAG, and not the potentially bidirected ones?
limitations: limitations: The authors do not adequately discuss the limitations of their method---the limitations section in the appendix focuses purely on the empirical study. The authors claim that AdaScore is adaptive in the sense of being "less reliant on prior assumptions which are often untestable", but this is only in the sense that it performs different algorithms depending on user specification, which hardly constitutes one single unifying adaptive algorithm.
paper_summary: paper_summary: This paper present novel theoretical results to identify causal effects in restricted ANMs even in case of unobserved confounders.
strength_weakness: strength_weakness: I have a few remarks on improving the flow of the paper; however, even the first four points are not considered major issues.

- Even though condition 1 is a well-known result in the causality literature, I suggest explaining why that admits linear models and including some description of **restricted ANMs** in the main text (at least for me, it is not evident, especially since the condition lacks intuition). To be clear, even this point does not diminish the main contribution, which I see regarding the results for confounders.
- As **inducing paths** are an important concept for the main contributions, please _include it in the main text_ if space permits (suggestion: you can reduce spacing in $\texttt{itemize}$ by setting $\texttt{\\\\begin\\{itemize\\}[nolistsep]}$ )
- I could **not find the definition of an active path** (not even in Def. 5, where it is said to be defined); I presume it is a path that is not blocked, but it would be better to state this explicitly. Maybe it would even be better to use "a path that is not blocked" instead of introducing new terminology (this is the first time I encountered "active paths"; I could be wrong about this)
- The text is sometimes difficult to follow, due to heavy reliance on notation. I'd consider delegating the not crucial part to the appendix (potential candidates in 2.1) and using the remaining space to explain the main quantities better, especially the residuals (e.g., Eq. 12)


## Minor points
- please specify what you are calculating the expectation with respect to (using a bold E is also unconventional, though it's clear from the context it's an expectation)
- as the mathematical objects for d-/m-separation are distinguishable, you might consider dropping the superscript to simplify notation; also, I'd suggest adding whitespace after $\perp^d_\mathcal{G}$ and the like to make it easier for the reader to attribute the indices to $\perp$ and not the not on its right
- I could not find the definition for $\dot{\cup}$ 
- in the explanation of Prop 3., the wording makes it a bit hard to discern that you also provide intuition for the second part; it would help if you refer to _Part (ii)_ explicitly
- _Score matching through the roof_ in the title does not have added value for me, I'd consider rephrasing it to convey the message that "we propose score-based causal discovery methods for confounded restricted ANMs"
questions: questions: Some methods in the literature [1-4] use the Jacobian of a learned neural network instead of relying on the score for causal discovery/identifiability claims. How do you see the relation of your contribution to those Jacobian-based methods?


- [1] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural DAG learning. In 8th International Conference on Learning Representations, ICLR 2020

- [2] Lazar Atanackovic, Alexander Tong, Jason Hartford, Leo J. Lee, Bo Wang, and Yoshua Bengio. DynGFN: Bayesian Dynamic Causal Discovery using Generative Flow Networks, February 2023.

- [3] Patrik Reizinger, Yash Sharma, Matthias Bethge, Bernhard Scholkopf, Ferenc Huszar, and Wieland Brendel. Jacobian-based Causal Discovery with Nonlinear ICA. Transactions on Machine Learning Research, April 2023.

- [4] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A Linear NonGaussian Acyclic Model for Causal Discovery. Journal of Machine Learning Research, 7(10), 2006.
limitations: limitations: The authors provide an **honest comparison** of their method in the experiments, clearly stating its limitations compared to other methods.

Comments: title: Tables of p-values for sparse graphs
Tables of p-values for dense graphs
Statistical tests
Score increased 7->8
comment: We thank the reviewer for the active participation in the discussion. We are surprised that the reviewer specifies only now that all experiments with less than ten variables and in sparse settings are not meaningful in analyzing the algorithm performance, as this point was never raised in the review or in the authors-reviewers discussion, where the requests were limited to a better analysis of scalability (without even mentioning density as a discriminating point). Regarding the claim that these are the setting common in the literature, we point to CAM-UV paper (the work closest to ours), where they have experiments on at most nine variable, and one semi synthetic experiment on ten variables. Further, we notice that the analysis of the reviewer on the goodness of AdaScore compared to other methods is not based on statistically significant statements: e.g. CAM-UV is better than AdaScore with statistical significance in one experimental setting out of 32. Yet, the reviewer concludes that CAM-UV is better than AdaScore on inference on linear latent variables models, which encompasses 8 out of 32 of our experimental settings. As the statistical tests were an explicit and well-motivated request from the reviewer to be able to make statistically significant conclusions, we will stick to the conclusions supported by the statistical tests. 

In this spirit, we provide the summary table with the level of the test $0.05$, as asked by the reviewer:

|  | Alternative: Less | Alternative: Greater |
| --- | --- | --- |
| LiNGAM | 13/32 | 6/32 |
| RCD | 12/32 | 2/32 |
| NoGAM | 0/32 | 0/32 |
| CAM-UV | 1/32 | 1/32 |
| Random | 26/32 | 0/32 |
| Fully Random | 29/32 | 0/32 |

These results contrast with the analysis of the superiority of certain methods made by the reviewer. Specifically, we respond to specific points made by the reviewer and not supported by the empirical evidence:

- “Linear fully observable on par with random”: we have the following list of p-values for *alternative : less* on linear fully observable models comparing with random: 0.0002, 6.49e-05, 0.0005, 0.31, 5.41e-6, 5.41e-6, 5.41e-6, 5.41e-6. The only case in which p-value is higher than 0.05 is 0.31. The *alterantive : greater* pvalues instead are never below 0.05 threshold for random, and not even close: 1, 1, 1, 1, 0.99, 0.99, 0.99, 0.71. AdaScore is not on par with random.
- “linear latent: worse than camuv”: as already discussed, there is only one time in which *alternative :great* in comparison with CAM-UV, in the linear latent setting, goes below 0.05. The claim is not supported by empirical evidence

Further, the reviewer writes that the linear latent case is the focus of our work: we disagree, both the linear and nonlinear settings are equally important, as the point of the paper is indeed relaxing linearity and nonlinearity assumptions on the causal mechanisms, in the context of latent variable models. If we suggested that linear latent variable cases were the focus of our work, we kindly ask the reviewer to point to the specific parts of the paper. 

Finally, concerning the claim that contributions need to be clarified, we again point to L48, the paragraph starting with “our main contributions” written in bold. Concerning the claim that assumptions need to be clarified, the only specific criticism expressed on this point was relative to the positioning of faithfulness assumption, that could be more upfront. We already addressed this and remark on our offer to make a separate paragraph.
We are surprised by the request for finite samples guarantees on the score matching estimation, as this request was not in the original review. Overall it is clear that we disagree with the reviewer about this paper. In any case, learning theory results are far beyond the scope of the paper: the majority of causal discovery algorithms in the literature do not come with finite sample guarantees.

Concerning the use of the word adaptivity: we agree with the reviewer that we could adopt a better wording; we will adopt the phrasing used in the rebuttal to better express the benefits of our method.
Thanks for providing the test results. They go in the direction that I feared.
This is what I read from the tables provided (I focus on the dense graphs and 9 variables as, in my view, 10 variables graphs are a minimum of usefulness\practicality and usually a minimum in the literature too, $\alpha=0.05$):
- linear fully observable: on par with random and the others, worse than lingam (expected)
- linear latent: worse than camuv, rcd and lingam
- nonlinear full: better than rcd and lingam (expected)
- nonlinear latent: better than lingam (expected)
- all the rest are not significantly different

Overall, the results could be summarised as follows: Adascore performs better than lingam in the nonlinear case but worse in the linear one. **In the linear latent case, the focus of the proposed extension, it performs worse than all baselines apart the very akin nogam.** Hence, although the stronger guarantees, empirical performance shows that Adascore is mostly on par with prior methods, only better than a method specialised in linear SCM, and worse than others in the scenario where it provides the added guarantees. 

Overall, the discussion highlighted the following points:
- Assumptions need to be clarified
- Contributions need to be clarified
- Scalability issues need to be clarified
- Experimental results do not show much added value

I believe the above constitutes enough rework to need additional reviews. I therefore intend to keep my score. I thank the authors for the good discussion.
I would like to thank the authors for their response. However, the rebuttal does not sufficiently address my main concerns with the paper, and I will opt to retain my score at this time.

## On the main contribution of the paper.

> the fact that 3 out of 4 possible situations are extension of existing works, does not mean that the 4th one is not there.

As mentioned in the original review, I acknowledge the novelty of using score matching to orient edges in a mixed graph. I think I am in agreement with the authors that this is the main contribution of the work. In view of this, the paper does not do a thorough enough investigation to warrant publication at this time. The specific example I mentioned pertains to the lack of an analysis of the proposed t-test as a conditional independence test. 

>  We do not state that the Jacobian of the score can be used as conditional independence testing or that we propose a statistical test. We say that we use the Jacobian of the score as an alternative to conditional independence testing (L9-10, L41-42): while constraint-based methods (e.g., PC) usually rely on conditional independence testing (e.g. Zhang et al., 2011), instead, we consider sparsity of the Jacobian of the score. 

If using the Jacobian of the score as an alternative to a conditional independence test relieves it from being analyzed as a statistical test, then additional theory should be developed, e.g., adapting Thm. 3.2 from Zhu et al. on the theoretical properties of the discovered graph. This may not be a trivial addition, as Zhu et al. deals only with the non-linear Gaussian model. 

> we don’t merely note the equivalence. We prove it as a corollary of Spantini et al., (see the 1st point of the general response).

As reviewer `mmww` also points out, this is a rather trivial application of the result when assuming faithfulness.

## On adaptiveness of the algorithm

> We remark the adaptivity of AdaScore in the sense described in the paper, meaning that it can perform with theoretical guarantees on latent variable models, still identifying the causal direction when this is theoretically possible, which is not done by other methods (FCI, NoGAM, LiNGAM), or it is done under more restrictive assumptions (CAM-UV, RCD). 

This is closer to my understanding of the contribution of the paper, but  not how I understood the authors interpretation of "adaptivity". For example, in l261: "The main strength of our approach is its adaptivity with respect to structural assumptions: based on the user’s belief about the plausibility of several modeling assumptions on the data".

> We remark that this does not require any user specification. The user specification simply allows to turn off the search of latent variables, if this is thought as unnecessary.

I do not understand how this is not user specification. For me, an algorithm that is adaptive in the sense of the paper would need a routine to detect the presence of (significant) latent confounders, then automatically apply the mixed graph version. Otherwise, the AdaScore algorithm reads more as a computer application that includes previous methods rather than a unifying algorithm. 

### More minor points.

> Allowing for linear and nonlinear mechanisms is weaker than allowing for nonlinear mechanisms only.

My apologies, thank you for the clarification.

> Point 3) of the review says that the only novel part “may only apply to few edges.” 

I acknowledge that the "few edges" comment was subjective. I recognize that this may have come off as dismissive and I apologize for that. This was however not a significant point for my score.
For completeness, we present the table of p-values of all the MW tests on experiments relative to __sparse graphs__. Please notice that adding the random and fully random baselines to the experiments changed the random state of our running relative to Fig. 3 in the main paper. Results are qualitatively the same, but we observe a few mismatches due to the variance resulting from different random states.

# Alternative: less

__linear Fully observable model__
|              |           3 |           5 |           7 |           9 |
|:-------------|------------:|------------:|------------:|------------:|
| camuv        | 0.5147    | 0.6020    | 0.2644    | 0.8236    |
| nogam        | 0.5147    | 0.4852    | 0.1575      | 0.05256   |
| rcd          | 0.3696    | 0.7820    | 0.0216   | 0.6578    |
| lingam       | 0.5147    | 0.8762    | 0.8425      | 0.9999    |
| random       | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |
| fully_random | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |


__nonlinear Fully observable model__
|              |           3 |           5 |           7 |           9 |
|:-------------|------------:|------------:|------------:|------------:|
| camuv        | 0.2893    | 0.2893    | 0.4267    | 0.09515   |
| nogam        | 0.4852    | 0.2179    | 0.1575      | 0.3152    |
| rcd          | 5.41254e-06 | 5.41254e-06 | 1.08251e-05 | 0.0003  |
| lingam       | 0.0034  | 0.0034  | 0.0177   | 0.0446   |
| random       | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |
| fully_random | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |


__linear Latent variables model__
|              |           3 |           5 |           7 |           9 |
|:-------------|------------:|------------:|------------:|------------:|
| camuv        | 0.5147    | 0.6847    | 0.4852   | 0.8912    |
| nogam        | 0.4852    | 0.6303    | 0.5732    | 0.2644    |
| rcd          | 0.3696    | 0.3979    | 0.1237    | 0.3979    |
| lingam       | 0.6578    | 0.8425      | 0.7593    | 0.9942    |
| random       | 2.16502e-05 | 2.16502e-05 | 0.000162376 | 5.41254e-06 |
| fully_random | 1.08251e-05 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |


__nonlinear Latent variables model__
|              |          3 |           5 |           7 |           9 |
|:-------------|-----------:|------------:|------------:|------------:|
| camuv        | 0.3421   | 0.1763    | 0.5441    | 0.2406    |
| nogam        | 0.1965   | 0.5147    | 0.217936    | 0.0987    |
| rcd          | 0.0044 | 0.0025  | 0.0525  | 0.0177   |
| lingam       | 0.0715  | 0.0116   | 0.1399    | 0.0057  |
| random       | 0.0003 | 1.08251e-05 | 6.49505e-05 | 0.0001 |
| fully_random | 0.0003 | 5.41254e-06 | 2.16502e-05 | 2.16502e-05 |

---
&nbsp;

# Alternative: greater

__linear Fully observable model__
|              |        3 |        5 |        7 |           9 |
|:-------------|---------:|---------:|---------:|------------:|
| camuv        | 0.5147 | 0.4267 | 0.7593 | 0.1965    |
| nogam        | 0.5147 | 0.5441 | 0.8762 | 0.9553    |
| rcd          | 0.6578 | 0.2644 | 0.9855 | 0.3979    |
| lingam       | 0.5147 | 0.1399 | 0.1763 | 0.0001 |
| random       | 1        | 1        | 1        | 1           |
| fully_random | 1        | 1        | 1        | 1           |


__nonlinear Fully observable model__
|              |        3 |        5 |        7 |        9 |
|:-------------|---------:|---------:|---------:|---------:|
| camuv        | 0.7355 | 0.7355 | 0.6020 | 0.9172 |
| nogam        | 0.5441 | 0.8236 | 0.8600 | 0.7106 |
| rcd          | 1        | 1        | 0.9999 | 0.9997 |
| lingam       | 0.9974 | 0.9974 | 0.9855 | 0.9684 |
| random       | 1        | 1        | 1        | 1        |
| fully_random | 1        | 1        | 1        | 1        |


__linear Latent variables model__
|              |        3 |        5 |        7 |          9 |
|:-------------|---------:|---------:|---------:|-----------:|
| camuv        | 0.5147 | 0.3696 | 0.5441 | 0.1399   |
| nogam        | 0.5441 | 0.4267 | 0.4558 | 0.7593   |
| rcd          | 0.6578 | 0.6578 | 0.9048 | 0.6303   |
| lingam       | 0.3696 | 0.1965 | 0.2893 | 0.0092 |
| random       | 0.9999 | 0.9999 | 0.9999 | 1          |
| fully_random | 0.9999 | 1        | 1        | 1          |


__nonlinear Latent variables model__
|              |        3 |        5 |        7 |        9 |
|:-------------|---------:|---------:|---------:|---------:|
| camuv        | 0.6847 | 0.8425   | 0.5147 | 0.8034 |
| nogam        | 0.8236 | 0.5147 | 0.8236 | 0.9048 |
| rcd          | 0.9965 | 0.998  | 0.9623 | 0.9855 |
| lingam       | 0.9384 | 0.9907 | 0.8912 | 0.9965 |
| random       | 0.9997 | 0.9999 | 0.9999 | 0.9999 |
| fully_random | 0.9997 | 1        | 0.9999 | 0.9999 |
For completeness, we present the table of p-values of all the MW tests on experiments relative to __dense graphs__. Please notice that adding the random and fully random baselines to the experiments changed the random state of our running relative to Fig. 1 in the main paper. Results are qualitatively the same, but we observe a few mismatches due to the variance resulting from different random states.

# Alternative: less

__Linear fully observable model__
|              |           3 |           5 |           7 |         9 |
|:-------------|------------:|------------:|------------:|----------:|
| camuv        | 0.657895    | 0.544101    | 0.485256    | 0.782064  |
| nogam        | 0.657895    | 0.342105    | 0.139931    | 0.455899  |
| rcd          | 0.514744    | 0.782064    | 0.514744    | 0.8425    |
| lingam       | 0.657895    | 0.803476    | 0.999756    | 1         |
| random       | 0.000243564 | 6.49505e-05 | 0.000525017 | 0.315264  |
| fully_random | 1.08251e-05 | 1.08251e-05 | 5.41254e-06 | 0.0216285 |

__nonlinear Fully observable model__
|              |           3 |           5 |           7 |           9 |
|:-------------|------------:|------------:|------------:|------------:|
| camuv        | 0.139931    | 0.240625    | 0.0177315   | 0.098781    |
| nogam        | 0.397968    | 0.369682    | 0.602032    | 0.630318    |
| rcd          | 3.78878e-05 | 0.000243564 | 2.16502e-05 | 5.41254e-06 |
| lingam       | 0.0177315   | 0.00927169  | 0.0315064   | 0.00446535  |
| random       | 0.000102838 | 1.08251e-05 | 1.08251e-05 | 1.08251e-05 |
| fully_random | 6.49505e-05 | 5.41254e-06 | 5.41254e-06 | 5.41254e-06 |

__linear Latent variables model__
|              |         3 |           5 |        7 |         9 |
|:-------------|----------:|------------:|---------:|----------:|
| camuv        | 0.485256  | 0.823659    | 0.544101 | 0.982269  |
| nogam        | 0.264424  | 0.782064    | 0.289371 | 0.1575    |
| rcd          | 0.426714  | 0.630318    | 0.735576 | 0.985597  |
| lingam       | 0.455899  | 0.139931    | 0.982269 | 0.998955  |
| random       | 0.0715701 | 0.000525017 | 0.139931 | 0.0715701 |
| fully_random | 0.037628  | 0.000525017 | 0.176341 | 0.426714  |

__nonlinear Latent variables model__
|              |        3 |          5 |          7 |          9 |
|:-------------|---------:|-----------:|-----------:|-----------:|
| camuv        | 0.735576 | 0.426714   | 0.684736   | 0.426714   |
| nogam        | 0.630318 | 0.217936   | 0.630318   | 0.455899   |
| rcd          | 0.973787 | 0.0525612  | 0.759375   | 0.426714   |
| lingam       | 0.289371 | 0.00574812 | 0.0216285  | 0.00927169 |
| random       | 0.123725 | 0.0019431  | 0.0116153  | 0.00342073 |
| fully_random | 0.176341 | 0.00036264 | 0.00927169 | 0.00342073 |

---
&nbsp;

# Alternative: greater

__linear Fully observable model__
|              |        3 |        5 |           7 |           9 |
|:-------------|---------:|---------:|------------:|------------:|
| camuv        | 0.369682 | 0.514744 | 0.544101    | 0.240625    |
| nogam        | 0.369682 | 0.684736 | 0.876275    | 0.602032    |
| rcd          | 0.514744 | 0.264424 | 0.544101    | 0.196524    |
| lingam       | 0.369682 | 0.240625 | 0.000525017 | 5.41254e-06 |
| random       | 0.999897 | 0.999978 | 0.999637    | 0.710629    |
| fully_random | 1        | 0.999995 | 1           | 0.982269    |


__nonlinear Fully observable model__
|              |        3 |        5 |        7 |        9 |
|:-------------|---------:|---------:|---------:|---------:|
| camuv        | 0.891219 | 0.782064 | 0.988385 | 0.917253 |
| nogam        | 0.630318 | 0.657895 | 0.455899 | 0.426714 |
| rcd          | 0.999989 | 0.999838 | 0.999995 | 1        |
| lingam       | 0.988385 | 0.992655 | 0.973787 | 0.997402 |
| random       | 0.999935 | 0.999995 | 1        | 0.999995 |
| fully_random | 0.999962 | 1        | 1        | 1        |


__linear Latent variables model__
|              |        3 |        5 |         7 |         9 |
|:-------------|---------:|---------:|----------:|----------:|
| camuv        | 0.544101 | 0.196524 | 0.514744  | 0.026213  |
| nogam        | 0.759375 | 0.264424 | 0.759375  | 0.876275  |
| rcd          | 0.602032 | 0.426714 | 0.315264  | 0.0177315 |
| lingam       | 0.573286 | 0.891219 | 0.0216285 | 0.0019431 |
| random       | 0.938497 | 0.999637 | 0.891219  | 0.938497  |
| fully_random | 0.968494 | 0.999756 | 0.860069  | 0.602032  |


__nonlinear Latent variables model__
|              |        3 |        5 |        7 |        9 |
|:-------------|---------:|---------:|---------:|---------:|
| camuv        | 0.315264 | 0.602032 | 0.342105 | 0.630318 |
| nogam        | 0.426714 | 0.823659 | 0.426714 | 0.602032 |
| rcd          | 0.037628 | 0.955395 | 0.289371 | 0.630318 |
| lingam       | 0.759375 | 0.996579 | 0.985597 | 0.992655 |
| random       | 0.904842 | 0.99856  | 0.992655 | 0.998057 |
| fully_random | 0.860069 | 0.999756 | 0.994252 | 0.997402 |
We ran WM statistical tests on the SHD relative to the 32 experimental settings presented in Figure 1 (dense graphs) and Figure 3 (sparse graphs) of the main paper. We summarize the statistical tests outcomes in the following table: in the “*Alternative: less”* column we provide the number of times we reject the null in favor of the alternative that AdaScore SHD is stochastically lower (i.e. AdaScore "better"), in the “*Alternative: greater”* column we provide the number of times we reject the null in favor of the alternative that AdaScore SHD is stochastically greater (i.e. AdaScore "worse"). As level of the test to reject the null, we consider 0.1. 

&nbsp;

|  | Alternative: Less | Alternative: Greater |
| --- | --- | --- |
| LiNGAM | 14/32 | 6/32 |
| RCD | 14/32 | 2/32 |
| NoGAM | 2/32 | 0/32 |
| CAM-UV | 3/32 | 1/32 |
| Random | 29/32 | 0/32 |
| Fully Random | 29/32 | 0/32 |
&nbsp;

These __results suggest the following interpretations__:

- None of the alternative methods is preferable to AdaScore
- In absence of prior knowledge about the SCM underlying data generation, AdaScore is preferable in terms of SHD accuracy to LiNGAM, RCD, and the random baselines
- In absence of prior knowledge about the SCM underlying data generation, we can not say much in comparison to NoGAM and CAM-UV.

These outcomes match our findings in the paper, where we summarise that AdaScore “accuracy is often comparable and sometimes better than competitors.” (L 311-312); findings that are now more substantiated, for which we credit the reviewer suggesting a better experimental analysis. Note that these results are achieved while having the least strict assumptions among all of the benchmarked methods. __We propose to add these results to the main paper__, to substantiate our claims on empirical performance. Additionally, we remark on our intention to add all the additional experiments already discussed with the reviewer and partially presented in the PDF of the rebuttal.

For the sake of completeness, we present the table of p-values of all the MW tests in separate comments (due to space limits). As these contain valuable information, we propose to include barplots with binned range to visualize the different p-values associated to different algorithms in the main paper, such that results interpretation is not dependent on the choise on the level of the test.
Thanks for the availability and commitment to provide clear and robust results for the evaluation of your work.

If you are willing to make the extra assumption (of equal distributions up to a shift in location) in order to interpret the rank WM test as a difference in medians, the same assumption would make a one-sided t-test on the difference in means also valid. In fact, WM is usually used to test differences for asymmetric distributions, which, if you assume gaussianity, do not apply. A t-test would also consider dispersion around the mean, rather than the rank only. 

In any case, I would find a WM test satisfactory if you already have it setup. I am not trying to make the authors' life difficult, and again, I appreciate the efforts.
We thank the reviewer for promptly engaging in the discussion.

- **Q3.** Fig 5: we apologise for the confusion. We will add a similar plot with the elapsed time on the y-axis.
- **Q6.** We did not intend that the reviewer needs to prove results about the goodness of Adascore, and we apologize for this confusion: our quantitative statements were aimed at rebutting that AdaScore “*is almost always at par or worse than other baselines”,* as the empirical evidence does not support this. To help the analysis of our empirical results, we are willing to provide statistical tests in the time frame of the authors reviewers discussion. In particular, we would consider Whitney Mann-U test, testing for each algorithm and node size
    1. the null-hypothesis that the distributions are identical against the alternative that adascore has a stochastically lower distribution (i.e. lower median, e.g. assuming the average SHD of adascore and the competitor algorithm and adascore have the same distribution up to shift in location, in our case ~ Gaussian by CLT)
    2. The null-hypothesis that the distributions are identical against the alternative that adascore has a stochastically greater distribution (i.e. greater median)
    
    If the reviewer finds this a satisfactory approach,  we will do our best to provide these results in time to discuss them with the reviewer. Otherwise, we kindly ask the reviewer to specify what statistical tests they had in mind for the comparison.
Q2 

I agree that my proposed strategy is less straightforward than I anticipated given the different meanings of the mixed graph output of your method. I frankly cannot think of a solution right now but I take your point that you did what you could in the timeframe. 


Q3 

I am well aware of the space limitation, I was not asking the authors to move C.3. I was pointing out that the current content of the paper regarding scalability is not satisfactory in my view. I still believe that. In fact, the table you provide about elapsed time (thanks) already shows that the actual performance has a very exponential trend. I was just asking for clarity.  

Fig 5 in the paper shows shd. Number of samples is different from numeber of nodes. 

 



Q6 (Line 305) 

I am afraid but the burden of proof falls on you here, not me. I am not claiming anything, I am just supposing, and trying to interpret some unclear results. Again, what I am trying to elicit here is clarity, given that the authors do not provide readers with the means (e.g. statistical tests) to validate their claims.  

The author's claim that "CAM-UV better than AdaScore 4 out of 16 times, NoGAM 2 out of 16, RCD 4 out 16, LiNGAM 4 out of 16" is unsubstantiated until they present tests to back up it up. The overlap in distributions shown in Figure 1 make me suppose that the statistical differences between pairs of methods will be very minimal. But again, it is a supposition since I do not have the means to judge.
We thank the reviewer for dedicating time to our rebuttal. 

---

**Q1** “*Indeed I would create CI around the numbers you provided.”* We agree with the reviewer relative to the need to produce confidence intervals: we commit to adding them in the camera-ready.

---

**Q2** 
- *”the way you explained in the second bullet of the general rebuttal is already much clearer”* We will adopt this wording and add explicit assumption paragraphs. 

- *“The way I would go about it is to transform the output mixed graphs from ADAscore to a PAG, then compare to FCI.”* We find the required experiment not as straightforward as it seems. The output of AdaScore a priori only makes statements about direct causal connections but not about d-separations (or m-separations). E.g. suppose our algorithm outputs $X \leftrightarrow Y \leftrightarrow Z$. This means we cannot orient the edge between X and Y and Y and Z, which could mean that either there is a hidden mediator between X and Y (and the edge could go in either way, i.e. $X \rightarrow L \rightarrow Y$ or $X \leftarrow L \leftarrow Y$) or a hidden confounder. In the former case, we could have both, $X \perp Z | Y$ and $X \not\perp Z | Y$. Therefore the PAG w.r.t. our output is not well defined. Thus, to test “the other half” of our work a more sensible approach is via experiments as in Figure 1 of the paper.

---

**Q3.** “*Re scalability …”*

- The pointer to Appendix C.3 is indeed to better discuss the scalability, as the content of the appendix would not fit in the main paper page limit.
- On scalability experiments
    1. Fig 5 of the paper already shows AdaScore scalability with the number of samples.
    2. We propose to add more experiments on AdaScore scalability with the number of nodes: we present preliminary results showing that AdaScore scales much better than CAM-UV, the method with best theoretical guarantees previous to ours. The values in the table are mean +- std in seconds.
    
    |  | 3 nodes | 5 nodes | 7 nodes | 9 nodes |
    | --- | --- | --- | --- | --- |
    | adascore | 1.555 ± 0.274 | 12.582 ± 6.568 | 43.596 ± 29.867 | 133.950 ± 81.151 |
    | camuv | 3.999 ± 2.340 | 17.703 ± 3.894 | 80.196 ± 31.040 | 198.977 ± 69.306 |
    | nogam | 2.374 ± 0.129 | 6.127 ± 0.631 | 12.042 ± 0.608 | 20.407 ± 0.796 |
    | rcd | 0.409 ± 0.330 | 2.507 ± 0.960 | 10.720 ± 3.774 | 22.798 ± 3.641 |
    | lingam | 0.015 ± 0.0 | 0.051 ± 0.001 | 0.124 ± 0.001 | 0.251 ± 0.001 |

---
**Q4 (Line 28)** *“great, now it is clear, thanks. Again, from a reader's perspective, the way the motivation is presented raises expectations that the paper does not fulfill.”* Reviewer originally asked whether we could identify more than FCI, at parity of assumptions: as we discussed, this is theoretically impossible (see Glymour et al.) Thus, stating “*The FCI algorithm [11] can only return an equivalence class from the data. Appealing to additional restrictions ensures the identifiability of some direct causal effects in the presence of latent variables*” as done in L28-30 does not raise unfulfilled expectations, as this is exactly what we do in our paper.

---
**Q5 (Line 44)** “*in my view, the paper reads like ANM is not an assumption worth noting.”*
 L44 we write “*we prove that the score function identifies the causal direction of ANMs, with minimal assumptions on the causal mechanisms”,* the *minimal* part refers to additivity of the noise: this is a minimal requirement to achieve identifiability in the light that most of existing works place linearity/nonlinearity assumptions on the mechanisms.

We will remove the word minimal, stating “*we prove that the score function identifies the causal direction under the assumption of additive noise models, without further requirements on the causal mechanisms.”* We thank the reviewer to help us clarifying this point. 

---

**Q6 (Line 305)**
- *“for linear systems with confounders it is basically the same as random and lingam is much better”* Out of the 16 experimental settings of Fig. 1 of the PDF rebuttal, AdaScore is comparable to random only in one or two - large-scale linear. LiNGAM is comparable to random in at least 5 out of 16 (see nonlinear).
- *“it is almost always at par or worse than other baselines”* Considering median performance (please, recall that we already committed to adding statistical tests in the main paper) we see CAM-UV better than AdaScore 4 out of 16 times, NoGAM 2 out of 16, RCD 4 out 16, LiNGAM 4 out of 16.

Reviewer's claim that other methods are better than AdaScore is not supported by empirical evidence, and it does not account that their performance is achieved via stricter assumptions. Additionally, we remark on the following points that emerged in the discussion agreed upon both by us and the reviewer

1. Conclusions from experiments must be drawn with care
2. Our baseline has better theoretical guarantees, which are very relevant as they are valid beyond the scope of our experiments, i.e. in the real world.
Many thanks for the clarifications and the extra experiments provided. Responding in order: 

- Thanks for providing the results on some of the bnlearn datasets. I agree on having to take them with a pinch of salt. Indeed I would create CI around the numbers you provided, by changing the seed of the BN DGP. By the looks of the plot, I believe that a lot of them could not be significantly different from each other. 

- I apologise if my comment about the assumption being "kind of hidden" came across the wrong way. I was in no way accusing. I was just pointing out that, given the importance of assumptions in your claim about the added value of your proposed method, the organisation and wording of the paper do not make clear enough their importance, strictness or benefits. In my view, the way you explained in the second bullet of the general rebuttal is already much clearer and effective, at a high level, than the paper. As for the comparison with FCI, the skeleton accuracy is good. The way I would go about it is to transform the output mixed graphs from ADAscore to a PAG, then compare to FCI. Otherwise you only compare less than half of the work done by the algorithm. 

- The point to stress is not that Eq 6 does not use faithfulness, it's that Proposition 1 does. 

- Re scalability, thanks to the pointer to the appendix, I had indeed seen that. I still think that the discussion is not satisfactory. The main paper just reads: "This way, we get an algorithm that is polynomial in the best case (Appendix C.3)." and, to me, that is not a great way to discuss potential limitations. Why only 9 variables? I would show a comparison of elapsed/computing time to back up your claim that the proposed algorithm is practical and with superior guarantees. 

- Line 28: great, now it is clear, thanks. Again, from a reader perspective, the way the motivation is presented raises expectations that the paper does not fulfill. 

- Line 44: in my view, the paper reads like ANM is not an assumption worth noting. To go back to my previous point about "kind of hidden", assumptions should be better organised and elaborated upon. 

- Line 279: Thanks for the pointer. I had not read that since you only reference it when talking about the pruning. I would separate the two discussions, title them more clearly and reference them accordingly in the main paper. 

- Line 280: interesting, the reference I provided show a large impact of pruning. 

- Line 305 and extra experiments provided: I still think that the experiments show that the method is not really good at doing what it is supposed to do better than the others algorithms: for linear systems with confounders it is basically the same as random and lingam is much better (though, of course, it is more specialised). For the rest of the scenarios it is almost always at par or worse than other baselines. Confounder detection is quite poor, as shown in Figure 4 in the additional results.
Thank you for your detailed and thoughtful response. My concerns are addressed, I increase my score 7->8.