Hybrid Top-Down Global Causal Discovery with Local Search for Linear and Nonlinear Additive Noise Models
Authors: Sujai Hiremath
Jacqueline R. M. A. Maasch
Mengxiao Gao
Promit Ghosal
Kyra Gan
Venue: NeurIPS 2024 poster (2024)
Keywords: global causal discovery
additive noise model
local structure
Primary area: probabilistic_methods
TL;DR: We present a novel hybrid method for causal discovery that leverages local causal relationships in SEMs to construct a compact hierarchical topological sort, followed by a novel nonparametric constraint-based method for efficient edge discovery.

Abstract: Learning the unique directed acyclic graph corresponding to an unknown causal model is a challenging task. Methods based on functional causal models can identify a unique graph, but either suffer from the curse of dimensionality or impose strong parametric assumptions. To address these challenges, we propose a novel hybrid approach for global causal discovery in observational data that leverages local causal substructures. We first present a topological sorting algorithm that leverages ancestral relationships in linear structural causal models to establish a compact top-down hierarchical ordering, encoding more causal information than linear orderings produced by existing methods. We demonstrate that this approach generalizes to nonlinear settings with arbitrary noise. We then introduce a nonparametric constraint-based algorithm that prunes spurious edges by searching for local conditioning sets, achieving greater accuracy than current methods. We provide theoretical guarantees for correctness and worst-case polynomial time complexities, with empirical validation on synthetic data.

Hybrid Top-Down Global Causal Discovery with Local Search for Linear and Nonlinear Additive Noise Models
Sujai Hiremath
Jacqueline Maasch Cornell Tech
Mengxiao Gao
Cornell Tech
Tsinghua University
sh2563@cormel1.com edu
jam887@corne11 edu
gaomx210mails. www.gtut..ncncmm. edu. cn
Promit Ghosal University Chicago
Kyra Gan
Cornell Tech
Promit@uchicago edu
kyragan@ccorelll.dde edu
Abstract
Learning the unique directed acyclic graph corresponding to an unknown causal model challenging task. Methods based on functional causal models identify a unique graph, but either suffer from the curse of dimensionality or impose strong parametric assumptions. To address these challenges, propose novel hybrid approach for global causal discovery in observational data that leverages local causal substructures. We first present topological sorting algorithm that leverages ancestral relationships in linear structural causal models to establish compact top-down hierarchical ordering, encoding more causal information than linear orderings produced by existing methods. We demonstrate that this approach generalizes nonlinear settings with arbitrary noise. We then introduce a nonparametric bonstraint-based algorithm that prunes spurious edges by searching for local conditioning sets, achieving greater accuracy than current methods. We provide theoretical guarantees for correctness and worst-case polynomial time complexities, with empirical validation on synthetic data.
Introduction
Causal graphical models compactly represent the data generating processes (DGP) of complex systems, including physical, biological, and social domains. Access to the true causal graph or its substructures can offer mechanistic insights [31, 13] and enable downstream causal inference. including effect estimation [11, 34, 2, 8, 16, 35]. In practice, the true causal graph is often unknown, and can challenging to assume using domain knowledge. such -ootteedd-owwedeeee settings, we can instead rely on causal discovery algorithms that learn the causal graph from observational data in a principled, automated manner [41, 7].
Traditional approaches to causal discovery infer causal relationships either through conditional independence relations (PC [40]) goodness-of-fit measures (GES [3], GRaSP [12]). While these discovery methods can provide flexibility by not requiring assumptions over the functional form of the DGP, they are generally worst-case exponential in time complexity and learn Markov equivalence
classes (MEC) rather than unique directed acyclic graphs (DAGs) [19]. Therefore, additional modeling assumptions are often necessary for time-efficient and accurate global discovery.
Certain parametric assumptions can enable recovery of the unique ground truth DAG, e.g., assuming particular functional causal model (FCM) [48]. Under the additive noise model (ANM), obtain unique identifiability by assuming linear causal mechanisms with non -on-Gassia noise distributions
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

[3] or nonlinear causal functions with arbitrary noise [10]. Under the independent additive noise assumption, the causal parents of variable are statistically independent of its noise term. For this class of models, discovery often entails regressing the variable of interest against its hypothesized parent set and testing for marginal independence between this set and the residual term [25].
Current FCM approaches to global causal discovery trade off between two main issues, suffering from either 1) strong parametric assumptions over the noise or functional form (or both) 2) the use of high-dimensional nonparametric regressions, which require large sample sizes for reliable estimation and do not scale to large graphs. In addition, current for causal discovery in sparse causal graphs. setting that characterizes many high- high-dimensional applications (e.g., analysis of genetic data in healthcare applications) [6, 47, 46, 15].
methods are ill-suited
Contributions We propose hybrid causal discovery approach to graph learning that combines functional causal modeling with constraint-based discovery. We depart from previous methods by characterizing conditions that allow us to search for and exploit local, rather than global, causal relationships between vertices. These local relationships stem from root vertices: this motivates a top-down, rather than bottom-up, approach. Thus, we learn the topological sort and discover true edges starting from the roots rather than the leaves, as in existing methods [25, 30, 21J. This approach leverages sparsity in both the ordering phase and edge discovery phase to reduce the size of conditioning sets, as well as the number of digenssinni regressions. We summarize our major
contributions as follows:
We introduce topological ordering algorithm LHTS for linear non-Gaussian ANMs that exploits local ancestor-descendent relationships to obtain a compact hierarchical sort.
We introduce topological ordering algorithm NHTS for nonlinear ANMs that exploits local parent-child relationships to run fewer high-dimensional regressions than traditional methods, achieving lower sample complexity.
We introduce a constraint- a discovered topological ordering, leveraging local properties of causal edges to use smaller conditioning sets than traditional sparse regression techniques.
algorithm ED that nonparametrically prunes spurious edges from
We achieve accurate causal discovery in synthetic data, outperforming baseline methods.
Organization After describing the preliminaries in Section 2, we introduce the linear problem setting in Section 3, establishing the connection between ancestral relationships and causal active paths and introducing a linear hierarchical topological sorting algorithm (LHTS). Next, we extend our method to the nonlinear setting in Section 4 by establishing the connection between parental relationships and active causal paths, introducing a nonlinear hierarchical topological sorting algorithm (NHTS). In section 5, we establish a sufficient conditioning set for determining edge relations and introduce an efficient edge discovery algorithm (ED). We then test LHTS, NHTS and ED in synthetic experiments in Section 6. To conclude, we discuss future work that might generalize our approach to full ANMs.
Related Work Our work is related to two kinds of discovery methods that explicitly leverage the topological structure of DAGs: 1) permutation-basee approaches, and 2) FCM-based approaches.
The original permutation- based approach SP [27] searches over the space of variable orderings find
permutations that induce DAGs with minimal edge counts. Authors in [39] introduce greedy variants of SP (such as GSP) that maintain asymptotic consistency; GRaSP [12] relaxes the assumptions of prior methods to obtain improvements in accuracy. These methods highlight the importance of using permutations for efficient causal discovery, but generally suffer from the need to bound search
runtime with heuristics, poor sample efficiency in high dimensional settings, and are unable to recover a unique topological ordering or DAG ([22))
On the other hand, the recent stream of FCM-based approaches decompose graph learning into two phases: 1) learning the topological sort, i.e., inferring a causal ordering of the variables; and 2) edge discovery, i.e., identifying edges consistent with the causal ordering [38, 25, 1, 30, 21, 32, 20).
The literature topological ordering algorithms for ANMs is organized along the types of parametric assumptions made on both the functional forms and noise distributions of the underlying DGP. Early approaches like ICA-LiNGAM [37] and DirectLiNGAM [38] focus on learning DAGs generated by linear functions and non-Gaussian noise terms. Recent work leverages score matching to obtain the causal ordering in settings with nonlinear functions and Gaussian noise: SCORE [30] and DAS [21] exploit particular variance properties, while DiffAN estimates the score function with a diffusion model [32]. NoGAM [20] generalizes the score moore-matching procedure of SCORE to nonlinear causal

---

mechanisms with arbitrary noise distributions. RESIT [25] leverages residual independence results in nonlinear ANMs to identify topological orderings when the noise distribution is arbitrary. NoGAM and RESIT both rely on high-dimensional nonparametric regression.
Once a topologiaal ordering is obtained, spurious edges are pruned. Works that agnostic the distribution of noise often use parametric approach, implementing either a form of sparse regression (e.g., Lasso regression [18]) or a version of additive hypothesis testing with generalized additive models (GAMs) [17] (e.g., CAM-pruning [1)). RESIT [25] provides another alternative edge pruning procedure for nonlinear ANMs, relying again on high-dimensional nonparametric regression.
2
Preliminaries
We focus on structural causal models (SCMs) represented as DAGs. These graphs describe the causal relationships between variables, where an edge r; + Ij implies that x; has direct causal influence æj. Let G = (V, E) be a DAG on = d vertices, where E represents directed edges. To define pairwise relationships between vertices, we let Ch(zs) denote the children of x; such that Ty E Ch(z;) if and only if x; + z;, and Pa(æ:) denote the parents of such that z. E Pa(a:) if and only if z2 + zg. Similarly, let An(x,) denote the ancestors of F; such that z; E An(x;) if and Ti, and De(xz,) denote the descendants of 2; such that I E De(x;) ) if and only if there exists a directed path T; -++ z.. Vertices can be classified based on the totality of their pairwise relationships: x; is root if and only if Pa(x;) = 0, a leaf if and only if Ch(x,) = 0. an isolated vertex if v; is both a root and a leaf, and an intermediate vertex otherwise. See illustrative DAG in Figure 1. Vertices
only if there exists directed path 2j
also be classified in terms of triadic relationships: m, is a confounder of æy. 2y if and only if z; E An(z,) n An(x)) ); a mediator of x to 2y. if and only if cry E De(rr) n An(r):;; and a collider between 2y and 37 if and only if x;; E De(w;) O De(rrg).
Undirected paths that transmit causal information between vertices Fj, Jy can be differentiated into frontdoor and backdoor paths [42]. A frontdoor path is directed path z, --+
I4. that starts Ik
with an edge out of æj. and ends with an edge into zg. A backdoor path is a path zy +-- that starts with edge into xj, and ends with edge into Zg. Paths that start and end with an edge out of z and Ty (Iig --+ 6-- 2y) do not transmit causal information between Ij. Zg.
Paths between two vertices are further classified, relative to a vertex set Z, as either active or inactive [42]. A path between vertices zy,, 2g is active relative
every node on the path is active relative to Z. Vertex Iy on a path is active if one of the following holds: 1) x is not a collider and x; & Z. 2) x; is collider and x; E Z, 3) xt; is a collider and x; & Z. but De(x;)) n 7. * 0. An inactive path is simply a path that is not active. Following convention, throughout the rest of the paper we will describe causal paths as active or inactive with respect to 0 unless otherwise specified.
Definition 2.1 (Topological Orderings). Consider given DAG G = (V, F).. A topological sort (linear order) is mapping π V + 1.. V 1}, such that if z; E Pa(x,)) then I; appears before Ij the sort π: π (z:) T TTæ ). A hierarchical sort (between a partial and linear order) is | 1}, such that if Pa(x;) = 0, then π (I.) = O, and if Pa(az) & 0,
mapping TL V + then Ty (x,) equals the maximum length of the longest directed path from each root vertex to Ty, i.e. T (I:) = + max (z,) Iz E P-(x:))-
We note that the hierarchical sort is unique, and that it coincides with a topological sort when the number of layers equals [V], i.e., the DAG is complete.
Definition 2.2 (ANMs). ANMs [9] are a popular general class of SCMs defined over a DAG G with
xt (=((a(a))) + Ei, Vx; E V,
where fis are functtons and EiS are independent arbitrary noise distributions.
This model implicitly assumes the causal Markov condition and acyclicity; we adopt the aforemen- tioned assumptions, as well as faithfulness [41].
3
Linear Setting
We first restrict our attention to ANMs that feature only linear causal functions f, known as Linear Non-Gaussian Acyclic causal Models (LiNGAMs). Following [38], note that LiNGAM can
3

---

T3
Figure 1: Illustraivve DAG, where 21 is root, 23 is leaf, 23 E Ch(r2), E De(r1).
API
AP2
AP3
AP4
Zi
Jj
Ik
24
I
F4
Ij
Ft
Figure 2: Enumeration of active causal path relation types between a pair of nodes I; and rj. Dashed arrows indicate ancestorship.
be represented as a d adjacency matrix {bıy]. }. where biy is the coefficient from 2y to Tü- Note that, for any topological ordering of LiNGAM, if π (22) > I (24). then bYИ = 0. Thus, each Iz admits the following compact representation: 2; (=()))<(())) bijay Eg.
Identifiability Identifiability conditions for LiNGAMs [37] primarily concern the distribution of errors E;: under Gaussianity, distinct linear DGPs can admit the same joint distribution, making them impossible distinguish. Shimizu al. [37] generalize this intuition with independent component analysis (ICA) [4] to provide a multivariate identifiability condition for LiNGAMs (see Appendix A.1). In this section, we adopt the aforementioned condition.
We first establish the connection between ancestral relationships and active causal paths. We exhaustively enumerate and define the potential pairwise causal ancestral path relations in Figure 2 and Lemma 3.1 (proof in Appendix A.2):
Ancestral Relations and Active Causal Paths
Lemma 3.1 (Active Causal Ancestral Path Relation Enumeration). Each pair of distinct nodes R. E V can be in one of four possible active causal ancestral path relations: API) no active path exists between T4, 2jj AP2) there exists an active backdoor path between 24, F.j, but there is no active frontdoor path between them; AP3) there exists an active frontdoor path between Iı, F,, but there is no active backdoor path between them; AP4) there exists an active backdoor path between Fe, Tj. and there exists active frontdoor path between them.
Next, in Lemma 3.2, we summarize the connection between causal paths and ancestral relationships (proof in Appendix A.3):
Lemma 3.2. The ancestral relationship between pair of distinct nodes I I4,23 E V can be expressed using active causal path relations: Ii are not ancestrally related if and only if they are in API or AP2 relation; and Ti, 3 are ancestrally related if and only if they are in AP3 or in AP4 relation.
The active causal ancestral path relation of a pair of nodes Ii, that are not ancestrally related can be determined through marginal independence testing and sequential univariate regressions as illustrated in Lemmas 3.3 and 3,4 (proofs in Appendices A.4, A.5):
Lemma 3.3 (AP1). Vertices z4, Iy are API relation if and only if xz TL 1..-
Lemma 3.4 (AP2). Let M be the set of mutual ancestors between a pair of vertices and æ.j, iee. M = An(xs) O An(zz)) Let x.M. 2 be the result of sequentially regressing all mutual ancestors in M out of æis Jy with univariate regressions, in any order. Then, let r! be the residual of x. regressed on x,M and r be the residual of a regressed on x.5 Suppose 1 Zj. Then, I4, It are in AP2 relation if and only ip 1L a.M and rš 11. xr/1.
If pair of nodes I2 T, is ancestrally related, fully ascertaining their ancestral relation involves discerning between the ancestor and descendent. As illustrated Lemmas 3.5 and 3.6 (proofs in Appendices A.6, A.7), this can be determined through marginal independence testing after sequential univariate regressions with respect to the mutual ancestor set.
Lemma 3.5 (AP3). Let r be the residual of the V, regressed on Ti. and r be the residual of xn regressed on Zj. Vertices Is, 2j are in AP3 relation if and only if 1L Ij and one of the following

---

Algorithm LHTSS: Linear Hierarchical Topological Sort
input: features .1,.....
.,, Fd E
2: output: hierarchical topological sort (V).
16: Stage 3: AP2 and AP4 Relations
3: initialize: ancestral relations set ARS.
17:
while 3. with unknown relations do for 24, with unknown relations do Sequentially regress Vi, 2y on mutual ancestors, store final residuals rf. if r IL I, and rr TT T3 then Store Ti are not related in ARS.
4: Stage I: API Relations
18:
5: for all pairs Vi, J. do
19:
6: if xt; 11. 2. then
7:
Store Tit J. are not related in ARS.
20:
8: Stage 2: AP3 Relations
21:
9:
for all Tı; Ry with unknown relations do
10:
Set r as residual of æj regressed on Set ri as residual of z; regressed on 2j.
22:
else if r/ IL 2; and r3 A xy then Store 2, An(a,) ) in ARS.
11:
23:
if r IL æ4 and r A zÿ then
24:
else if rf 4 2; and rj IL Iy then Store 2y E (An(zs) in ARS. 26: Stage 4: Obtain sort by subroutine AS
Store xy E An(zy) ( in ARS.
25:
14:
rt AL x4 and r IL Zi then
15:
Store n. E An(x;) in ARS.
27: return: (L(V) + AS(ARS)
holds: 1) 23 IL ri and 2j 1 riy, corresponding to Ha E An(zz), or 2) 1 rf and zj TL rj. corresponding to Iz E An(a:))
Lemma 3.6 (AP4). Let M, Tl ",,3," M .IA be as defined in Lemma 3.4. Suppose Is, Tı,T; not in AP3 relation. Then, De, Ty are AP4 relation if and only if a4 V Fj and one of the following holds: rf II. x:f and rj V aM corresponding to Is E An(x;), or 2) ri H xnd and r² 11 M corresponding to xy E An(xa).
Linear Hierarchical Topological Sort We propose LHTS in Algorithm 1 leverage the above results to discover a hierarchical topological ordering. Marginal independence tests and pairwise regressions are used identify active causal ancestral path relations between pairs of vertices.
LHTS discovers all pairwise active causal ancestral path relations: Stage discovers all API relations using marginal independence tests; Stage discovers all AP3 relations by detecting when pairwise regressions yield an independent residual in only one direction; Stage 3 iteratively discovers all AP2 and AP4 relations through marginal independence tests and pairwise sequential regressions involving mutual ancestors. Note that Stage can be viewed as special case Stage 3, where the mutual ancestor is empty. The process utilizes proofs provided in Appendices A.4, A.6, A.5, and A.7, respectivl,,..l.. Stage 4 uses the complete set of ancestral relations to build the final hierarchical topological sort via the subroutine AS (Algorithm in Appendix A.8) by recursively peeling off vertices with no unsorted ancestors. We show the correctness of Algorithm in Theorem 3.7 (proof in Appendix A.9) and subroutine AS (proof in Appendix A.8), as well the worst case time complexity in Theorem 3.8 (proof in Appendix A.10). We provide a walk-through of LHTS in Appendix A. 1..
Theorem 3.7. Given a graph G, Algorithm I asymptotically finds a correct hierarchical sort of G.
Theorem 3.8. Given n samples of d vertices generated by a LiNGAM, the worst case runtime complexity of Algorithm I is upper bounded by O(dn??..
LHTS can be seen as generalization of DirectLiNGAM [38], where they recursively identify root nodes vertices that have either API AP3 relations with every remaining vertex each step. The next lemma, Lemma 3.9 (proof in Appendix A.12) implies that LHTS discovers all root vertices in Stage 2. obtaining the first layer in the hierarchical sort. LHTS extends DirectLiNGAM by discovering AP2 and AP4 relations, allowing us to recover compact hierarchical sort.
Lemma 3.9. A vertex is a root vertex if and only if it has API or AP3 relations with all other vertices.
4 Nonlinear Setting
In this section, we develop a version of Algorithm for ANMs that feature only nonlinear causal functions f. Rather than detecting ancestor-descendant relationships, outline the connection
5

---

PP1
PP2
PP3
PP4
Z2
I2
C
JJ
J
Jjj
Jjj
Figure 3: Enumeration of the potential active causal paths among a fixed variable zj, one of its potential parents Tú, and C = PA(I,) \ IG. Solid arrows denote parenthood relations, and undirected dashed connections indicate the existence of active paths.
between active causal paths and parent-child relationships. We assume the unique identifiability of the nonlinear ANM as described by Peters et al. [25], and provide the conditions in Appendix B.I.
Nonlinear Topological Sort In the linear setting, we determined ancestral relations through sequence of pairwise regressions that led to independent residuals. However, a naive extension of this method into the nonlinear setting would fail, as regressions yield independent residuals under different conditions in the nonlinear case. For clarity, we demonstrate how LHTS fails to correctly recover causal relationships in exemplary 3-node DAG with nonlinear causal mechanisms.
Consider DAG G with three vertices æ2; 23, where I1 + 33 T2 + 33g. The functional causal relationships are nonlinear, given by 21 = E1, 22 E2, D3 = 22 + Es. where are mutually independent noise variables. We focus on whether LHTS can recover the parent-child relationship between I1 and F3. LHTS finds that the relationship between I1, is unknown in Stage 1. In Stage LHTS runs pairwise regressions between æ1, 4g but incorrectly concludes that I 43 are not in AP3 relation because neither pairwise regression provides an independent residual; both parents of Ty must be included the covariate set for an independent residual to be recovered.
To handle nonlinear causal relationships, we shift our focus to searching for different set of local substructures: the connection between active causal parental paths and the existence of parent-child relationships. We will use the existence of specific parent-child relationships to first obtain superset of root vertices, then prune away non-roots. Once all root vertices are identified, we build the hierarchical topological sort through nonparametric regression, layer by layer.
Given vertex zj and one of its potential parents 2;, first provide enumeration of all potential active casual parental path types between them with respect to C = PA(x;) \ de (which could potentially be the empty set) in Figure 3 and Lemma 4.1 (proof in Appendix B.2).
Lemma 4.1 (Active Causal Parental Path Relation Enumeration). Let x, J E V be a pair of distinct nodes, where I; is of the potential parents of z.J. Let C = PA(zz) \ 24. There are in total four possible types of active causal parental path relations between 2, and x, with respect to C: PP1) Is and zj are not directly causally related, and active path exists between D, and C; PP2) x; directly causes 3i (x: - Rj), and no active path exists between I; and C; PP3) xa directly causes Iz (at -> æzjh, and there exists an active path between Ia and C; PP4) V; and æ.j are not directly causally related, and there exists an active path between I; and C.
Next, for pair of distinct vertices F4, 2z E V, we establish the connection between pairwise inde- pendence properties and active causal parental path relations in Lemma 4.2 (proof in Appendix B.3). This allows us to reduce the cardinality of the potential pairs of vertices under consideration in the later stages of the algorithm.
Lemma 4.2 (Non-PP1). Vertices I,, E V are not in PPI relation if and only if x;; M. Tz-
In Lemma 4.3, we show that all pairs of vertices that are in PP2 relation can be identified through local nonparametric regressions (proof in B.4).
Lemma 4.3 (PP2). Let Ti, E P.3 = (Ir V ZER TL Ii, PL x333, r be the residual of xt nonparametrically regressed on I,, and ri,P be the residual of x nonparametrically regressed on I; and all x E Poj- Suppose x; and x; are not in PPI relation. Then, æs and x are in PP2 relation if and only if one of the following holds: 1) xg 1I. r! or 2) xz 11 r..p-
Condition 1) of Lemma 4.3 relevant when C Ø: in this case, pairwise regression identifies x, as a parent of Zy. Condition 2) is relevant when C' * 0: we leverage the independence of I; from the rest

---

Algorithm 2 NHTS: Nonlinear Hierarchical Topological Sort
,,...,,,, input: vertices ......, Fdd E V.
15:
PRS: : 24, are in PP2 relations, 3; E Pa(x)))
2: output: hierarchical sort (V).
3: initialize: parent relations set PRS.
16: Stage 3: Root Identification
4: Stage I: Not-PP1 Relations
17:
for 2; E PP2 relation do
5: for all pairs Viz I.j do
18: Check if z; does not d-separate any child from all marginally dependent vertices. If so, then update PRS xs is a root, sort x; into layer O, TT, = 0.
6: if' x;; IL z, then
7:
PRS: 2,, not in PPI relations.
8: if ar; is in PP1 relation with all vertices then
9: PRS: xy is isolated, sort ze. NT. (x;(); =
19: Stage 4: Obtain Sort
10: Stage 2: PP2 Relations
20: for l E (1....,d -ddd-1-- 1 do for all unsorted z; do
11: for all F4, By not in PPI relations do
21:
12:
Set ri as residual of 2 regressed on Ti. Set r p as residual of I, regressed on
22:
Set r, as the residual of z; regressed on sorted features in Tи. If ri _ll x, Vyy лH, add into the current layer (x;) k.
13:
23:
24, Pij.
14:
if I. IL r or Dj A rl. then
24:
return (V)
of zzys 's parents to generate Pij. set that contains C, but does not contain any of 2y descendants. an independent residual were to be recovered by nonparametrically regressing 2. onto z; and Pij. we identify 2; as parent of Tz.
Let W be the set of all parent vertices that are in PP2 relation with at least one vertex, i.e., the union of x; satisfying either condition of Lemma 4.3. We now show that all non-isolated root vertices contained in W, and they can be differentiated from non-roots in W that are also PP2 relations.
Lemma 4.4 (Roots). All non-isolated root vertices are contained in W. addition, x4 W is rool vertex if and only if 1) T4 is not a known descendant of any T.j E W, and 2) for each x. E W, either a) F2 11. F. b) Fy is PP2 relation to Fs, i.e., Fy E Ch(xx),, or e) there exists a child of xs. denoted by Ik. that cannot be d-separated from 3j given Fú, i.e., Ijj 1L 24/24
Lemma 4.4 relies on the following intuition: for any non-isolated root vertex o, and descendant Ih E De(x;). there exists child of r; that 1) lies on directed path between r; and zg. and 2) is in PP2 relation with z;. Vertex Ij will fail to d-separate this specific child from the marginally dependent non-root descendant Tg. On the other hand, non-roots in W must d-separate the children they are in PP2 relation to from all marginally dependent roots: this asymmetry allows non-roots to be pruned. See Appendix B.5 for detailed proof.
We propose a method, Algorithm 2, that leverages the above results discover a hierarchical topological ordering: we first use the above lemmas to identify the roots, then use nonparametric regression to discover the rest of the hierarchical sort.
In Stage 1, we discover all pairs 2;, Ty not in PP1 relation by testing for marginal dependence; in Stage 2, we leverage Lemma 4.3 to find the vertex pairs that are in PP2 relations, a superset of the root vertices; in Stage we prune non-roots by finding they d-separate their children from at least one dependent vertex in W (Lemma 4.4); in Stage we identify vertices in the closest unknown layer by regressing them on sorted nodes and finding independent residuals. We show the correctness of Algorithm 2 in Theorem 4.5 (proof Appendix B.6), and the worst case time complexity in Theorem 4.6 (proof in Appendix B.7). We provide walk- through of Algorithm 2 in Appendix B.8.
Theorem 4.5. Given a graph G, Algorithum 2 asymptotically finds a correct hierarchical sort of G.
Theorem 4.6. Given n samples of d vertices generated by identifiable nonlinear ANM, the worst case runtime complexity of Algorithm 2 is upper bounded by O(²nn))
The number of nonparametric regressions run by NHTS in each step is actually inversely related to the size of the covariate sets, while the number of regressions in each step of RESIT and NoGAM are directly proportionate to the covariate set size. We provide formal analysis of the reduction in complexity for the worst case (fully connected DAG) in Theorem 4.7 (proof in Appendix B.9):
Vertices Ha, æj are said to be d-separated by a set Z iff there is active path between æi relative Z.
7

---

M
Figure 4: DAG corresponding to Lemma 5.1, which tests whether 2; Pa(zy) (i.e., whether the red arrow exists).
Algorithm 3 ED: Edge Discovery
1: input: features 21,.... za E V, topologiaal
7:
Set x; equal to ith index of n.
order g.
8:
Set Zij Pa(x;) U Pa(æ;).
2: initialize: empty parent sets (Pa(z;)]
9:
Step 2: Test for Parenthood
3: Step 1: Loop over
10:
if x; 14 x, ZZ then Update 2; E Pa(zz). return (Pa(æ:))
4: for j E (2..... d) do
11:
5: Set equal to jth index of
12:
6: for i {j- .,..,, 1 do
Theorem 4.7. Consider a fully connected DAG = (V, F) with nonlinear ANM. Let d [V. Let "ik
NHTS be the number of nonparametric regressions with covariate set size k E [d 2] run by NHIS RESIT and RESIT NoGAM = k: 1. This implies that for all k >
when sorting V; we similarly define nik nik NoGAM respectively. Then, B2 NHTS d k, and
d NHTS < <nk
RESIT ni NoGAM
2. 1kk
5 Edge Discovery
Parent selection from a topological ordering v via regression traditionally a straightforward task in the infinite sample setting: for each vertex 24, establishes J; = {>x <(xx) < T(zx)), }. superset of x 's parents that contains none of x; 's descendants. A general strategy for pruning Pa(x,) from Jj is to regress H; on Jj and check which 2 E J; are relevant predictors. The key issuc is that J, grows large in high-dimensional graphs: current edge pruning methods either make strong parametric assumptions or suffer in sample complexity. Lasso and GAM methods impose linear and additive models, failing to correctly identify parents in highly nonlinear settings. RESIT assumes a more general nonlinear ANM, but requires huge sample sizes and oracle independence tests for accurate parent set identification: authors in [25] confirm this, saying "despite [our] theoretical guarantee[s] [s], RESIT does not scale well high number of nodes."
We propose an entirely nonparametric constraint- based method that uses a local conditioning set Zi to discover whether x, Pa(xy), rather than Jr outperforming previous methods by relaxing parametric assumptions and conditioning on fewer variables. The following lemma outlines a sufficient condition for determining whether an edge exists between two vertices (proof in Appendix C.I), D), visualized in Figure 4.
Lemma 5.1 (Parent Discovery). Let π be a topological ordering. xz. T, such that (xx:) < T(2)). Let Zi = Cij U M3j. where Cī = {zx T (#x) (Iz))). Then, 2; > x, > Ii A x, Zaj.
Pa(xx)), 2x A.. 2z }. Mij = {ar 2 E Pa(x;)) a <
The intuition is that to determine whether x; - xj. instead of conditioning on all potential ancestors of zj. suffices condition on potential confounders of x,, vy (Cty) and potential mediators between z; and x; (Mz)). This renders all backdoor and frontdoor paths inactive, except the frontdoor path corresponding to a potential direct edge from 2; to Zz-
Edge Discovery We propose an algorithm that leverages the above results to discover the true edges admitted by any topological ordering by running the described conditional independence test in a specific ordering. We give the implementation for pruning a topological sort here, but our approach can be generalized to a hierarchical version (see Appendix C.2). We show the correctness of Algorithm 3 in Theorem 5.2 (proof in Appendix C.3).
Theorem 5.2. Given correct topological ordering of G, Algorithm 3 asymptotically finds correct parent sets PA(z.), Vx, G.
The key insight to check each potential parent-offspring relation using the conditional independence test 2; 1L 2z) Zig such that previous steps in the algorithm obtain all vertices in both Cis and Myj.
8

---

We first fix a vertex 2j whose parent set we want to discover. We then check if vertices ordered before 2j are parents of z, reverse order, starting with the vertex immediately previous to 2; in This process starts at the beginning of meaning we discover parent-offspring relations from root to leaf (see Appendix C.3 for detailed walk-through and proof). We show the worst case time complexity of Algorithm 3 in Theorem 5.3 (proof in Appendix C.4).
Theorem 5.3. Given samples of d. vertices generated by model corresponding DAG G, the runtime complexity of ED upper bounded by O(dPn)).
6 Experiments
Setup Methods are evaluated on 20 DAGs in each trial. The DAGs randomly generated with the Erdos-Renyi model [5]; the probability of an edge is set such that the average number of edges in each d-dimensional DAG is d. Gaussian, Uniform, Laplace noise is used the exogenous In experiments for linear topological sorting methods (Figure 5), we use linear causal mechanisms to generate the data; in experiments for nonlinear topological sorting methods (Figure 6) and edge pruning algorithms (Figure 7), we use quadratic causal mechanisms to generate the data. Existing ANM methods are prone to exploiting artifacts that are more common in simulated ANMs than real-world data [28, 29], inflating their performance on synthetic DAGs and leaving real-world applicability an open question. To reduce concerns about such artifacts, data were generated with reduced R -orttabbiy [29], and standardized to zero mean and unit variance [28].
Metrics Atop is equal to the percentage of edges that can be recovered by the returned topological child is sorted before a parent). We note that Atop is normalized version of the topological ordering divergence Dtop defined in [30]. Edge pruning algorithms return list of predicted parent sets for each vertex: Fi 2 Precision Recall measures the performance of these predictions.
ordering (an edge cannot be recovered
Precision Recall
Linear Topological Sorts Figure 5 demonstrates the performance of our linear topological ordering algorithm, LHTS, in comparison with the benchmark algorithms, DirectLiNGAM [38] and R [29]. RR-Sort is a heuristic sorting algorithm that exploits artifacts common in simulated ANMs; both benchmarks are agnostic to the noise distribution. We observe that both DirectLiNGAM and LHTS significantly outperform -ort. LHTS demonstrates asymptotic correctness in Figure 5(c), achiev- ing near near-perfer sear-perfect Atop at n = 2000. However, LHTS has consistently lower Atop than DirectLiNGAMM in Figure 5(a). On the other hand, LHTS encodes more causal information: the orderings produced by LHTS in Figure 5(b) had roughly ~ 70% fewer layers than the orderings produced by DirectLiNGAM, reducing the size of potential parent sets J by identifying many non-causal relationships.
1.
08
06
Ontere
04
-4
UPS
02
-4-
Divet.B..GMM
1"2.Sort
04
Orderno SCOVTST
comsacmess
1.5
Grsertng
62
H"" 2-225
200
1200
300
300
50
1000
Ordenng accunacy
Ordernn contantinoo
Figure 5: Performance of LHTS on synthetic data. Top row: 500 with varying dimension d.
Bottom row: 10 with varying sample size n. See Appendix D.I for runtime results.
"Code:

---

Nonlinear Topological Sorts Figure 6 illustrates the performance of our nonlinear topological sorting algorithm. We take GES [3]. GRaSP [12], GSP [39], DirectLiNGAM [38], NoGAM [20] and Re-SSort [29] baseline comparators that are all agnostic to the noise distribution. We excluded and RESIT since general they perform much worse than baseline methods [20]. We note that as DirectLiNGAM, NoGAM, and NHTS are FCM-based methods, they each return unique topological ordering; however, as GES, GRaSP, and GSP are scoring-base -base methods [7]. they return only a MEC. All topological orderings contained within an MEC satisfy every conditional independence constraint in the data, and therefore are all equally valid. To enable a fair comparison, we randomly select ordering permitted by outputted MEC for evaluation. We note that NHTS outperformed all baselines, achieving the highest median Atop in all trials. Furthermore, as expected from Theorem 4.6, NHTS up to xx faster than NoGAM (see Appendix D.2). We provide additional experiments over DAGs with increasing edge density in Appendix D.3.
1.00
340
1.00
0.30
5.00
0.50
0225
0.00
KITS D.NGH KoGAN
Gludh
125T
Mins OUNGH NelaM
GES
Chudt
227
R2ST
WHTS
DUKCM NoGAH
GReSP
GSP
R2ST
Figure 6: Performance of NHTS on synthetic data, n 300, dimension d = 10, with varying error distributions: Gaussian, Laplace, Uniform (left, middle, right). See Appendix D.2 for runtime results.
Edge Pruning Figure illustrates the performance of our edge pruning algorithm ED. We take covariate hypothesis testing with GAMs (CAM-pruning [1]), Lasso regression, and RESIT as baseline comparators that are all agnostic to the noise distribution. All algorithms were given correct topological sorts: ED significantly outperformed all baselines, with the highest median F1 score in all trials. ED was slower than Lasso, but was significantly faster than the other nonlinear edge pruning algorithms, CAM-pruning and RESIT. RESIT was excluded from higher-dimensional tests due to runtime issues. The poor performance of baseline methods highlights the need for sample efficient nonparametric method for accurate causal discovery of nonlinear DGPs. We provide additional experiments in settings with increasing density and varying noise distributions in Appendix D.4.
04
03
7.5
2000
2200
Figure 7: Performance of ED on synthetic data, uniform noise. Left, middle: n 300 with varying dimension d. Right: d 10 with varying sample size n. See Appendix D.5 for runtime results.
Discussion In this paper we developed novel global causal discovery algorithms by searching for and leveraging local causal relationships. We improved on previous topological ordering methods by running fewer regressions, each with lower dimensionality, producing hierarchical topological sorts. Additionally, we improved on previous edge pruning procedures by introducing a nonparametric constraint-basee method that conditions far fewer variables achieve greater recovery of parent sets. We tested our methods on robustly generated synthetic data, and found that both our nonlinear sort NHTS and edge pruning algorithm ED significantly outperformed baselines. Future work includes extending the topological sorting algorithms to the full ANM setting with both linear and nonlinear functions. simultaneously exploiting both ancestor-descendent and parent-child relations, as well as adapting our approach to handle various forms of unmeasured confounding. Additionally, we aim to develop statistical guarantees of sample complexity for our methods, extending previous results [50] derived in the setting of nonlinear ANMs with Gaussian noise.
10

---

References
[1] Bühlmann, P., Peters, J., and Ernest, CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 2014. URL httt:://axxxii org/ abs/1310.. 1533.
[2] Cheng, D., J., Liu, L., Yu, K., Duy Le. and Liu, Toward Unique and Unbiased Causal Effect Estimation From Data With Hidden Variables. IEEE Transactions on Neural Networks and Learning Systems, 2023. URL https::/ieeeepploreee ieee. org/ /ttp///mmee////19996ee
[3] Chickering, D. M. Learning Equivalence Classes of Bayesian Network Structures. Journal of Machine Learning Research, 2002. URL httpps//arxiv..or org/ 2ats//332.3350
[4] Comon, P. Independent component analysis, a new concept? Signal Processing, 1994. URL
http://aa0SSceeeeeeaa/o000113/3 doccueent.
[5] Erdos, P. and Renyi, A. On the evolution of random graphs. Publicationes Mathemati- cae, 1960. URL /tttpp:://ubiiiiee math. unideb. hw/whhhulddae wwppppoii-10548611111211111
1959_6_3_4_12
[6] Gan, K., Jia, S., and Li, A. Greedy approximation algorithms for active sequential hypothesis testing. Advances in Neural Information Processing Systems, 2021. URL pesttss://dl. acm.
http///do/ppoceeinggs////// /10 35402654
[7] Glymour, C., Zhang, K., and Spirtes, P. Review of Causal Discovery Methods Based on Graphical Models. Frontiers in Genetics, 2019. URL https://www, http://wwwwoootrrriiin.n org/ article/10. 3389 fggee//eeee...e 2019. 00524/ful1.
[8] Gupta, S., Childers, D., and Lipton, 7. C. Local Causal Discovery for Estimating Causal Effects. In Proceedings the 2nd Conference Causal Learning and Reasoning (CLeaR), 2023. URL
http: //tp//aaxvvv////////v///ee//aaav org/ argg/b///aa ab////bss2222... 08070.
[9] Hoyer, P. O. and Hyttinen, A. Bayesian discovery of linear acyclic causal models. In Proceedings of the 25th Conference Uncertainty in Artificial Intelligence, 2009. URL https://arxit: org/ 120511255222 2641.
[10] Hoyer, P. O., Janzing, D., Mooij, J. M., Peters, J., and Schölkopf, B. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, 2008. URL (5555/29817800 2981866.
https: / /dl. acm. www.aa//......................a.
[11] Hoyer, O., Shimizu, S., Kerminen, A. J., and Palviainen, M. Estimation of causal effects using linear non-Gann Gan-Gassiinn causal models with hidden variables. International Journal of Ap- proximate Reasoning, 2008 http::
elsevier. com/retrieve/pii.
[12] Lam, W.-Y., Andrews, B., and Ramsey, J. Greedy Relaxations of the Sparsest Permutation Algorithm. Proceedings of the 38th Conference on Uncertainty in Artificial Intelligence, 2022.
URL. https::/
mlr. www.prees// Prressvv80/// v180/1am22a////// 2aamm222aamm2.....mm pdf.
[13] Lee, J. J., Srinivasan, R., Ong, C. S., Alejo, G. J., and Malinsky, D. Causal determinants of postoperative length of stay in cardiac surgery using causal graphical learning. The Journal of Thoracic and Cardiovascular Surgery, 2022. URL
D., Schena, S., Shpitser, I., Sussman, M., Whitman,
ott.epcc...cc..eerrrreeoo//////p elsevier. com www.///trieee/eeeeeeooceeeeeeeee pttt:////0//0022222000000
[14] Li, A., Lee, J., Montagna, F., Trevino, C.,, and Ness, R. Dodiscover: Causal discovery algorithms in Python., 2022. URL
[15] Linder, J., Koplik, E., Kundaje, A., and Seelig, G. Deciphering the impact of genetic variation human polyadenylation using aparent2. Genome biology, 2022. URL https: //pubmed. ncbi nlm. nih, gov /0v/33335397/.
[16] Maasch, J., Pan, W., Gupta, S., Kuleshov, V., Gan, K., and Wang, Local discovery by partitioning: Polynomial-time causal discovery around exposure-outcome pairs. In Proceedings of the 40th Conference on Uncertainy Artificial Intelligence, 2024. URL https: //axxxvx.. org/ aa//ab////// wtwp///.......... 17816.
11

---

[17] Marra, G. and Wood, S. Practical variable selection for generalized additive models. Com- putational Statistics and Data Analysis, 2023. URL https://wwx http://wwwwwcccccccrcettct URLntppss//owweeeeeeeddeeee science/ hammc//article/aaeeeeaa httpa////p.....10177711000000000
com/
[18] Marra, G. and Wood, S. Regression Shrinkage and Selection via the Lasso. Computational Statistics and Data Analysis. 2023. URL https://wwt sciencedirect, ctt/ sts/ssiiencceee article/abs (ericc//abb/pii//0166///33 pttpp//1////7900000000000
[19] Montagna, F., Mastakouri, A. A., Eulig, E., Noceti, N., Rosasco, L., Janzing, D., Aragam, B., and Locatello, F. Assumption violations in causal discovery and the robustness of score matching. In 37th Conference on Neural Information Processing Systems, 2023. URL http: //arxiv. orgx aaaabb/s1111 13387.
[20] Montagna, F., Noceti, N., Rosasco, L., Zhang, K., and Locatello, Causal Discovery with Score Matching on Additive Models with Arbitrary Noise. Proceedings of the 2nd Conference on Causal Learning and Reasoning, 2023. URL http:://azzii org/ httg//bb//30///////////// 03265.
[21] Montagna, F., Noceti, N., Rosasco, L.,, Zhang, K., and Locatello, F. Scalable Causal Discovery with Score Matching. Proceedings of the 2nd Conference on Causal Learning and Reasoning. 2023. URL http: //arxiv...iii///aa , org/ aba/ 2304/ 03382.
[22] Niu, W., Gao, Z., Song, L.,, and Li, L. Comprehensive Review and Empirical Evaluation of Causal Discovery Algorithms for Numerical Data, 2024. URL https://arxiv orgooogg/aa abs/ 2407. 13054.
[23] Pearl, J., Glymour, M., and Jewell, N. P. Causal inference in statisties: primer. Wi- ley, 2016. URL https://www dwwasciencasssnnnnnnnnennnnnnnne www.nnaa...org/.......e... /aitts//efffut//tee///////// 201IFFEEECCE2220002000TTTTTS
pdf.
[24] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. Scikit-leam: Machine learning in Python. Journal of Ma- chine Learning Research, 2011 URL https://jmlr csail mit. edu/paperes ceaperssvolumeet volumm12// wtwpeggssaaa/ppdreegggeeee 1aa.iaaa.... pdf
(25) Peters, J., Mooij, J., Janzing, D., and Schölkopf, B. Causal Discovery with Continuous Additive Noise Models. Journal of Machine Learning Research, 2014. URL http: htttt::/aaxxx..rrree org/ abs/ 1309.6773
[26]
C. and Torrecilla, J. dcor: Distance correlation and energy statisties in python. SoftwareX, 2023. URL https: //www. sciencedirect. wttp////c.//eecccccccccccccce.ec ecience/articce/ pi1//2352711023000222
[27] Raskutti, G. and Uhler, C. Learning directed acyelic graphs based on sparsest permutations. STAT, 2013. URL httts:://arti: wtw.aaiwww......aaa.aooooooo /ttp//11700000
[28] Reisach, A., Seiler, C., and Weichwald, S. Beware of the simulated dag! causal discovery benchmarks may be easy to game, Advances in Neural Information Processing Systems, 2021. URL https: h/ta:://axiiviae abs//b////// awtpa///..aa.e.ee/aaa/aeeaaaa 13647.
[29] Reisach, A. G., Tami, M., Seiler, C., Chambaz, A., and Weichwald, Scale-Invariant Sorting Criterion to Find a Causal Order in Additive Noise Models. 37th Conference on Neural Information Processing Systems, 2023. URL http:://ariiii http::///xxxv.oo//aa/saasse /3033 18211.
[30] Rolland, P., Cevher, V., Kleindessner, M., Russel, C., Scholkopf, B., Janzing, D., and Lo- catello, F. Score Matching Enables Causal Discovery of Nonlinear Additive Noise Mod- els. Proceedings of the 39th International Conference on Machine Learning, 2022. URL https: //tp/aaxvvv////////////////aaava wwwxi..rg////eeee.. 2203. 220..2441...
[31] Runge, J., Bathiany, S., Bollt, E., Camps- G., Coumou, D., Deyle, E., Glymour, C., Kretschmer, M., Mahecha, M. D., Muñoz-Marí, J., van Nes, E. H., Peters, J., Quax, R., Re- ichstein, M., Scheffer, M., Schölkopf, B., Spirtes, P., Sugihara, G., Sun, J., Zhang, K., and Zscheischler, J. Inferring causation from time series in Earth system sciences. Nature Commu- nications, 2019. URL http:/ httt:://wwww .wwura. com/ http///icceee//1477-011/1111
12

---

[32] Sanchez, P., Liu, X., O'Neil, A. Q., and Tsaftaris, S. A. Diffusion models for causal discovery via topological ordering. Proceedings of the 39th International Conference on Machine Learning. 2023. URL https: h //txxx/xvvvv///////////aa/./.av org/ h/////b///// /tts/220...o 2220022211.02.....22122222111111
[33] Seabold, S. and Perktold, J. statsmodels: Econometric and statistical modeling with python. In 9th Python in Science Conference, 2010. URL //wp::///wwwwww https://www.stattmmmdee........ wrggstablas index. html.
[34] Shah, A., Shanmugam, K., and Ahuja, K. Finding Valid Adjustments under Non-ignorability with Minimal DAG Knowledge. Proceedings of the 25th International Conference on Artificial Intelligence and Statistics, 2022. URL http: arxxiv. org/abs/ /0bb//216.... 11560.
[35] Shah, A., Shanmugam, K., and Kocaoglu, M. Front-door adjustment beyond markov equivalence with limited graph knowledge. Advances in Neural Information Processing Systems, 2023. URL
https://aree org/ abs/2366. 11008.
[36] Shimizu, Lingam: Non-Gaussian Methods for Estimating Causal Structures. Behaviormetrika.
2014. URL http: //link. springer. com/10 //3k///kkk.... 41.44 65.
[37] Shimizu, S., Hoyer, P. O., Hyvarinen, A., and Kerminen, A. A Linear
Acyclic Model for Causal Discovery. Journal of Machine Learning Research, 2006. URL https:
//dl. /cml.acm.e.. www//oi//////ooooooooo / 10. 5555/ (1655/1248544 1248547 1248619.
[38] Shimizu, S., Inazumi, T., Sogawa, Y., Hyvarinen, A., Kawahara, Y., Washio, T., Hoyer, P. O., Bollen, K., and Hoyer, Directlingam: direct method for learning a linear non-gaussian structural equation model. Journal of Machine Learning Research, 2011. URL h/tps://dl.
acm. bogg/oo//ii 5555/15553118 1953049. 2021040.
[39] Solus, L., Wang, Y., and Uhler, C. Consistency guarantees for greedy permutation-bawee causal inference algorithms. Biometrika, 2021. URL httpttt//ccadeeeee oup. com/biometi article- secccleebbeeee
[40] Spirtes, P. An Anytime Algorithm for Causal Inference. Proceedings of Machine Learning Research, 2001. URL https::/prooceedings mlr. httpaa/////ppprtestaappppppppeee pdf.
[41] Spirtes, P. and Zhang, K. Causal discovery and inference: concepts and recent method- ological advances. Applied Informatics, 2016. URL https http//apppppdiinoeeeeeeeeeeeeeep sondd-ifoommaiccs-1 springeropen. cwt.eeococarrraa cow//article///i/aocccoocooaeeee
[42] Spirtes, P., Glymour, C., and Scheines, R. Causation, Prediction, and Search. Springer New York, 2000. http:/tp://ink springer. 1001101100////1114412222277
[43] Squires, C. Graphical Model Learning, 2021. URL graphical_model__ poodeellearnig.
[44] Steeg, G. Non Naraperamt Entropy Estimation Toolbox, 2014. URL https://///ttt com/ Progversteog/MPEE:
[45] Takashi, I,, Mayumi, I., Yan, Z., Takashi Nicholas, M., and Shohei, S. Python package for causal discovery based on lingam. Journal of Machine Learning Research, 2023. URL https://www.
org/ wwtggpapers//eeeee www.///4/211100022111111
pdf..
[46] Wang. S. K., Nair, S.. Li, R., Kraft, K. Pampari, A., Patel, A., Kang, J. B., Luong, C., Kundaje, A., and Chang. H. Y. Single-cell multiome of the human retina and deep learning nominate causal variants in complex eye diseases. Cell genomics, 2022. URL https://wwet sciencedirect.coom
htwm:/ccenneeeen wwt.///cinnee//rtccccceeeeeeeeee wtthee/rrtillleeee pttpaaa.....eee/////66672000000.
[47] Wong, A. K., Sealfon, R. S., Theesfeld, C. L., and Troyanskaya, O. G. Decoding disease: from genomes to networks to phenotypes. Nature Reviews Genetics, 2021. URL https: //www. nature. http///rtcceee//1156--222222
[48] Zhang, K. and Hyvarinen, A. On the Identifiability of the Post-Nonlinear Causal Model. Uncertainty in Artificial Intelligence, 2009, URL https: / /arxivv. weerivi.og///e. org/pdf/ 120// 2599.
13

---

[49] Zheng, Y., Huang, B., Chen, W., Ramsey, J., Gong, M., Cai, R., Shimizu, S.. Spirtes, P., and Zhang, K. Causal-leam: Causal discovery in python. Journal of Machine Learning Research, 2024. URL www.h:tp://te //axxvv. (errg/bbs/33000 16405.
[50] Zhu, Z., Locatello, F. and Cevher, V. Sample Complexity Bounds for Score-Matching:: Causal Discovery and Generative Modeling, 2023. URL
attps::/prroceedingg..t neurips wwtppppprriieeeeeeeeeeee wtt._illle//ppper/2////eeee 0BB8:55229ccao09944633344dd Coneadd--pppeeeeeeeeeeeeee pdf.
14

---

Appendix
A Assumptions, Proofs, and Walkthrough for LHTS
A.1 Identifiability in Linear Setting
As an example [25]: suppose X, N are normally distributed, X TT N, and ax + N. If ä x'Var(X)+a
xV)r(X) * à and N = X aY, then N TT and the following DGP also fits the same
distribution of X, Y. N: X = aY N. To generalize the intuition that non-Gaussia could lead to identifiability, [37] use independent component analysis [4] to provide the following theorem:
Theorem A.1. Assume a linear SEM with graph G, where I; = EugyeMee bizTj + Eis Vj 1.....d d, where all E are jointly independent and non-Gaussian distributed. Additionally, for each j E (l.....d} we require bi for all Pa(x:). Then, G is identifiable from the distribution.
We assume that the identifiability conditions described in Theorem A. hold throughout Section 3.
A.2 Proof of Lemma 3.1
Lemma 3.1 (Active Causal Ancestral Path Relation Enumeration). Each pair of distinct nodes Dj, E V can be in one of four possible active causal ancestral path relations: API) )) no active path exists between T4, AP2) there exists active backdoor path between F4, Z. but there is active frontdoor path between them; AP3) there exists an active frontdoor path between Is, Ej but there is no active backdoor path between them; AP4) there exists active backdoor path between Ze, Zz. and there exists an active frontdoor path between them.
Proof of Lemma 3.1. For pair of nodes x; and Ij. either x; JL Ty or x; 1/. Tj. If the former, then T; Tj are in AP1 relation. Suppose the latter is true: di ML. Vj implies 3 at least one active path between ni; Myj. Active paths either be backdoor frontdoor paths: therefore either frontdoor path exists, a backdoor path exists, or both a frontdoor and backdoor path exist. Thus, 24, Iy are either in AP2, AP3, or AP4 relation.
A.3 Proof of Lemma 3.2
Lemma 3.2. The ancestral relationship between pair of distinct nodes xı, Ii,35 E V can be expressed using active causal path relations: TG are not ancestrally related if and only if they are in API or AP2 relation; and Ty are ancestrally related if and only if they are in AP3 or in AP4 relation.
Proof of Lemma 3.2. Note, Zi, Ty are ancestrally related if and only if there exists an active frontdoor path between them. The conclusion follows immediately.
A.4 Proof of Lemma 3.3
Lemma 3.3 (AP1). Vertices Ie, Fy are API relation if and only if x; Tz--
Proof of Lemma 3.3. Note, JL xy if and only if there does not exist active causal path between them. The conclusion follows immediately.
A.5 Proof of Lemma 3.4
Lemma 3.4 (AP2), Let M be the set of mutual ancestors between a pair of vertices and æjj i.e. M = An(xx,) An(zz) Let x:M. 2 be the result of sequentially regressing all mutual ancestors in M out of 24, H. with univariate regressions, in any order. Then, let r be the residual of aM regressed M. on x:M and r² be the residual of x.M regressed on x Suppose H, V x... Then, I4,T, are in AP2 relation if and only if 1L xM and r 11 NM
15

---

Proof of Lemma 3.4. Under an identifiable LiNGAM, and given a linear topological ordering k V H; and øj admit the following representation:
Σ <k())
24
Fm + Eis
Σ
2j
ajmEm + Ey:
k(m)<<(f)
where Em are jointly independent noise terms and non-Gaussian distributed, following from Theo- rem A.I. Additionally, for each m k(m) < k(i) we require dim. + 0 for all m E Pa(x;). Similarly for jm
We first show the forward direction. Suppose Ti, 2y are in AP2 relation. Note, z; M zz. Let M be complement M, i.e. M V M. Then, the following decomposition Ta, 2zz:
T2 Σ
Σ
EmimE +
OsmEm + Ess
Jmm E AM
Ame MrAn(za)
Σ
Tj
jm +
Cjjm + Es.
Pm EM
Then, the x,M, _ have the general form
n.
Σ
Em Exi
Tne FY
Σ
BimEm + +
j
Fm EZ
where Y Mn)n()))) Z S M n An(a))) Note, Y Z U. As e, are all mutually independent, this implies that rf 11 rj.
We now show the reverse direction. Suppose x; 1 Zj. 7. 1 xz and r TL M " Note that I; IL.. Ty Ti, 2y are not in API relation. Suppose for contradiction that x,, Zj are in AP3 AP4 relation, where x; E An(z;) (WLOG). Then 3 frontdoor path between 2;, Ty. WLOG x; E An(r,). Regressing the descendent
the ancestor will always result in dependent residual, even after projecting out the influence of the mutual ancestors. Therefore, 2. If. rj. leading to a contradiction. Therefore, Ty must be in AP2 relation.
A.6 Proof of Lemma 3.5
Lemma 3,5 (AP3). Let ri be the residual the 2, regressed on H, and r² be the residual of xy regressed on Zz.. Vertices Ii, 2j are in AP3 relation if and only if; 1/. Ij and one of the following holds: Fj IL r and æj 14 ri corresponding to I; E An((z)), or 2) x4 14 ri and xz 1L r corresponding to xz An(x:;).
Proof of Lemma 3.5. We first show the forwards direction. Suppose 24, are in AP3 relation. Note an active path exists between Fe, Ij. so x; P1 zy. WLOG, let x; An(zy). Then, 2i, Zy admit the following representation:
22 = Z;
2, bystt +
ajm + Ej.
Tic A.n(
Note, there does exist a backdoor path between x; æy, we have F; IL I' Vxm An(x;) xz Therefore, r IL z;- Note that H, E De(r;): pairwise regression will leave rf H xy as this an example of reverse causality, special case of endogeneity [23].
16

---

We now show the reverse direction. WLOG, suppose 24, 2 satisfy condition 1) in Lemma 3.5. As I; fi y,, they cannot be in API relation. Suppose for contradiction that they are in AP2 or AP4 relation. Then, F backdoor path between 4;;; 4.j and confounding 22 on that backdoor path. As zz confounds 3i and is not adjusted for in pairwise regression, have n; V/. r/ 2j x). contradiction. Therefore, Ti, 2j must be in AP2 relation.
A.7 Proof of Lemma 3.6
Lemma 3.6 (AP4), Let M. ry F3. 2i , be as defined in Lemma 3.4. Suppose 24, 24,2; not in AP3 relation. Then, De, F. are AP4 relation if and only if 24 V Tj and one of the following holds: ri IL z and rj 4 23 corresponding to H, E An(xz), or 2) ri + xa and ry 11 23 corresponding to xy E An(aa)..
xM
Proof of Lemma 3.6. We first show the forward direction. Suppose n,, 2y are in AP4 relation. Note that active path exists between 24, Zy, 2; 1i zz. WLOG, suppose zı. An(zz). Note, x,M xx) are the result of projecting mutual ancestors M out of n,, By. Therefore, z; aay"" admit the following representation:
M 2%
Σ
Oea em
Um EMrAn(an)
2j
M
ajiE +
ajjmE
4m MTAn((z)))
Note that M n Ann(n; (MAn(( = Therefore, rf IL x.A. Note that 6; E both x. ". therefore rj HL z.M.
M
We now show the reverse direction. Suppose condition 1) holds, i.e., T, M Iy and WLOG of a.M .r H4. =n- Suppose for contradiction that T4 not in AP4 relation. They cannot be in AP1 relation because x; V. xy.- They cannot be AP2 relation because x. M. rj. By assumption they
not in AP3 relation. Therefore, we have a contradiction, therefore they are in AP4 relation. C
A.8
Ancestor Sort
Algorithm 4 AS: Ancestor Sort
1: input: set of ancestral relations ARS.
2: initialize: hierarchical topological order O []. remaining variables RV [.1,...; Fa], layer []
3: while RV is not empty do
4: Stage Determine Roots
5: for vi RV do
6: if xx has no ancestors E RV then Append 24 to T..
7:
8: Stage Update Sort and Remaining Variables
9:
Append L to O. Remove all vertices in from RV.
11:
return O
Overview In each iteration, this algorithm uses the set of ancestral relations to identify which vertices have no ancestors amongst the vertices that have yet to sorted. peels those nodes off, adding them as a layer to the hierarchical topological sort. First the roots are peeled off, then the next layer, and so on.
Proof of Correctness
Proof. The input to the algorithm is an ancestor table ARS cataloging all ancestral relations for every pair of nodes in an unknown DAG G. Let H be the minimum number of layers (H 1 the length of the longest directed path in the graph) necessary in valid hierarchical topological sort of G.
17

---

Base Cases (1,2)
1. H = 1. None of the nodes have ancestors, all nodes are added to layer 1 in Stage 1.
2. H 2. Nodes with ancestors are added to layer in Stage I, and then removed from the remaining variables in Stage 2. As H = 2, the longest directed path is 1, so all nodes in the remaining variables must have ancestors in the remaining variables (otherwise, this would make the longest directed path 2). The nodes are added to layer 2 in Stage 1, and removed from the remaining variables in Stage 2. Remaining variables is empty, so now the order is correctly returned.
Inductive Assumption For any graph G with H < k, Hierarchical Topological Sort will return minimal hierarchical topological ordering.
We now show that hierarchical topological sort now yields minimal hierarchical topological ordering for G where
1. In the first iteration of Stage 1, nodes with no ancestors will be added to layer 1, then removed from the remaining variables in Stage
2. At this point, note that we can consider the induced subgraph formed by the nodes left in remaining, G'. The minimal hierarchical topological ordering that represents G' must have k k- 1 layers. It cannot be greater than k 1, because G' is subgraph of G and we removed nodes in layer of, reducing the maximal path length by 1. It cannot be less than - 1, because that would imply that not the minimal number of layers needed represent G.
3. By Inductive Assumption, Hierarchical Topological sort will return the a correct minimal hierarchical topological ordering for the induced subgraph G'.
4. Appending the layer and the sort produced by G' yields the full hierarchical topological sort for G.
Inductive assumption is satisfied for H = k, so for graph G with arbitrary H, the algorithm recovers the correct hierarchical topological sort.
A.9 Proof of Correctness for Linear Hierarchical Topological Sort
Theorem A.2. 3.7 Given graph G. Algorrthhm asymptotically finds correct hierarchical sort of G.
Proof. The goal is show that all ancestral relations between distinct pair of nodes Wi, By V determined by LHTS in Stages 1, 2, and 3 are correct. By A.8, Stage 4 will return a valid hierarchical topological sort given fully identified ancestral relations.
Stage identifies 24, Ty as in API relaton if and only xz IL Øy: follows by Lemma 3.3 that vertices in AP1 relation are correctly identified, and vertices in AP2, AP3 and AP4 relation are not identified.
Stage 2 identifies Ty unidentified in Stage as in AP3 relation if and only if after pairwise regression, one of the residuals is dependent, and the other independent: it follows by Lemma 3.5 vertices in AP3 relation are correctly identified, and vertices in AP2 or AP4 relation are not identified.
Consider hierarchical topological sort of DAG G TH, with h layers. Note, Layer 1 of मh consists of root nodes. Ancestral relations between root nodes and all other nodes are of type AP1 and AP3, and were discovered in Stage and Stage 2.
We induct on layers to show that Stage recovers all ancestral relations of type AP2 and AP4, one layer at time in each iteration.
18

---

Base Iterations (1,2)
1. All ancestral relationships are known for nodes in layer of тHe therefore all mutual ancestors of layer 2 and every lower layer are known. Then, by Lemma 3.4 and Lemma 3.6 ancestral relations of type AP2 and AP3 between vertices in layer 2 and all lower layers discovered in iteration 1.
2. All ancestral relationships are known for nodes in layer and 2 of THe therefore all mutual ancestors of layer's 3 and every lower layer are known. Then, by Lemma 3.4 and Lemma 3.6 ancestral relations of type AP2 and AP3 between vertices in layer 2 and all lower layers are discovered in iteration 2.
Iteration it 1, Inductive Assumption We have recovered all ancestral relationships nodes in layers 0
1. As ancestral relationships of nodes in layers k are known, mutual ancestors of vertices in layer relations of type AP2 and AP4 between vertices in layer and every lower layer are discovered.
k and every lower layer are known, so by Lemma 3.4 and 3.4 all ancestral
Iteration k: Inductive Assumption is satisfied for iteration k. therefore recover all ancestral relations of type AP2 and AP4 and 4) for nodes in layers O to k. Thus, for a DAG with an arbitrary number of layers, Stage 3 recovers all ancestral relations of type AP2 and AP4.
A.10 Time Complexity Proof for Linear Hierarchical Topological Sort
Proof. In Stage 1, LHTS runs O(d)) marginal independence tests that each have O(n²) complexity, In each step for Stage and Stage 3, L.TSS runs (dd)) marginal independence tests each with (nn)) complexity. In the worse case of fully connected DAG, there are d² steps in total, across Stage and Stage 3: this is because in each step one layer of the layered sort DAG is identified, and a fully connected DAG has d layers. Therefore, the overall sample complexity of LHTS is Odd
A.11 Walk-through of LHTS
The following diagram illustrates each stage of LHTS on an exemplary 5-node DAG;
Stage 1
Stage 2
Stage 3
Stage
Stage
A
A
C
C
B
C
D
D
D
D
E
E
E
All vertices are dependent on each other, so no AP1 relations are discovered in Stage 1. In Stage 2, vertex A is discovered in AP2 relation to all vertices, and vertex D is discovered in AP2 relation to vertex E. In Stage 3, vertices B and are discovered to be in AP3 relation to each other; then, we discover vertex B in AP4 relation to vertex D, and vertex C in AP4 relation vertex D. Therefore, in Stage 4, we recover the topological sort after running the AS subroutine.
19

---

A.12 Proof of Lemma 3.9: Root Active Causal Ancestral Path Relations
Lemma 3.9, vertex is root vertex if and only if it has API or AP3 relations with all other vertices.
Proof of Lemma 3.9. We first show the forward direction. If z; is root, then De(r;) = 0. Therefore, cannot have backdoor path between and any other vertex, therefore must be either AP1 or AP3 relation with other vertices.
We now show the reverse direction. If xj in API AP3 (as an ancestor) relations with all other nodes, cannot have any parents. Therefore, dq is root.
20

---

Pe.I...) satisfies Condition B.I, and fj are B Assumptions, Proofs, and Walkthrough for NHTS
B.1 Identifiability in Nonlinear Setting
Following the style of [19], we first observe that the following condition guarantees that the observed distribution of pair of variables 2y can only be generated by unique ANM:
Condition B.1 (Hoyer& & Hyttinen (9]). Given bivariate model I; Eis J.j = fz(24) + generated according to (1), we call the SEM an identifiable bivariate ANM the triple Pe, does not
solve the differential equation K" = k" 8 8²' + +D) -2g" f" f + g pr + s' 9 g' all wi, 2ij such that ('(aaag" [](]]) + where Pe,: are the density of E,, -- jj,k= log Pe.. Pe,, The arguments 2j - (z:), 24 and x; of g, k and f respectively, are removed for readability.
for
8
There is generalization of this condition to the multivariate nonlinear ANM proved by [25]: Theorem B.2. (Peters et al. [25]). An ANM corresponding to DAG G is identifiable if E V, E Pa(a;) and all sets S S V with Pa(x;) {i} C S De(i) \ {xi, 3 Xs with positive joint density such that the triple (1,(Pa()) non-constant in all arguments.
We assume that the identifiability conditions described in Theorem B.2 hold throughout Section 4.
B.2 Proof of Lemma 4.1
Lemma 4.1 (Active Causal Parental Path Relation Enumeration). Let Iq J V be a pair of distinct nodes, where I; is of the potential parents of z.j. Let C = \ 24. There are in total four possible types of active causal parental path relations between z, and x, with respect to C: PPI) I; and zj not directly causally related, and active path exists between n, and C; PP2) x directly causes zj (æ; + Z5), and no active path exists between I; and C; PP3) æs directly causes nj (x; -> æz2h. and there exists an active path between a; and C; PP4) Y; and zj are not directly causally related, and there exists an active path between I; and C.
Proof. Either T; II. My or 24 Zyz if the former, then 22, in PP1 relation. Suppose the latter true. Then, either x; C or 24 V. C: if the former, then Tj are in PP2 relation. Suppose the latter is true. Then, either I, PA(x;) or æ4 & PA(x)). If the former, then 24, Ty are in PP3 relation, and if the latter, then 24, Ij are in PP4 relation.
B.3 Proof of Lemma 4.2
Lemma 4.2 (Non-PPI). Vertices 2., V are not in PPI relation if and only if x; 11 xz.
Proof. We first show the forward direction. If v,, Ty are not in PP1 relation, then either o; E Pa(x)) or there exists active causal path between 24, 2y therefore, x; 14.4 zz--
We now show the reverse direction. If vx JK æy there exists an active causal path between By Ty therefore they cannot be in PP1 relation.
B.4 Proof of Lemma 4.3
Lemma 4.3 (PP2). Let 2i, J.j E
CLAKEV: nonparametrically regressed on Ii, and Th.p p be the residual of x.j nonparametrically regressed on x;; and all x E Pej. Suppose x, and x, are not in PPI relation. Then, æs and x are in PP2 relation if
T,, Tk f1 x.3,, r be the residual of az
and only if one of the following holds: 1) 1L r or 2) H₁ 1L r..pp
Proof. Note, there exists two sub cases of PP2 relation: PP2a) z; is the only parent of æ.,, i.e. |C| =0 0 and PP2b) Vi not the only parent of x,, i.e., [C| 0.
21

---

We first show the forward direction. Suppose 2;, are in PP2a relation. Then, as the
= ر (ma) 6
Tj TT r.. Suppose that Ti, are in a PP2b relation. Consider P.y: note that æ; TT C, and æj 4 C. C' Pij. Note that De(r;) J 24, therefore De(r;) Piz. Therefore, Pij contains all parents of zj. and excludes all descendent 2y. Then, as
2j fz C) t
we have 2i IL rjp
We now show the reverse direction. Suppose z; T1 r/. Note by assumption Zq. 2. are not in PP1 relation. Suppose for contradiction that 24, 2y are in PP3 or PP4 relation. Then, there must exist at least one active path between 25 that goes through C. If 3 backdoor path from 2,, 2j, then at least one vertex in C' is confounder of 24; Zy, and thus x; 14. r/: therefore, there must exist frontdoor path between 24, Iy through C, with at least one mediator 2g.- Note that Ty is dependent on both x; and zj,, and therefore is excluded from Pej. The exclusion of Iy introduces omitted variable bias [23], and thus F, V. r leading to contradiction. Thus, Fy; 2y are in PP2 relation.
B.5 Proof of Lemma 4.4
Lemma 4.4 (Roots). All non-isolated root vertices contained in W. In addition. z; W is rool vertex if and only if 1) xa is not a known descendant of any 3. E W, and 2) for each Ij E W, either a) F2 II. F.j, b) x in PP2 relation to Vú, i.e., Fj Ch(ar)), or c) there exists child of xa. denoted by Tk, that cannot be d-separated from Tj given Tú ie., Tz V. F4/34
Proof. We first show that all non-isolated root vertices are contained in W. Note, if a root x; is not isolated, it has at least one child zy. By definition, 2; has no parents, so there cannot exist a backdoor path between 24, Ty. Consider 2, E Ch(x;) that is in a hierarchical layer closest to zy. Note that there cannot exist an active path from 2; to Pa(zz) \ Ig, because that would imply that v, is not the child of 2, that is in the closest hierarchical layer. Therefore, 2 must be in PP2 relation, and a; E W.
We now show the forward direction. Suppose x; is root vertex. Then, I; has parents, so it cannot be the descendent of any vertex, let alone any vertex in W: Condition 1) is satisfied. For each Iy W, there either exists an active path between IY, T2j there does not. If there does not exist an active path, then I; IL Tj. satisfying Condition a). If there does exist an active path, then as I; is a root, x, De(r;)): the active path must be a frontdoor path. Consider the vertex Ty E Ch(rr)) that is in hierarchical layer closest to zy. AND is along the active frontdoor path from 2, to Iy.. Note that F; is in PP2 relation with Zg. If æ Iy. then Condition b) is satisfied. If n3 * Zj. then we note that there is an active frontdoor path between and z. Therefore, Tk, Rj are dependent even after conditioning on the mutual ancestor 24, i.e., Iz V. zh ry Condition e) satisfied. Therefore, Condition 2) is always satisfied.
We now show the reverse direction. Suppose Condition 1) and Condition 2) are satisfied for xy E W. Suppose for contradiction that xj is not a root. Note that Dr, E W such that x; is a root and Iy E De(x;). Consider any vertex x E Ch(zr) such that æj is in PP2 relation to Tg. Note that all active paths between F; and Ty must be frontdoor paths, and xy must lie along all these active frontdoor paths: otherwise, this would contradict that I, and Iy are in PP2 relation. Therefore. I; IL Z&) zz. This contradicts Condition 2c), therefore Condition 1) and Condition 2) cannot hold true for x4 E W that not a root.
B.6 Proof of Correctness for Nonlinear Hierarchical Topological Sort
Theorem 4.5. Given graph G, Algorithm 2 asymptotically finds a correct hierarchical sort of G.
22

---

Proof. This proof can be broken into two parts: we first show 1) all root vertices are identified after Stage 3, then we show that Stage 4 recovers a correct hierarchical topological sort given the root vertices.
Note that Stage identifies 2e, my not in PP1 relation if and only if x; xyy it follows by Lemma 4.2 that these identifications are correct.
Note that Stage 2 identifies the set æ; E W if and only if Bay such that either 1) 24 11 r! 2) Fy TT ri.p: follows by Lemma 4.3 that these identifications are correct.
Note that Stage explicitly uses condition from Lemma 4.4 to find all non-isolated root vertices it therefore correct. Note that all isolated root vertices were identified in Stage 1. Therefore, roots are identified.
Consider hierarchical sort of DAG T.A suppose the maximal directed path in G has size h. Note that layer O of TT consists of root nodes, and therefore has been identified. We note the following observation: when F; nonparametrically regressed on all Ihk E TL to obtain residual rt. if an ancestor of a; not in TL, then r; will be dependent on at least one ar due to omitted variable bias [23].
We now induct on the layers of Ty, to show that Stage 4 correctly recovers all layers.
Base Iterations (1,2)
1. All roots are identified, so layer of æy. is identified.
2. By observation, only z; in layer 1 of TT will have independent residuals r; after nonpara- metric regression on E TT, so layer 1 is correctly identified.
Iteration k Inductive Assumption We have recovered all layers of a up to layer k 1.
By observation, only x, in layer k of n will have independent residuals r; after nonpara- metric regression 2; E TT, layer l is correctly identified.
Iteration & Inductive Assumption is satisfied for iteration kt, therefore we recover all layers of n from 0 to k. Thus, for a DAG with an arbitrary number of layers, Stage 4 recovers the full hierarchical sort.
B.7 Time Complexity for NHTS
Theorem 4.6. Given n samples of d vertices generated by identifiable nonlinear ANM, the worst case runtime complexity of Algorithm is upper bounded by (dd
Proof. In Stage 1, NHTS runs (df) marginal independence tests that each have O(n²) complexity. In Stage NHTS runs ((d²) nonparametric regressions and (d2) marginal independence tests, each of which have O complexity. In Stage 3 NHTS runs at most ((d)) conditional independence tests, each of which has O(n") complexity. In the worst case of a fully connected DAG, NHTS goes through (d) steps in Stage 4: in each step of Stage 4, NHTS runs ((d) nonparametric regressions and (d2² marginal independence tests, each of which has (n)) complexity. Therefore, the overall sample complexity of NHTS is O(d'n3))
B.8 Walk-through of NHTS
The following diagram illustrates each stage of NHTS an exemplary 5-node DAG:
In Stage 1, we discover that none of the vertices are PPI relations. In Stage 2, we discover that vertex A is in PP2 relation to vertices B and C, and vertex D is in PP2 relation to vertex E: therefore, and D are our potential roots. In Stage 3, we find that vertex D d-separates its child, vertex E, from vertex A: therefore A is the root vertex. In Stage 4, we regress the unsorted vertices onto A, finding that B,C. are independent of the residual in the first round, then D then E in the last round.
23

---

Stage
Stage 2
Stage 3
Stage 4 [[A][[C]]]
Stage
Stage [[A]].C]ОН]]]
DAG
0
[[]]]
[[][]cC][[]]
A
B
C
8
B
B
B
C
D
D
D
B.9 Time Complexity of NHTS vs RESIT, NoGAM
Theorem 4.7. Consider fully connected DAG (V, E) with nonlinear ANM. Let d [V Let NHTS be the number of nonparametric regressions with covariate set size k E [d 2] run by NHTS NHTS
when sorting V; we similarly define nkk RESIT and rig NoNGAM respectively. Then, 1ig
k; and
NHTS < <1k RESIT NoGAM = k + 1. This implies that for all k R2 k
RESIT = NoGAM
Proof. RESIT and NoGAM both identify leaf vertices in an iterative fashion, regressing each unsorted vertex on the rest of the unsorted vertices; RESIT tests the residual for independence with the covariate set while NoGAM uses the residual for score matching. Therefore, the number of regressions run in both methods in each step equals one plus the covariate set size. Therefore, when the covariate set size is k 4. there are + regressions run.
In the case of a fully directed graph, the first stage of NHTS only runs pairwise regressions with empty conditioning sets. After the first stage, NHTS regresses each unsorted vertex onto all
sorted vertices, finding vertices with independent residuals. Therefore, the number of regressions run is equal to d minus the size of the covariate set. Therefore, when the covariate set size is k > 2 there d regressions run.
24

---

C Proofs and Walkthrough for ED
C.1 Proof of Lemma
Lemma 5.1 (Parent Discovery). Let π be a topological ordering. Zi, 2, such that x(x;)) < (T(x)). Let
7aj Cij U M3j. where Cij { Iy Pa(x::) ak A Fz). Mij {Ix mg E Pa(æ;)) (x:s) T (Tt) < π (x;))) Then, 34 - æj - Is A J. Zzj.
Proof. Note, Ci, blocks all backdoor paths between By and Mij blocks all frontdoor paths between Fa; y, except the direct edge.
We first show the forward direction. If a4 fl. 2. Zan. there must be an direct edge between Tjj and as k(i) < k()), we have T; -> Jzj-
We now show the reverse direction. If + 2j. then there does not exist conditioning set that makes 2; IL Ty, which implies r; Ail 2z [Z
C.2 Hierarchical Version of Edge Discovery
The current implementation of Edge Discovery takes linear hierarchical topological sort as input; this is equivalent to hierarchical topological sort where every layer contains only one vertex. To generalize ED to hierarchical sort, we simply adjust how the algorithm loops over the sort, and which vertices included in which conditioning sets Cij. Maj- We give example of how the latter would change: suppose the algorithm checking whether x; E Pa(æ,) where x; is in layer (I;) 2) and zj is layer 6 (2,) = 6). Mij would be equal the vertices who are parents of 2j that are in layers that are between layer 2 and 6: i.e. Mi equals 24 E Pa(xy) such that 2 Tz (#x) 6. Cij. would stay the same, just being equal to Pa(æ;). Now, generalize the looping part, notice that the linear version of ED loops over indices of the linear topological sort. The hierarchical version of Edge Discovery would loop sort. We give an example of this: suppose the algorithm has found all parents of vertices in layer any
and within layers of the topological
4 lower: the next step would be to find the parents of vertices in layer 5. The algorithm vertex layer 5 as 2. and first finds all parents of xy in the immediately preceding layer, layer 4. The algorithm then finds all parents of zy in the layer preceding layer layer 3. This is essentially the same approach as the linear version, except the size of the conditioning sets and number of conditional tests run is far lower, as the layers provide extra knowledge, limiting which variables can be parents or confounders of other variables.
C.3 Edge Discovery
DAG
Sort
Iter 1-2
Iter 3
Iter 4
A
A
B
B
B
A
C
C
C
C
D
D
D
E
E
E
We first provide thrkuth of ED exemplary 5-node DAG:
25

---

We first use topological ordering algorithm to obtain sort from the DAG: the sort corresponds to fully connected DAG, from which we will prune spurious edges. the first two iterations we find that A B, but B +f C, as B II. C|A. In the third iteration we find that B - D.C D. but A A D, as BI C. In the final iteration, we find that D - F, but A A E, B A C,C A E, given A IL EID, B IL D C IL D. Therefore, we recover the correct set causal edges, removing all spurious edges. We now present a statement of correctness and proof for ED:
Theorem 5.2. Given correct topological ordering of G, Algorithm 3 asymptotically finds correct parent sets (A,),), Vx; E G.
Proof of Correctness
Background We are given a linear topological sort π (V) = [TL, 22,, Fal (where I has no parents). For Za that appears before Tb. let Za,b C..b Ma... where Co,b is the set of potential confounders of aa Tb. and Ma.b is the set of known mediators of I Ty (defined precisely in Lemma 5.1). Let 11 y/- be the value of conditional independence test between x and y.
Proof.
Base Case Iterations (1,2,3,4)
1. Finding parents of x1: by the topological sort, 21 has no parents.
2. Finding parents of Ø2: As #1 is the only possible parent of 22, there are no possible confounders or mediators between and 22, we initialize Zı 0, then by L3
21 A Z = 21 A A D2 22].
2 + 22. Thus, we recover all possible edges between
3. Finding parents of 23: first we check whether 22 + 33, then we check whether I1 23. From iteration 2, we know all edges between [2].
Case 1: if 21 + x2. then it is possible confounder. Therefore, we initialize C2,3 = 21. By topological sort there no mediator between I2 and 23, therefore we initialize M2.3
Then by Lemma 5.1 22 A1 I3l Z2.3 I2 LL 2z|21 < 22 + Z3. If an edge between 22 and 2iz exists, initialize M1,3 accordingly. Note, xi has no parents, so we initialize C1,3
Then, by Lemma 5.1 n azl 71.3 [x1, 22,, x3]]
21 23. We recover all possible edges between
Case 2: if œ1 A x2. there are no possible confounders between 22 and T3. Therefore, we initialize C2.3
By topological sort there is mediator between 22 and 23, therefore we initialize M2,3 = 0. Then by Lemma 5.1 x2 1 Tzl Z2.3 I2 p. V3 < I2 - Z3- Note, 21 has no parents, so we initialize C1.3 = Ø. Note, as D1 A F2, there are no possible mediators between and Zz: we initialize M1.3 = 0. Then, by Lemma 5.1 21 + 33. We recover all possible edges between
21 1 23 Z1. I1 I3 [21, T3)
4. Finding parents of r4: first we check whether whether 21 + Z4. From iteration 3, we know all edges between [21, xz).
24. then whether 22 + then
Case 1: if [Ii, x2, x3] has no edges, then no node causes Ig directly or indirectly, therefore we initialize C2,4 0. There are no possible mediators between 23 and v4, M3.4
Therefore, by Lemma 5.1 23 A sa Za.4 = :3 A J4 O F3 + T4- As [2], 22 33 has no edges,
initialize C2,1 0. Note, x2 A Ø3, so there are no possible mediators between and I4, so we initialize M2,4 0. Then, by Lemma 5.1 Ø2 yLl x4 22,4 22 A I4 + I2 + E4. Note, I has no parents, so we initialize C1,4 = 0. As, does not cause or 23 there are no possible mediators between I1 and F4, therefore we initialize M2,4 0. Then, by Lemma 5.1 1 I4 ZZ,4 21 [21, F3, 22).
node causes 22 directly or indirectly, therefore
1 IA + + Eg. We recover all possible edges between
Case 2: . + B2 is the only edge between the nodes [1), 22, va). then no node causes 23 directly or indirectly, therefore Ca,4 0. By topological sort there is no mediator
26

---

between 23, so we initialize M3,1 0. Then, by Lemma 5.1 23 A In Zs.1 = y
24 = 23 + 24. As 21 - Z2. we initialize C2 = 21. As x2 A 23. there are no possible mediators between 22 and T4, so we initialize M2.4 0. Then, by Lemma 5.1 22 X. x4 22,4 = 22 f I4 [1 O 12 - 24. If an edge exists between 12 and 24, we initialize M1,4 accordingly. As 1 has no parents, we initialize C1, = 0. Then, by Lemma 5.1 #1 A zz4 12,4 + 21 + I4. We recover all possible edges between [x1 22, T3, T/4.
Case 3: if Φ1 + I3 is the only edge between the nodes [T1, 22, 23]. then Ij is the only potential confounder of ny and 24 so we initialize C3,4 21. There are no possible mediators between 23 and 24, so we initialize M3,4 0. Then, by Lemma 5.1 23 y 1 z1| Z3.1 I3 A Fa < 23 + Z4. Note, 22 has no parents, so we initialize 2,4 = 9. Note, 22 A T3, so we initialize M2.4 = A. Then, by Lemma 5.1 22 A 34 72,4 = T2 Al 24 E 22 + T4- If edge exists between 23 and 24, we initialize M14 accordingly. As xı has no parents, we initialize C1,4 D. Then, by Lemma 5.1
31 A F4 11,4 < 21 + 24. We recover all possible edges between T2, 43, Ta].
Case 4: if + 22, + I3 are the only edges between nodes [Iı, 42, 23). then I and 22 cause 23 either indirectly or directly. Therefore, we initialize C3.4 { 22). There are no possible mediators between T3 and T4. so we initialize M3.4 = 0. Then, by Lemma 5.1 23 A 24 Z3.4 = I3 A I4 [1, 22 + 33 + 24. If an edge exists between I3 and 21, we initialize M2,1 accordingly. Note, x1 is parent of 22, so we initialize C2,1 21- Then, by Lemma 5.1 Z2 A xa 22,4 + Z2 + I4. If an edge exists between Z3 and Z4. and/or between I2 and 24, initialize M1,4 accordingly. Note, has no parents, so initialize
C1,4 0. Then by Lemma 5.1 x1 f na Z1,4 + x1 + F4. We recover all possible edges between [21, 23,
Case 5: x1 + I3, x2 + zgare the only edges between nodes 12, ra], then 2ı and I2 cause 23 directly. Therefore, we initialize C3,4 = {II, [2]. There are no possible mediators between T3 and F.g so initialize M3,4 6. Then, by L3 xy A. x. Z,/25,1
T3 AL 21|21, D2 + 23 + I4. If an edge exists between 23 and 24. we initialize M2,4 accordingly. Note, Z2 has no parents, so we initialize C2.1 0. Then, L3 x2 A Ia 22,1 + Z2 + 24. If an edge exists between 23 and I', initialize MI.1 accordingly. Note, 21 has no parents, so initialize C1 = Then by L.3 x fl I.s Z1,4 = 21 T4. We recover all possible edges between [1, 22, T3, ..a.
Case 6: if x1 + I1 + T3, + I3 the only edges between nodes [x 22 ma]. then 2 and x2 cause 23 directly. Therefore, we initialize C3,4 {21, J. There are no possible mediators between i3 and 24, so we initialize M3,4 0. Then, by Lemma 5.1 I' A. x, Z3.1 = I3 A Iy 12 < I3 + I4. If an edge exists between V3 and 21, we initialize M2,4 accordingly. Note, Ti is a parent of ø2, so we initialize C2 = Ø1. Then, by Lemma 5.1 x2 A xz Z2,1 +
+ 24. If an edge exists between 23 and Z1. and/or between za and 21, initialize M1,1 accordingly. Note, Tı has no parents, so initialize
C1,4 = 0. Then by Lemma 5.1 xı JL za Z1, F.).
In + F4. We recover all possible
edges between [x1, T2, T3.
5. Finding parents of r5: first we check whether D4 + Øs, then whether 23 + Fs. then whether 22 + 2s, then whether 21 - 2z.
Iteration A: Inductive Assumption We have recovered the edges between [T] [r,,,..., I.- I]-
Iteration is We now find all nodes in [ar,,..,, [k-]] that cause ær (which yields all edges between
(Ei....,7p)
Base Case Sub-Iteration (1,2)
1. We first check whether ft-1 - TA- We initialize the potential confounders Ck using(r1,...,.t J.)1) The set of mediators M3-1,6 1,k is empty by the topological sort. Then, by Lemma 5.1 Fx-1 > Th b FR-1 11. xnl 7
27

---

2. We now check whether Σk-2 + We initialize the potential confounders Ck-1.k , Only -g-1 can be a mediator: we know whether 21-2 causes 11 and in our previous step we found whether Tk-1 causes Jk- Thus, we initialize M accordingly. Thus, by Lemma 5.1 F + 2k = - 2 A Nz 2. 2,k-
Sub-Iteration Inductive Assumption
We have recovered edges between [z [z...... Tik
Sub-Iteration j - We now find 5 causes T/--
1. For node zj-1 where < j - < k. we obtain Cj-1.& by Iteration k Inductive Assumption and M by Sub-Iteration j Inductive Assumption. Then, by Lemma 5.1 zj-1 - Ig + x--1 _ myl Z
Sub-Iteration j Inductive Assumption is satisfied for j therefore we recover all nodes in [x1,,... .el,,..., -k-- that cause Jy. This satisfies Iteration Inductive Assumption for k, which means we recover all edges between Tu]] Thus, for topological sort of arbitrary length, the algorithm recovers all possible edges.
C.4 Time Complexity for Edge Discovery
Theorem 5.3. Given samples of d vertices generated by model corresponding DAG G, the runlime complexity of ED upper bounded by ((d)n))
Proof. ED checks for the existence of every edge permitted by topological sort ग by running one conditional independence test that has complexity O(n²). In the worst case, there are (dd) possible edges, so the overall complexity is (ddn)).
28

---

D Additional Experiments and Runtimes
D.1 Runtimes for Linear Topological Sort
LHTS
DirectLiNGAM
R Sortabbilliyyy
0.06
0.00
LHTS
DirectLINGAM
R~2-Sortabiiity
100
0.28 0.12
0.05 0.01
0.00 = 0.00
1.37 z 0.49
0.42 = 0.08
0.00 = 0.00
500
0.58 x: 0.21
0.06 : 0.00
0.00 x 0.00
3.34 x 0.86
3.51 0.22
0.02 0.02
700
1.96 0.46
0.10 = 0.02
0.00 = 0.00
0.01
0.00 0.00
1.29
0.03
0.00 = 0.00
20.39
0.04
2000
8.71
0.27 0.04
0.00 = 0.00
Figure 10: Runtimes for linear topological sorts, left: top row; right: bottom row, Figure 5.
D.2 Runtimes for Nonlinear Topological Sort
Noise
NHTS
DLNGM
AcGAM
GES
GRaSP
GSP
R2ST
Uni forro
26.35 12.35
0.10 0.04
03.55 1.42
0.34 0.27
0.03 = 0.02
0.01 = 0.00
0.00 = 0.00
Laplace
22,01 6.39
0.17 0.03
04.14 0.40
0.87 0.10
0.03 1 0.02
0.01 1 0.01
0.00 0.30
Figure 11: Runtimes for nonlinear topological sorts, Figure 6.
D.3 Topological Sorts on Nonlinear Data
3.
3.9
1.1
6.11
Oolen BUDA: 3PA1 DJ00 Suaen SPaP
Figure 12: Performance of NHTS data generated with Gaussian, Laplace, or Uniform noise (left, middle, right columns), the average number of edges set to 2d, 3d, or 4d (top, middle, bottom rows).
In figure 12 we provide additional experiments for nonlinear topological sorts (d = 10, 300); we see that NHTS maintains superior performance even as the density of the underlying graph increases. although the performance gap decreases, especially for laplacian noise, as the graph becomes denser.
29

---

D.4 Edge Pruning on Nonlinear Data
370
Figure 13:
Performance ED on data generated with Gaussian, Laplace, or Uniform noise (left, middle, right columns), the average number of edges set to 2d, 3d, or 4d (top, middle, bottom rows).
In figure 13 we provide additional experiments for edge pruning methods (d = 20, n = 300); ED generally maintains superior performance as the noise distribution is varied and density is increased, although the performance gap decreases for all noise distributions as the graph becomes denser.
D.5 Runtimes for Edge Pruning
Lets
CAM
RESIT
L2090
8.00
Low
CAR
300
8.18 0.24
0.23 9.02
738
1.20
0.02
142
203
281
3 26.0
25.
2:14
1300
3:.6 23.77
043 0.33
0.22
Figure 14: Runtimes for edge pruning: tables correspond to the graphs in Figure 7, from left to right.
30

---

E
Implementation Details
All tests were done in Python. All runtimes were computed locally on Apple M2 Pro Chip, 16 Gb of RAM, with no parallelization.
E.1 Topological Sort for LINGAM
Direct imported from the lingam [45] package, R Rorsort sort was imported from the CausalDisco [28, 29] package, and LHTS was implemented using the Sklearn [24) package. All assets used have a CC-BY 4.0 license. We follow [38] to generate the data for Figure using linear causal mechanisms with randomly drawn coefficient values, plus independent uniform noise. Data is standardized to remove shortcuts [28]. See Github repository for more details, Cutoff values for independence tests were set to a 0.05 for all methods.
E.2 Topological Sort for Nonlinear ANM
GES and GRaSP were imported from the causal-lears [49] package; GSP was imported from the graphical_ wwl.mmelelaeninnnnn lllarniinng [43] package. Direct .NMA was imported from the lingam package [45]. NoGAM was imported from the dodi [14] package. Rt= sort was imported from the CausalDisco package. NHTS and LoSAM were implemented using the kernel ridge regression (KRR) function from the Sklearn package, used independence tests from either the causal-learn package or the dcor [26] package, and a mutual information estimator from the npeet [44] package. All assets used have CC-BY 4.0 license. We follow [16] generate the data used for Figure 6 and Figure 12, using quadratic causal mechanisms with randomly drawn coefficient values, plus independent gaussian, laplace or uniform noise. Features generated with quadratic mechanisms were standardized after being generated to remove shortcuts [28] and to prevent the quadratic mechanisms from driving all values close to (ensuring stability). See Github repository for more details. Note that, to enable fair comparison between NHTS and other topological ordering methods, we implement a version of NHTS that returns a linear topological sort, rather than hierarchical topological sort, by adding only one vertex to the sort in each iteration of its sorting procedure. For NHTS, Stage 2 we used KRR with polynomial kernel, a = degree = 3, coef0 0.1, y 0.01. Cutoff values for independence
and in Stage 4 we used KRR with RBF kernel, tests were set to a = 0.05 for all methods, no cross validation was allowed for any method. Otherwise, default settings were used for all baselines.
E.3
Edge Pruning
Lasso and RESIT were implemented using the sklearn package, hypothesis testing with GAMs was implemented using Bsplines and GLMGam from the statsmodel [33] package. Independence tests used either the causal-learn package or the dcor package, All assets used have a CC-BY 4.0 license. We follow [16] to generate the data used for Figure 7 and Figure 13, using quadratic causal mechanisms with randomly drawn coefficient values, plus independent uniform, gaussian, or laplace noise. Features generated with quadratic mechanisms were standardized after being generated to remove shortcuts [28] and to prevent the quadratic mechanisms from driving all values close (ensuring stability). See Github repository for more details. Cutoff values for independence tests were set to a 0.05 for all methods.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We address all claims made in the abstract and contributions section throughout the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We clearly outline the assumptions and identifiaiility conditions needed for our methods to hold, only claiming aysymptotic correctness when appropriate. See Definition 2.2, Appendix A.I and Appendix B.I.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Each theoretical result includes the necessary assumptions, and links to correct proof in the Appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Experimental procedures are described both Section 6 and Appendix E, and code is released on github.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: Code for the algorithms and data generation is provided in the supplemental material and github link.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: See Section Appendix and the supplemental material github link.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Error bars for experiments for linear topological sorts and edge pruning methods in the main text are standard deviations. Error bars for nonlinear topological sorts and edge pruning experiments not in the main text are whiskers boxplot (1.5 times the IQR).
8. Experiments Compute Resources
32

---

For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Appendix E.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: All parts of the Code of Ethics were followed.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This paper addresses foundational problem in the field of causal discovery, an issue not fundamentally tied to any specific set of applications.
Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: Not applicable.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: See Appendix E.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: No new assets introduced in the paper.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Not applicable.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Not applicable.
33


Reviews:
[{"review_id": "DGgJ8Qi71v", "reviewer": "Reviewer_LEV3", "summary": "The paper mainly focuses on proposing efficient search algorithms for finding the hierarchical sort ordering (linear topological sort) of variables. As mentioned in Section 5, finding such hierarchical orders can significantly improve the efficiency of causal discovery of edges, making the algorithm tractable (traditional algorithms such as PC are exponential). The paper studies two cases: linear (LiNGAMs) and non-parametric, where a complete algorithm based only on path analysis is developed for the linear case, and a combination of path analysis and layer-wise search is developed for the non-parametric case. Both algorithms improve the discovery of hierarchical order.", "questions": "- I'm curious if it's all the results can also be explained using independencies (d-separations) instead of regressions? This allows us to think only in terms of graphs. I guess regressions in the non-parametric setting are equivalent to d-separations, how about the linear case? Are there any independencies that hold in the linear case but not in the non-parametric case?\n- For algorithm 1 (LHTS), is stage 2 really needed? It seems that stage 2 is a special case of stage 3 when mutual ancestors = $\\emptyset$.\n- For experiments, the paper mentions a tradeoff between accuracy and encoded causal information. Would it be more fair to restrict the ordering length (say limit it to some length $k$) and compare the ordering accuracy?", "limitations": "OK", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is well structured and clearly written. The theoretical contributions, including the causal path analysis and corresponding algorithms, are interesting and also important in practice as can be told from the analysis of computational complexity. All results are properly formulated as definitions and theorems and proofs are included in the appendix. Experiments are also conducted and their results are discussed in depth in Section 6. In general, I enjoyed reading it.", "weaknesses": "- In general, I suggest adding more examples to demonstrate the procedure of algorithms, probably for NHTS (Algorithm 2) so that we can see a clear cut between the two stages (root-identification and layer identification).\n- While the authors touched a bit at the beginning of Section 4, non-experts may benefit more if the paper could include additional details about the difference between the linear and non-linear cases (especially how they affect conditional independencies if any). \n- For definition 2.1, it will be great to provide a hierarchical topological sort that cannot be trivially converted to a linear topological sort; that is, we cannot simply add more layers to a hierarchical sort to obtain a linear topological sort.\n- Lemma 4.1 is a little confusing to me: if $x_i$ is a parent of $x_j$, how are PP1 and PP4 possible? $x_i$ must be a direct cause of $x_j$, right? Also, when you say \"$x_i$ and $x_j$\" are not causally related, does it mean that there is no directed edge from $x_i$ to $x_j$ or no directed path? Does \"active path\" mean any unblocked dependency path (backdoor or frontdoor)?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their insightful questions about how our methods work, and suggestions for how to clarify our explanations and results. We respond to Weakness 2 in the general response and address everything else below (following the order of the review).\n\n**Explanatory Examples**\n\nWe thank the Reviewer for their suggestion; in the appendix of our paper, we have added a figure for each algorithm (LHTS, NHTS, ED) with corresponding text that walks through the steps of each algorithm on an exemplary DAG. We hope that this makes it clearer to readers which vertices or edges are being identified in which step, and and easier for readers to understand how the identification occurs. _We include the figure corresponding to NHTS in the PDF attached to the general response and include the corresponding text here:_\n\n\"In Stage 1, NHTS finds that A, B, C, D, and E are all mutually not in PP1 relations. In Stage 2, NHTS discovers that A is in PP2 relation with B and C, while D is in PP2 relation with E: therefore, {A, E} are the candidate root set. In Stage 3, NHTS confirms that A is indeed the root vertex, adding it to the first layer of the hierarchical sort. In the first iteration of Stage 4, NHTS recovers an independent residual when regressing B or C on A; B, C are added to the next layer in the hierarchical sort. Similarly, NHTS recovers an independent residual when regressing D and E in the second and third iteration of Stage 4, respectively: therefore, NHTS recovers the hierarchical topological sort [[A],[B, C],[D],[E]] by the end of Stage 4.\"\n\n\n**Definition 2.1**\n\nTheoretically, any hierarchical topological sort can be converted into a linear topological sort by flattening the hierarchical structure: vertices in the same layer can be placed in an arbitrary order to each other, while vertices in lower layers are placed before \nvertices in upper layers.\n_By \"trivially converted,\" we take the Reviewer to be asking\nwhether given a DAG, the hierarchical topological ordering is unique._\nIndeed, we meant that the hierarchical topological sort in our paper is unique. We have corrected Definition 2.1 to reflect this change:\n\n\"Definition 2.1: A hierarchical topological sort (between a partial and linear order) \nis  a mapping\n$\\pi_L: V \\rightarrow$ {$0,1,\\ldots, |V|-1$}, \nsuch that if $\\text{Pa}(x_i)=\\emptyset$, then $\\pi_L(x_i)=0$, and if $\\text{Pa}(x_i)\\neq \\emptyset$, then $\\pi_L(x_i)$ equals the maximum length of the longest directed path from each root vertex to $x_i$, i.e., $\\pi_L(x_i)= 1 + \\max$ {$\\{\\{\\pi_L(x_j): x_j \\in \\text{Pa}(x_i)\\}}$}.''\n\nWe thank the Reviewer for catching this, and _we are happy to clarify further if this interpretation is incorrect._\n\n**Lemma 4.1**\n\nWe thank the Reviewer for catching the typo in line 217: the sentence now reads \"...where $x_i$ is one of the _potential_ parents of $x_j$...\". Therefore, $x_i$ is not necessarily a parent of $x_j,$: PP1 and PP4 are possible.\n\n We thank the Reviewer for pointing out the unclear notation: $x_i,x_j$ are not causally related if and only if no active path exists between $x_i,x_j$, which includes both directed edges and directed paths. For clarity, the sentence in line 218 now reads: \"PP1) _no active path exists between_ $\\mathbf{x_i,x_j}$...''.\n\nThe Reviewer is correct: we use the term \"active path\" to refer to any unblocked dependency path, which can either be a backdoor path or a frontdoor path.\n\n**Relationship between Results, D-Separations and Regressions**\n\nWe agree that _some of our results may be explained entirely by d-separations, but not all_. We have further clarified the following points in Section 4 of our paper.\n\n- Our edge pruning method ED is entirely constraint-based, and therefore is explainable via d-separations in both the linear and nonlinear case. We refer the Reviewer to our proof of correctness for ED in Appendix C.2 for more detail.\n- Our topological ordering algorithms (LHTS, NHTS) are not fully explainable by d-separations: d-separations only characterize a MEC of DAGs, and different topological sorts may exist within the same MEC. To identify a unique topological sort, we exploit the additive noise assumption: the causal parents of a variable $x_i$ are independent of the noise $\\varepsilon_i$ corresponding $x_i$. In practice, we use this property by running regressions and checking if the residual is independent of the covariates. \n- Further, we note that the conditions under which the residual is independent of the covariates are different in the linear and nonlinear case; we refer the Reviewer to the general response for more detail.\n\n**Stage 2 of LHTS** \n\nWe agree with the Reviewer.\nWe distinguish between\nStages 2 and 3 in Algorithm 1 for explanatory purposes, \nemphasizing that all AP3 relations are discovered in Stage 2, and all AP2 and AP4 relations are discovered in Stage 3.\nThese relations constitute distinct local causal substructures, motivating a conceptual separation between Stages 2 and 3.\nWe have further clarified this point in Line 181 of our revised paper. \n\n**Tradeoff between Accuracy and Encoded Information**\n\nRestricting the ordering length does not impact ordering accuracy. \n_The loss of accuracy of our method comes from classifying vertices into layers according to a cutoff_, instead of pinpointing the 'best' candidate to add to a linear topological sort, which involves selecting the vertex that achieves the best test statistic among the remaining vertices. Both kinds of algorithms are asymptotically consistent. While cutoff-based methods may be slightly more error-prone in low-sample settings, they inherently capture more information about the underlying DAG, reducing the complexity of the edge-pruning stage and making the overall process more efficient when the graph is sparse.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "g3Ypf1qVWQ", "reviewer": "Reviewer_rTf8", "summary": "The paper presents theoretical results about extensions of the partial order induced by a causal DAG and uses these results to propose new constraint-based algorithms for ANMs.\n\n**Edit**: increased rating from 3 to 5, soundness from 1 to 2, and contribution from 2 to 3.\n\n**Edit 2**: increase rating from 5 to 6, after the authors fixed $A_{top}$ calculation; solid paper, but I think the impact of a hybrid causal discovery method like this is limited", "questions": "1. Aren't there just $d \\choose 2$ (i.e., number of entries above the diagonal of a corresponding adjacency matrix) possible edges in a DAG for a given linear order, rather than the $d^2$ claimed on line 303?\n2. Suggestion: Include more explicit theorem statements and proofs for the complexity results.", "limitations": "The authors have adequately described what assumptions their methods require (and hence to which settings the results are limited).", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- Takes a simple idea (which seems original but also somewhat obvious from an order theory perspective) and turns it (creatively, originally) into causal discovery algorithms with contisency guarantees, broad applicability (ANMs), and good identifiability (specific DAG instead of MEC)\n- Very clearly written, as far as grammar, organization, motivation (but importantly, not mathematical notation)\n- Based on the theoretical results, the algorithms have potential to be very significant to the field of causal discovery", "weaknesses": "1. The main (and fatal) weakness is the claims of strong performance in the abstract combined with the inadequate experimental results:\n    1. the abstract makes a claim of \"achieving greater accuracy than current methods\", but then the limited experiments only compare on simulated data (rather than real) to a few closely related algorithms (as opposed to a selection of classic or state-of-the-art methods, such as PC or GRaSP) in settings the authors have already explained are challenging for existing algorithms (very sparse DAGs, rather than a range of sparsities), and even then the proposed algorithm doesn't seem to do especially well. It also seems the NHTS algorithm is missing from the experiments.\n2. A smaller but nonetheless important weekness is notation that contradicts mathematical conventions, making the writing unecessarily difficult:\n    1. consulting introductory texts on partial orders and order theory would help clear up some of the confusion. For example, a topological sort is conventionally a linear extension of a partial order, making the introduced terms \"linear topological sort\" and \"hierarchical topological sort\" confusing. Replacing the former introduced term with just \"topological sort\", \"linear order\", or \"total order\", and replacing the latter introduced term with something that more clearly indicates it is 'between' a partial order and a linear order (i.e., it extends the partial order, but not completely into linear order), would be more natural/conventional and easier to understand.\n    2. the authors seem to use $\\mapsto$ to indicate the domain and image of the ordering functions, but $\\mapsto$ conventionally denotes how a specific element in the domain is mapped to a specific element of the image, hindering precise and easy comprehension.\n    3. other notation in Definition 2.1, such as inconsistent/unexplained indexing of $\\pi$ make the definition harder to understand/not rigorous\n    4.  it's unclear what the difference between $x_j \\dashrightarrow x_i$ (called a directed path) and $x_j \\dashrightarrow \\ldots \\dashrightarrow x_k$ (called a front door path) is.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their detailed comments on how we might improve our notation, and better support our theoretical results. We respond to Weakness 1 in the general response, providing many new experiments that compare our algorithms against many state-of-the-art baselines (CAM, NoGAM, GES, GRaSP, GSP) across a range of sparsities. We find that both our nonlinear topological sorting algorithm NHTS and our edge pruning algorithm ED outperform baselines across all settings, consistently achieving the highest median $A_{top}$ and $F1$. _These strong empirical results provide convincing support for the gains in sample efficiency promised by our theoretical results_. We reply to all other comments below (following the order of the review).\n\n**Unclear Mathematical Notation**\n\nWe thank the Reviewer for identifying unclear notation; we are committed to improving the accessibility of our work, as this allows for communication between disparate academic communities. _We have replaced Definition 2.1 in the main text with the following definitions, incorporating the Reviewer's subpoints 1, 2, and 3:_\n\n- \"Definition 2.1: Consider a given DAG $G = (V,E)$. A topological sort (linear order) is a mapping $\\pi: V \\rightarrow $ {$0,1,\\ldots, |V|-1$}, \nsuch that if $x_i \\in \\text{Pa}(x_j)$, then $x_i$ appears before $x_j$ in the sort: $\\pi(x_i) < \\pi(x_j)$.\"\n- \"Definition 2.2: A layered sort (between a partial and linear order) \nis  a mapping $\\pi_L: V \\rightarrow ${$0,1,\\ldots, |V|-1$}, \nsuch that if $\\text{Pa}(x_i)=\\emptyset$, then $\\pi_L(x_i)=0$, and if $\\text{Pa}(x_i)\\neq \\emptyset$, then $\\pi_L(x_i)$ equals the maximum length of the longest directed path from each root vertex to $x_i$, i.e., $\\pi_L(x_i)= 1 + \\max${$\\pi_L(x_j): x_j \\in \\text{Pa}(x_i)$}.\"\n\n Following a request from Reviewer LEV3, our updated definition of the layered sort enforces uniqueness.  Additionally, we respectfully note that previous work in the literature [1] has used the term \"hierarchical topological ordering\"; however, if the Reviewer strongly believes that this would be confusing to a general audience, we are ok with changing the name to a \"layered sort\" to indicate that the sort is between a partial and linear order.\n \nRegarding subpoint 4: a directed path $x_i \\dashrightarrow x_j$ is the same as a front door path $x_i \\dashrightarrow \\ldots \\dashrightarrow x_j$. For clarity, line 105 now reads \"A frontdoor path is a _directed_ path...''.\n\n**Possible Edges given a Linear Order**\n\nWe thank the Reviewer for catching a typo in Line 304.\nThe sentence now reads \"Edge discovery checks _$O(d^2)$_ possible edges allowed by $\\pi_L$\".\n\n**Explicit Theorems for Complexity Results**\n\nWe thank the Reviewer for the suggestion: _we have added the following theorems for our complexity results to the main text_ (Sections 4, 5, and 6), with the corresponding formal proofs in the appendix. We provide proof sketches below for clarity.\n- _Theorem 1_: Given $n$ samples of $d$ vertices generated by a LiNGAM, the worst case runtime complexity of LHTS is upper bounded by $O(d^3n^2)$. \n  - Proof sketch: in Stage 1, LHTS runs $O(d^2)$ marginal independence tests that each have $O(n^2)$ complexity. In each step of Stage 2 and Stage 3, LHTS runs $O(d^2)$ marginal independence tests, each with $O(n^2)$ complexity. In the worse case of a fully connected DAG, there are $O(d)$ steps in total, across Stage 2 and Stage 3: this is because in each step one layer of the layered sort DAG is identified, and a fully connected DAG has $d$ layers. Therefore, the overall sample complexity of LHTS is $O(d^3n^2)$.\n\n- _Theorem 2_: Given $n$ samples of $d$ vertices generated by an identifiable nonlinear ANM, the worst case runtime complexity of NHTS is upper bounded by $O(d^3n^3)$.\n  - Proof Sketch: In Stage 1, NHTS runs $O(d^2)$ marginal independence tests that each have $O(n^3)$ complexity. In Stage 2, NHTS runs $O(d^2)$ nonparametric regressions and $O(d^2)$ marginal independence tests, each of which has $O(n^3)$ complexity. In Stage 3 NHTS runs at most $O(d^2)$ conditional independence tests, each of which has $O(n^3)$ complexity. In the worst case of a fully connected DAG, NHTS goes through $O(d)$ steps in Stage 4: in each step of Stage 4, NHTS runs $O(d)$ nonparametric regressions and $O(d^2)$ marginal independence tests, each of which has $O(n^3)$ complexity. Therefore, the overall sample complexity of NHTS is $O(d^3n^3)$.\n    \n- _Theorem 3_: Given $n$ samples of $d$ vertices generated by a model corresponding to a DAG $G$, the runtime complexity of ED is upper bounded by $O(d^2n^3)$.\n  - Proof Sketch: ED checks for the existence of every edge permitted by a topological sort $\\pi$ by running one conditional independence test that has complexity $O(n^3)$. In the worst case, there are $O(d^2)$ possible edges, so the overall complexity is $O(d^2n^3)$.\n\n**Citations** \n\n- [1] Wu et. al, Hierarchical Topological Ordering with Conditional\nIndependence Test for Limited Time Series, (2023).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "othmcmjSNW", "reviewer": "Reviewer_LHcu", "summary": "In this paper, the authors present a causal discovery method by firstly determining the order of the causal variables, then determining the existence of edges between any two variables. The experimental results demonstrate the superiority of the proposed method compared to relevant methods.", "questions": "Could the authors further elaborate on Line 205 - 207? It is not quite clear to me why Alg. 1 cannot be used for the non-linear case.\n\nI am happy to adjust my score according to the authors' rebuttal.", "limitations": "No.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "I thank the authors for their detailed clarifications, which address most of my concerns. I increase my score to 5.\n\n------------------------------------\n\nDespite the theoretical results simple, the idea and the method are interesting and somewhat novel. \n\nThe theoretical results seem sensible.", "weaknesses": "1. Lack of necessary discussions: I think there are some similar idea in the literature, such as [1], where they maintain the order of the variables. What is the advantage of this method on [1]? The proposed method should be compared to [1] as well. \n[1] L. Solus, Y. Wang, C. Uhler, and L. Matejovicova. Consistency guarantees for permutation based causal inference algorithms. ArXiv preprint arXiv: 1702.03530 (2017)\n\n2. Lemma 4.1 is confusing. In the condition, it is required that $x_i$ is one of the parents of $x_j$. Why is it possible that $x_i$ and $x_j$ are not causally related?\n\ntypo: \nLine 202: the", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for bringing [1] to our attention. We respond to Question 1 in the general response. We provide clarification on Weakness 1 and 2 below.\n\n**Necessary Discussion of Related Literature**\n\nMethods in [1] tackle a different setting than our paper: they obtain an MEC under the Sparsest Markov Representation assumption, while our methods obtain a unique DAG under the assumption of an ANM. The main contributions of our paper are twofold: 1) our methods _scale to high dimensional datasets by directly building a valid topological sort_, rather than searching across a large permutation space, 2) our methods _improve sample complexity by using smaller, local conditioning sets_ for edge identification. We clarify further below. \n\nMethods in [1] belong to a literature of \"permutation-based\" causal discovery algorithms that search over the space of permutations to find a DAG that achieves an optimal score. These methods build upon the original Sparsest Permutation (SP) algorithm [2], which follows a two step procedure. In the first step, SP constructs a DAG $G_{\\pi} = (V, E)$ for each permutation $\\pi$, where $\\pi(i)\\rightarrow \\pi(j) \\in E \\iff i < j$ and $\\pi(i)$ is not independent of $\\pi(j)$ conditioned on $\\{\\pi_k\\}$ for all $k < j, k \\neq i$. In the second step, SP outputs the set of $G_{\\pi}$ that obtain the minimal number of edges, which corresponds to the MEC. SP exploits the fact that all DAGs in the MEC have the same skeleton structure, implying all valid permutations $\\pi_v$ produce DAGs $G_{\\pi_v}$ that achieve the minimal number of edges; permutations $\\pi_i$ not in the MEC are guaranteed to produce $G_{\\pi_i}$ with strictly larger edge counts [2]. Unfortunately, SP requires a search over all $d!$ permutations to obtain the MEC; _the key contribution of [1] is to introduce a consistent Greedy Sparsest Permutation algorithm (GSP)_ that uses DFS to efficiently explore the space of permutations, again returning the set of $G_{\\pi}$ with minimal edge count (MEC).\n \nWhile GSP improves on SP, it suffers from 1) the need to use heuristics to bound runtime when run on high-dimensional DAGs, and 2) a loss in sample efficiency as conditioning sets grow large. Authors in [1] clarify that it is often necessary in practice to ``bound the search depth $d$ and number of runs $r$ allowed before the algorithm terminates'', reducing the accuracy of GSP as the permutation search space grows. In contrast, our ordering methods do not utilize a heuristic for high-dimensional discovery: _they directly build the correct topological sort, vertex by vertex, entirely avoiding a search procedure._ Additionally, when computing each $G_{\\pi}$, to determine whether the edge $\\pi(i)\\rightarrow \\pi(j)$ exists, GSP runs a CI test that conditions on all vertices prior to $\\pi(j)$ in $\\pi$: in contrast, our edge pruning algorithm ED is sensitive to the sparsity of the underlying DAG, identifying $\\pi(i)\\rightarrow \\pi(j)$ with _conditioning sets that are bounded above by $|\\text{Pa}(\\pi(i))| + |\\text{Pa}(\\pi(j))|$_ (see Step 2 of ED). By leveraging local conditioning sets, ED avoids sample complexity issues suffered by GSP [3].\n\nTo accommodate the Reviewer's suggestion, we have inserted the following two paragraphs at the beginning of our Related Work section:\n\n``Our work is related to two kinds of methods that explicitly leverage the topological structure of DAGs in their discovery procedures: 1) permutation-based approaches and 2) functional causal model-based approaches. \n\nRecent permutation-based approaches, such as SP [2], GSP [1], and GRaSP [4], constitute a family of score-based approaches that utilize permutations for efficient discovery of a MEC. SP searches over the space of variable orderings to find permutations that induce DAGs with minimal edge counts. GSP introduces greedy variants of SP that maintain asymptotic consistency; GRaSP relaxes the assumptions of prior methods to obtain improvements in accuracy. These methods highlight the importance of using permutations for efficient causal discovery, but suffer from poor sample efficiency in high dimensional settings [5].\"\n\nWe have also included various permutation-based algorithms (GSP, GRaSP) as baselines in our new experimental results (see PDF in global response), in order to better evaluate the efficacy of our methods against the discovery literature. _We find that our ordering algorithm NHTS systematically outperforms permutation-based methods._\n    \n**Lemma 4.1**\n\nWe thank the Reviewer for catching a typo in line 217: it is possible that $x_i,x_j$ are not causally related. The sentence now reads \"... where $x_i$ is one of the _potential_ parents of $x_j$ ...\". \n\n**Citations**\n\n- [1] Solus et. al, Consistency guarantees for permutation based causal inference algorithms, (2017).\n- [2] Raskutti et. al, Learning directed acyclic graphs based on sparsest permutations , (2013). \n- [3] Lu et. al, Improving Causal Discovery By Optimal Bayesian Network Learning, (2021). \n- [4] Lam et. al, Greedy Relaxations of the Sparsest Permutation Algorithm, (2022). \n- [5] Niu et. al, Comprehensive Review and Empirical Evaluation of Causal Discovery Algorithms for Numerical Data, (2024).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "33nLX89P0I", "reviewer": "Reviewer_fR7W", "summary": "This paper studies structure learning problem for additive noise model (ANM) in both linear and nonlinear settings. It proposes a hybrid constraint based approach to learn the DAG by leveraging the local ancestral relationships. The algorithm consists of ordering search and edge discovery these two steps. Correctness is shown and simulation is conducted to compare with other approaches.", "questions": "- As the main contribution of the paper, why are the runtime results in the appendix? It is also spurious that the runtime for linear case is slower than the benchmarks in Figure 7; runtime for d=12 is faster than d=8 in Figure 8. There does not seem to be significant improvement empirically;\n- Since theoretically the number of nonparametric regression and CI tests are reduced, is it possible to establish some statistical guarantee and sample complexity dependence on the sparsity, e.g. in-degree?", "limitations": "No limitation is discussed in the paper.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Though ANM is shown to be identifiable for a long time, e.g. by RESIT, the high computational complexity and hardness in nonparametric regression and CI tests stand as roadblock. The finer analysis and exploitation of local structure in the proposal show potential to tackle this task efficiently;\n- The introducing of the proposed method is well-written and easy to follow for researchers working in relevant area.", "weaknesses": "- The main contribution of this work is the exploitation of local structure to reduce the number of nonparametric regression and CI tests. However, despite of the quick discussion below thm 3.7 and 4.5, there is no explicit and formal statement on these to emphasize the contribution, and also comparison with others, e.g. RESIT. \n- The experiments are preliminary. More setups should be considered to demonstrate the superiority of proposal: e.g. different graph types like scale-free graphs, different number of edges, different noise, recovery criterion like F1 for linear setting, more benchmarks like CAM, GSGES, etc. \n- See Questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their suggestions on how we might better emphasize the contributions of our work. We respond to Weakness 2 in the general response, adding to our preliminary experiments by comparing our algorithms against many state-of-the-art baselines (CAM, NoGAM, GES, GRaSP, GSP, etc.). We find that both our nonlinear topological sorting algorithm NHTS and our edge pruning algorithm ED outperform baselines across all settings, consistently achieving the highest median $A_{top}$ and $F1$. _These strong empirical results provide convincing support for the superiority of our proposal_. We give further clarification on all other comments below (following the order of the review).\n\n**Formalizing Computational Complexity Contribution** \n\nWe thank the Reviewer for their suggestion: we now include a formal analysis to compare the number of high-dimensional nonparametric regressions run by our method with baselines such as RESIT and NoGAM in the worst case, i.e. the DAG is fully connected._In particular, we formally show that NHTS runs significantly fewer high-dimensional nonparametric regressions_ to identify the correct topological ordering. We have included the following Theorem in Section 4 (with proof in appendix) of our paper:\n\n\"Theorem  4.6: Consider a fully connected DAG $G=(V,E)$ with nonlinear ANM. Let $d:= |V|$. \nLet $n^\\mathrm{NHTS}_k$ be the number of nonparametric regressions with covariate set size $k\\in[d-2]$ run by NHTS when sorting $V$; we similarly define $n^\\mathrm{RESIT}_k$ and $n^\\mathrm{NoGAM}_k$ respectively. Then, $n^\\mathrm{NHTS}_k = d - k$, and $ n^\\mathrm{RESIT}_k = n^\\mathrm{NoGAM}_k = k+1$. This implies that for all  $k > \\frac{d}{2}$, $n^\\mathrm{NHTS}_k<n^\\mathrm{RESIT}_k = n^\\mathrm{NoGAM}_k$.\"\n\nWe would like to note that it is difficult to analyze the reduction in computational complexity in the general case, as NHTS is sensitive to the specific local graph structure of the underlying data generating process. We leave such a characterization to future work: for now, we demonstrate our method's improved runtime and increased sample efficiency in new experimental results; _NHTS ran $9\\times$ faster than NoGAM with higher median accuracy on a set of 20 dense ER4 graphs (see PDF attached to general response)_.\n\nBelow, we provide a proof sketch to Theorem 4.6:\n\nRESIT [1] and NoGAM [2] both identify leaf vertices in an iterative fashion, regressing each unsorted vertex on the rest of the unsorted vertices; RESIT tests the residual for independence with the covariate set, while NoGAM uses the residual for score matching. Therefore, the number of regressions run in both methods in each step equals\none plus the covariate set size. Therefore, when the covariate set size is $k$, there are $k+1$ regressions run. \n\nIn the case of a fully directed graph, the first stage of NHTS only runs pairwise regressions with empty conditioning sets. After the first stage, NHTS regresses each unsorted vertex onto all sorted vertices, finding vertices with independent residuals. Therefore, the number of regressions run is equal to $d$ minus the size of the covariate set. Therefore, when the covariate set size is $k$, there are $d - k$ regressions run.\n\n**F1 Score** \n\nWe thank the Reviewer for the suggestion. However, we respectfully point out that the Precision subcomponent of the F1 score requires a well-defined notion of false positives and false negatives: _it is unclear how to translate these definitions into the setting of topological sorts_. We instead follow a stream of causal discovery papers [1, 2, 3] that validate their methods using topological divergence measures.\n\n**Runtime Results**\n\nWhile the runtime results of our algorithms are an important contribution, _our work focuses on achieving greater sample efficiency by leveraging local causal substructures to obtain smaller conditioning sets and run fewer high-dimensional nonparametric regressions_, as demonstrated in our new experimental results (see PDF attached to general response). Given the limited space in the manuscript at the time of submission, we chose to use the main text to highlight gains in accuracy demonstrated in preliminary experiments. Overall, our methods were faster than some baselines, while slower than others; we believe that the mentioned statistical anomalies were a result of the chance occurrence of extremely dense DAGs in some trials. We will make these points clear in the additional page afforded to camera-ready submissions.\n \n**Formal Sample Complexity Result** \n\nWe thank the Reviewer for pointing out the need for statistical guarantees of sample complexity dependent on sparsity, and agree that this is a crucial next step. Previous works [4] have successfully demonstrated such statistical guarantees of sample complexity under the assumption of a nonlinear ANM with Gaussian noise, utilizing concentration inequalities for Gaussian random matrices. However, extending these results to nonparametric and non-Gaussian models requires the development of an entirely new set of analytical tools for non-Gaussian matrices, which is _beyond the scope of the present paper._ We \nwill address this challenge in future research.     \n\n**Citations** \n- [1] Peters et. al, Causal Discovery with Continuous Additive Noise Models, (2014).\n- [2] Montagna et. al, Causal Discovery with Score Matching on Additive Models with Arbitrary Noise, (2023). \n- [3] Rolland et. al, Score Matching Enables Causal Discovery of Nonlinear Additive Noise Models, (2022).\n- [4] Zhu et. al, Sample Complexity Bounds for Score-Matching: Causal Discovery and Generative Modeling, (2023).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank the Reviewers for their insightful comments and questions, as they have helped improve the clarity of our paper. We have addressed all raised concerns in this rebuttal, and incorporated the feedback into our manuscript. 

We thank the Reviewers for unanimously acknowledging the novelty and potential utility of our proposed local search approach. In particular, Reviewer rTf8 appreciated how we leveraged a simple idea to develop "causal discovery algorithms with consistency guarantees, broad applicability (ANMs) and good identifiability (specific DAG instead of MEC)." Reviewer fR7W appreciated how "the finer analysis and exploitation of local structure" shows potential to efficiently tackle "the high computational complexity and hardness in nonparametric regression and CI tests." In terms of writing, Reviewers LEV3 and rTf8 both found the paper generally "clearly written", and Reviewer fR7W said it was "easy to follow for researchers working in relevant area."

In this general response, we will provide the following clarifications 
to all Reviewers:

 - Contextualize our paper's **importance to the field of causal discovery** (Reviewer LHcu Weakness 1, Reviewer fR7W Question 1).
 - Provide **additional experimental results** (Reviewer rTf8 Weakness 1, Reviewer fR7W Weakness 2).
- Specify the **difference between linear and nonlinear case** (Reviewer LHcu Question 1, Reviewer LEV3 Weakness 2, Reviewer LEV3 Question 1).

**Importance to the Field of Causal Discovery**

Traditional causal discovery algorithms (PC, GES, GSP, etc.) utilize conditional independence relations to recover a Markov Equivalence Class (MEC) of causal models. To enable the identification of a unique DAG, our methods operate under the assumption of an additive noise model (ANM), leveraging independent residuals from regressions to find a topological ordering, and CI tests to prune edges. Our main contributions are twofold: 1) our nonlinear topological ordering algorithm NHTS _exploits local search to run fewer high-dimensional regressions_ than traditional ANM methods, achieving lower sample complexity than baselines, 2) our constraint-based algorithm ED uses the discovered topological ordering to _nonparametrically prune edges with smaller conditioning sets_ than traditional sparse regression methods, outperforming baselines by improving sample efficiency.

**Additional Experimental Results**

In response to Reviewers fR7W and rTf8 we provide new experimental results (see the attached PDF) that test our nonlinear topological sorting algorithm NHTS and our edge pruning algorithm ED against more benchmarks and in many different settings.
_We find that both NHTS and ED systematically outperform baseline methods, providing 
strong evidence for the superiority of our methods_. 
We are happy to provide additional experimental results upon request. We provide two representative examples of runtime results in the attached PDF. We test our methods in 24 different settings:  two choices of dimensionality ($d=10,d=20$), three types of noise distribution (uniform, gaussian, laplacian),
and
four types of Erdos-Renyi [1] graphs with different sparsities; the expected number of edges either equals to $d$ (ER1), $2d$ (ER2), $3d$ (ER3), or $4d$ (ER4). In each experiment, methods are evaluated on 20 DAGs generated by nonlinear causal mechanisms, with $n=300$. We include NoGAM for comparison as a state-of-the-art ANM-based method [2]. As suggested by Reviewers rTf8, LHcu, and fR7W, we include CAM, GES, GRaSP, and GSP as additional baselines. GES, GRaSP and GSP return only a MEC; to enable a fair comparison, we randomly select one topological ordering permitted by the outputted MEC for evaluation. CAM is excluded from trials with $d=20$ due to
extremely long runtimes, taking over 23 minutes to complete a single trial.
We exclude PC and RESIT since in general they perform much worse than baseline methods [2].

**Difference between Linear and Nonlinear Case**

The reason why LHTS (Algorithm 1) cannot be naively applied is that regressions yield independent residuals under different conditions in the nonlinear case. _LHTS leverages the assumption of linear causal mechanisms when running pairwise regressions in Stage 2, which is insufficient in the nonlinear setting_. For clarity, we demonstrate how LHTS fails to correctly recover causal relationships in an exemplary 3-node DAG with nonlinear causal mechanisms, whereas NHTS (Algorithm 2) succeeds.

Consider DAG $G$ with three vertices $x_1,x_2,x_3$, where $x_1 \rightarrow x_3, x_2 \rightarrow x_3$. The functional causal relationships are nonlinear, given by $x_1 = \varepsilon_1, x_2 = \varepsilon_2, x_3 = x_1x_2 +\varepsilon_3$, where the $\varepsilon$'s are mutually independent noise variables. We focus on whether LHTS or NHTS can recover the parent-child relationship between $x_1$ and $x_3$. Both LHTS and NHTS find that the relationship between $x_1,x_3$ is unknown in Stage 1. 

 In Stage 2, LHTS runs pairwise regressions between $x_1,x_3$ but _incorrectly concludes that $x_1,x_3$ are not in AP3 relation_ because neither pairwise regression provides an independent residual; both parents of $x_3$ must be included in the covariate set for an independent residual to be recovered. In contrast, NHTS correctly constructs the conditioning set $P_{13} = \{x_2\}$ for the pairwise regression of $x_3$ on $x_1$; NHTS recovers an independent residual, _identifying that $x_1,x_3$ are in PP2 relation_, i.e, $x_3$ is a child of $x_1$.

We have incorporated the above explanation into the main text in Section 4. Additionally, as suggested by Reviewer LEV3, for reader clarity we have added figures and text to our paper's appendix that walk through the steps of LHTS, NHTS and ED (see attached PDF for the NHTS figure).

**Citations**
- [1] Erdos et. al, On the evolution of random graphs, (1960).
- [2] Montagna et. al, Causal Discovery with Score Matching on Additive Models with Arbitrary Noise, (2023).

Meta-review: The authors conducted a very detailed rebuttal process and it greatly helped clarify and convince the reviewers. All reviewers are recommending acceptance, although with varying enthusiasm levels. 

I recommend the authors include the "additive noise model" explicitly in the abstract, if not in the title, which currently is not mentioned openly. I would like to especially thank R# rTf8 for closely parsing the rebuttal and notifying the authors of the issue with their experiments on extending PO and also authors for correcting this issue. I am happy to see the review process has made the paper much stronger. The same reviewer first only increased their score to a 5, but after further back and forth with the authors, they increase their score to a 6. Most issues seem to have been resolved assuming that the authors will incorporate all the recommended changes in camera ready. In that sense, another round of review cycle may not be necessary since most of the content is already there.
Decision: Accept (poster)