A polynomial-time algorithm for learning nonparametric causal graphs

Ming Gao Yi Ding Bryon Aragam

Abstract

We establish finite-sample guarantees for a polynomial-time algorithm for learning a nonlinear, nonparametric directed acyclic graphical (DAG) model from data. The analysis is model-free and does not assume linearity, additivity, independent noise, or faithfulness. Instead, we impose a condition on the residual variances that is closely related to previous work on linear models with equal variances. Compared to an optimal algorithm with oracle knowledge of the variable ordering, the additional cost of the algorithm is linear in the dimension $d$ and the number of samples $n$ . Finally, we compare the proposed algorithm to existing approaches in a simulation study.

^†^†Contact: [email protected], [email protected], [email protected]

1 Introduction

Modern machine learning (ML) methods are driven by complex, high-dimensional, and nonparametric models that can capture highly nonlinear phenomena. These models have proven useful in wide-ranging applications including vision, robotics, medicine, and natural language. At the same time, the complexity of these methods often obscure their decisions and in many cases can lead to wrong decisions by failing to properly account for—among other things—spurious correlations, adversarial vulnerability, and invariances (Bottou, 2015; Schölkopf, 2019; Bühlmann, 2018). This has led to a growing literature on correcting these problems in ML systems. A particular example of this that has received widespread attention in recent years is the problem of causal inference, which is closely related to these issues. While substantial methodological progress has been made towards embedding complex methods such as deep neural networks and RKHS embeddings into learning causal graphical models (Huang et al., 2018; Mitrovic et al., 2018; Zheng et al., 2020; Yu et al., 2019; Lachapelle et al., 2019; Zhu and Chen, 2019; Ng et al., 2019), theoretical progress has been slower and typically reserved for particular parametric models such as linear (Chen et al., 2019; Wang and Drton, 2020; Ghoshal and Honorio, 2017, 2018; Loh and Bühlmann, 2014; van de Geer and Bühlmann, 2013; Aragam and Zhou, 2015; Aragam et al., 2015, 2019) and generalized linear models (Park and Raskutti, 2017; Park, 2018).

In this paper, we study the problem of learning directed acyclic graphs (DAGs) from data in a nonparametric setting. Unlike existing work on this problem, we do not require linearity, additivity, independent noise, or faithfulness. Our approach is model-free and nonparametric, and uses nonparametric estimators (kernel smoothers, neural networks, splines, etc.) as “plug-in” estimators. As such, it is agnostic to the choice of nonparametric estimator chosen. Unlike existing consistency theory in the nonparametric setting (Peters et al., 2014; Hoyer et al., 2009; Bühlmann et al., 2014; Rothenhäusler et al., 2018; Huang et al., 2018; Tagasovska et al., 2018; Nowzohour and Bühlmann, 2016), we provide explicit (nonasymptotic) finite sample complexity bounds and show that the resulting method has polynomial time complexity. The method we study is closely related to existing algorithms that first construct a variable ordering (Ghoshal and Honorio, 2017; Chen et al., 2019; Ghoshal and Honorio, 2018; Park and Kim, 2020). Despite this being a well-studied problem, to the best of our knowledge our analysis is the first to provide explicit, simultaneous statistical and computational guarantees for learning nonparametric DAGs.

Contributions

Figure 1(a) illustrates a key motivation for our work: While there exist methods that obtain various statistical guarantees, they lack provably efficient algorithms, or vice versa. As a result, these methods can fail in simple settings. Our focus is on simultaneous computational and statistical guarantees that are explicit and nonasymptotic in a model-free setting. More specifically, our main contributions are as follows:

•

We show that the algorithms of Ghoshal and Honorio (2017) and Chen et al. (2019) rigourously extend to a model-free setting, and provide a method-agnostic analysis of the resulting extension (Theorem 4.1). That is, the time and sample complexity bounds depend on the choice of estimator used, and this dependence is made explicit in the bounds (Section 3.2, Section 4).
•

We prove that this algorithm runs in at most $O(nd^{5})$ time and needs at most $\Omega((d^{2}/\varepsilon)^{1+0pt/2})$ samples (Corollary 4.2). Moreover, the exponential dependence on $d$ can be improved by imposing additional sparsity or smoothness assumptions, and can even be made polynomial (see Section 4 for discussion). This is an expected consequence of our estimator-agnostic approach.
•

We show how existing identifiability results based on ordering variances can be unified and generalized to include model-free families (Theorem 3.1, Section 3.1).
•

We show that greedy algorithms such as those used in the CAM algorithm (Bühlmann et al., 2014) can provably fail to recover an identifiable DAG (Example 5), as shown in Figure 1(a) (Section 3.3).
•

Finally, we run a simulation study to evaluate the resulting algorithm in a variety of settings against seven state-of-the-art algorithms (Section 5).

Our simulation results can be summarized as follows: When implemented using generalized additive models (Hastie and Tibshirani, 1990), our method outperforms most state-of-the-art methods, particularly on denser graphs with hub nodes. We emphasize here, however, that our main contributions lay in the theoretical analysis, specifically providing a polynomial-time algorithm with sample complexity guarantees.

Related work

The literature on learning DAGs is vast, so we focus only on related work in the nonparametric setting. The most closely related line work considers additive noise models (ANMs) (Peters et al., 2014; Hoyer et al., 2009; Bühlmann et al., 2014; Chicharro et al., 2019), and prove a variety of identifiability and consistency guarantees. Compared to our work, the identifiability results proved in these papers require that the structural equations are (a) nonlinear with (b) additive, independent noise. Crucially, these papers focus on (generally asymptotic) statistical guarantees without any computational or algorithmic guarantees. There is also a closely related line of work for bivariate models (Mooij et al., 2016; Monti et al., 2019; Wu and Fukumizu, 2020; Mitrovic et al., 2018) as well as the post-nonlinear model (Zhang and Hyvärinen, 2009). Huang et al. (2018) proposed a greedy search algorithm using an RKHS-based generalized score, and proves its consistency assuming faithfulness. Rothenhäusler et al. (2018) study identifiability of a general family of partially linear models and prove consistency of a score-based search procedure in finding an equivalence class of structures. There is also a recent line of work on embedding neural networks and other nonparametric estimators into causal search algorithms (Lachapelle et al., 2019; Zheng et al., 2020; Yu et al., 2019; Ng et al., 2019; Zhu and Chen, 2019) without theoretical guarantees. While this work was in preparation, we were made aware of the recent work (Park, 2020) that proposes an algorithm that is similar to ours—also based on (Ghoshal and Honorio, 2017) and (Chen et al., 2019)—and establishes its sample complexity for linear Gaussian models. In comparison to these existing lines of work, our focus is on simultaneous computational and statistical guarantees that are explicit and nonasymptotic (i.e. valid for all finite $d$ and $n$ ), for the fully nonlinear, nonparametric, and model-free setting.

Notation

Subscripts (e.g. $X_{j}$ ) will always be used to index random variables and superscripts (e.g. $X_{j}^{(i)}$ ) to index observations. For a matrix $W=(w_{kj})$ , $w_{\cdot j}\in\mathbb{R}^{d}$ is the $j$ th column of $W$ . We denote the indices by $[d]=\{1,\ldots,d\}$ , and frequently abuse notation by identifying the indices $[d]$ with the random vector $X=(X_{1},\ldots,X_{d})$ . For example, nodes $X_{j}$ are interchangeable with their indices $j$ (and subsets thereof), so e.g. $\operatorname{var}(j\,|\,A)$ is the same as $\operatorname{var}(X_{j}\,|\,X_{A})$ .

2 Background

Let $X=(X_{1},\ldots,X_{d})$ be a $d$ -dimensional random vector and $\mathsf{G}=(V,E)$ a DAG where we implicitly assume $V=X$ . The parent set of a node is defined as $\operatorname{pa}_{\mathsf{G}}(X_{j})=\{i\mathrel{\mathop{\mathchar 58\relax}}(i,j)\in E\}$ , or simply $\operatorname{pa}(j)$ for short. A source node is any node $X_{j}$ such that $\operatorname{pa}(j)=\emptyset$ and an ancestral set is any set $A\subset V$ such that $X_{j}\in A\implies\operatorname{pa}(j)\subset A$ . The graph $\mathsf{G}$ is called a Bayesian network (BN) for $X$ if it satisfies the Markov condition, i.e. that each variable is conditionally independent of its non-descendants given its parents. Intuitively, a BN for $X$ can be interpreted as a representation of the direct and indirect relationships between the $X_{j}$ , e.g. an edge $X_{i}\to X_{j}$ indicates that $X_{j}$ depends directly on $X_{i}$ , and not vice versa. Under additional assumptions such as causal minimality and no unmeasured confounding, these arrows may be interpreted causally; for more details, see the surveys (Bühlmann, 2018; Schölkopf, 2019) or the textbooks (Lauritzen, 1996; Koller and Friedman, 2009; Spirtes et al., 2000; Pearl, 2009; Peters et al., 2017).

The goal of structure learning is to learn a DAG $\mathsf{G}$ from i.i.d. observations $X^{(i)}\overset{\text{iid}}{\sim}\mathbb{P}(X)$ . Throughout this paper, we shall exploit the following well-known fact: To learn $\mathsf{G}$ , it suffices to learn a topological sort of $\mathsf{G}$ , i.e. an ordering $\prec$ such that $X_{i}\to X_{j}\implies X_{i}\prec X_{j}$ . A brief review of this material can be found in the supplement.

Equal variances

Recently, a new approach has emerged which was originally cast as an approach to learn equal variance DAGs (Ghoshal and Honorio, 2017; Chen et al., 2019), although it has since been generalized beyond the equal variance case (Ghoshal and Honorio, 2018; Park and Kim, 2020; Park, 2020). An equal variance DAG is a linear structural equation model (SEM) that satisfies

\displaystyle X_{j}=\langle w_{\cdot j},X\rangle+z_{j},\quad\operatorname{var}(z_{j})=\sigma^{2},\quad z_{j}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 4.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 4.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 4.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 4.0mu{\scriptscriptstyle\perp}}}\operatorname{pa}(j),\quad w_{kj}=0\iff k\notin\operatorname{pa}(j)

(1)

for some weights $w_{kj}\in\mathbb{R}$ . Under the model (1), a simple algorithm can learn the graph $\mathsf{G}$ by first learning a topological sort $\prec$ . For these models, we have the following decomposition of the variance:

\displaystyle\operatorname{var}(X_{j})

\displaystyle=\operatorname{var}(\langle w_{\cdot j},X\rangle)+\operatorname{var}(z_{j}).

(2)

Thus, as long as $\operatorname{var}(\langle w_{\cdot j},X\rangle)>0$ , we have $\operatorname{var}(X_{j})>\operatorname{var}(z_{j})$ . It follows that as long as $\operatorname{var}(z_{j})$ does not depend on $j$ , it is possible to identify a source in $\mathsf{G}$ by simply minimizing the residual variances. This is the essential idea behind algorithms based on equal variances in the linear setting (Ghoshal and Honorio, 2017; Chen et al., 2019). Alternatively, it is possible to iteratively identify best sinks by minimizing marginal precisions. Moreover, this argument shows that the assumption of linearity is not crucial, and this idea can readily be extended to ANMs, as in (Park, 2020). Indeed, the crucial assumption in this argument is the independence of the noise $z_{j}$ and the parents $\operatorname{pa}(X_{j})$ ; in the next section we show how these assumptions can be removed altogether.

Layer decomposition of a DAG

Given a DAG $\mathsf{G}$ , define a collection of sets as follows: $L_{0}\mathrel{\mathop{\mathchar 58\relax}}=\emptyset$ , $A_{j}=\cup_{m=0}^{j}L_{m}$ and for $j>0$ , $L_{j}$ is the set of all source nodes in the subgraph $\mathsf{G}[V-A_{j-1}]$ formed by removing the nodes in $A_{j-1}$ . So, e.g., $L_{1}$ is the set of source nodes in $\mathsf{G}$ and $A_{1}=L_{1}$ . This decomposes $\mathsf{G}$ into layers, where each layer $L_{j}$ consists of nodes that are sources in the subgraph $\mathsf{G}[V-A_{j-1}]$ , and $A_{j}$ is an ancestral set for each $j$ . Let $r$ denote the number of “layers” in $\mathsf{G}$ , $L(\mathsf{G})\mathrel{\mathop{\mathchar 58\relax}}=(L_{1},\ldots,L_{r})$ be the corresponding layers. The quantity $r$ effectively measure the depth of a DAG. See Figure 1(b) for an illustration.

Learning $\mathsf{G}$ is equivalent to learning the sets $L_{1},\ldots,L_{r}$ , since any topological sort $\pi$ of $\mathsf{G}$ can be determined from $L(\mathsf{G})$ , and from any sort $\pi$ , the graph $\mathsf{G}$ can be recovered via variable selection. Unlike a topological sort of $\mathsf{G}$ , which may not be unique, the layer decomposition $L(\mathsf{G})$ is always unique. Therefore, without loss of generality, in the sequel we consider the problem of identifying and learning $L(\mathsf{G})$ .

3 Identifiability and algorithmic consequences

This section sets the stage for our main results on learning nonparametric DAGs: First, we show that existing identifiability results for equal variances generalize to a family of model-free, nonparametric distributions. Second, we show that this motivates an algorithm very similar to existing algorithms in the equal variance case. We emphasize that various incarnations of these ideas have appeared in previous work (Ghoshal and Honorio, 2017; Chen et al., 2019; Ghoshal and Honorio, 2018; Park and Kim, 2020; Park, 2020), and our effort in this section is to unify these ideas and show that the same ideas can be applied in more general settings without linearity or independent noise. Once this has been done, our main sample complexity result is presented in Section 4.

3.1 Nonparametric identifiability

In general, a BN for $X$ need not be unique, i.e. $\mathsf{G}$ is not necessarily identifiable from $\mathbb{P}(X)$ . A common strategy in the literature to enforce identifiability is to impose structural assumptions on the conditional distributions $\mathbb{P}(X_{j}\,|\,\operatorname{pa}(j))$ , for which there is a broad literature on identifiability. Our first result shows that identifiability is guaranteed as long as the residual variances $\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))$ do not depend on $j$ . This is a natural generalization of the notion of equality of variances for linear models (Peters and Bühlmann, 2013; Ghoshal and Honorio, 2017; Chen et al., 2019).

Theorem 3.1.

If $\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))\equiv\sigma^{2}$ does not depend on $j$ , then $\mathsf{G}$ is identifiable from $\mathbb{P}(X)$ .

The proof of Theorem 3.1 can be found in the supplement. This result makes no structural assumptions on the local conditional probabilities $\mathbb{P}(X_{j}\,|\,\operatorname{pa}(j))$ . To illustrate, we consider some examples below.

Example 1 (Causal pairs, (Mooij et al., 2016)).

Consider a simple model on two variables: $X\to Y$ with $\mathbb{E}\operatorname{var}(Y\,|\,X)=\operatorname{var}(X)$ . Then as long as $\mathbb{E}[Y\,|\,X]$ is nonconstant, Theorem 3.1 implies the causal order is identifiable. No additional assumptions on the noise or functional relationships are necessary.

Example 2 (Binomial models, (Park and Raskutti, 2017)).

Assume $X_{j}\in\{0,1\}$ and $X_{j}=\operatorname{Ber}(f_{j}(\operatorname{pa}(j)))$ with $f_{j}(\operatorname{pa}(j))\in[0,1]$ . Then Theorem 3.1 implies that if $\mathbb{E}f_{j}(\operatorname{pa}(j))(1-f_{j}(\operatorname{pa}(j)))\equiv\sigma^{2}$ does not depend on $j$ , then $\mathsf{G}$ is identifiable.

Example 3 (Generalized linear models).

The previous example can of course be generalized to arbitrary generalized linear models: Assume $\mathbb{P}[X_{j}\,|\,\operatorname{pa}(j)]\propto\exp(X_{j}\theta_{j}-K(\theta_{j}))$ , where $\theta_{j}=f_{j}(\operatorname{pa}(j))$ and $K(\theta_{j})$ is the partition function. Then Theorem 3.1 implies that if $\mathbb{E}[K^{\prime\prime}(f_{j}(\operatorname{pa}(j)))]\equiv\sigma^{2}$ does not depend on $j$ , then $\mathsf{G}$ is identifiable.

Example 4 (Additive noise models, (Peters et al., 2014)).

Finally, we observe that Theorem 3.1 generalizes existing results for ANMs: In an ANM, we have $X_{j}=f_{j}(\operatorname{pa}(j))+z_{j}$ with $z_{j}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 4.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 4.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 4.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 4.0mu{\scriptscriptstyle\perp}}}\operatorname{pa}(j)$ . If $\operatorname{var}(z_{j})=\sigma^{2}$ , then an argument similar to (2) shows that ANMs with equal variances are identifiable. Theorem 3.1 applies to more general additive noise models $X_{j}=f_{j}(\operatorname{pa}(j))+g_{j}(\operatorname{pa}(j))^{1/2}z_{j}$ with heteroskedastic, uncorrelated (i.e. not necessarily independent) noise.

Unequal variances

Early work on this problem focused on the case of equal variances (Ghoshal and Honorio, 2017; Chen et al., 2019), as we have done here. This assumption illustrates the main technical difficulties in proving identifiability, and it is well-known by now that equality of variances is not necessary, and a weaker assumption that allows for heterogeneous residual variances suffices in special cases (Ghoshal and Honorio, 2018; Park and Kim, 2020). Similarly, the extension of Theorem 3.1 to such heterogeneous models is straightforward, and omitted for brevity; see Appendix B.1 in the supplement for additional discussion and simulations. In the sequel, we focus on the case of equality for simplicity and ease of interpretation.

3.2 A polynomial-time algorithm

The basic idea behind the top-down algorithm proposed in (Chen et al., 2019) can easily be extended to the setting of Theorem 3.1, and is outlined in Algorithm 1. The only modification is to replace the error variances $\operatorname{var}(z_{j})=\sigma^{2}$ from the linear model (1) with the corresponding residual variances (i.e. $\mathbb{E}\operatorname{var}(X_{\ell}\,|\,S_{j})$ ), which are well-defined for any $\mathbb{P}(X)$ with finite second moments.

A natural idea to translate Algorithm 1 into an empirical algorithm is to replace the residual variances with an estimate based on the data. One might then hope to use similar arguments as in the linear setting to establish consistency and bound the sample complexity. Perhaps surprisingly, this does not work unless the topological sort of $\mathsf{G}$ is unique. When there is more than one topological sort, it becomes necessary to uniformly bound the errors of all possible residual variances—and in the worst case there are exponentially many ( $d2^{d-1}$ to be precise) possible residual variances. The key issue is that the sets $S_{j}$ in Algorithm 1 are random (i.e. data-dependent), and hence unknown in advance. This highlights a key difference between our algorithm and existing work for linear models such as (Ghoshal and Honorio, 2017, 2018; Chen et al., 2019; Park, 2020): In our setting, the residual variances cannot be written as simple functions of the covariance matrix $\Sigma\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}X\!X^{T}$ , which simplifies the analysis for linear models considerably. Indeed, although the same exponential blowup arises for linear models, in that case consistent estimation of the covariance matrix $\Sigma\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}X\!X^{T}$ provides uniform control over all possible residual variances (e.g., see Lemma 6 in (Chen et al., 2019)). In the nonparametric setting, this reduction no longer applies.

To get around this technical issue, we modify Algorithm 1 to learn $\mathsf{G}$ one layer $L_{j}$ at a time, as outlined in Algorithm 2 (see Section 2 for details on $L_{j}$ ). As a result, we need only estimate $\sigma_{\ell j}^{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})$ , which involves regression problems with at most $|A_{j}|$ nodes. We use the plug-in estimator (3) for this, although more sophisticated estimators are available (Doksum et al., 1995; Robins et al., 2008). This also necessitates the use of sample splitting in Step 3(a) of Algorithm 2, which is necessary for the theoretical arguments but not needed in practice.

The overall computational complexity of Algorithm 2, which we call NPVAR, is $O(ndrT)$ , where $T$ is the complexity of computing each nonparametric regression function $\widehat{f}_{\ell j}$ . For example, if a kernel smoother is used, $T=O(d^{3})$ and thus the overall complexity is $O(nrd^{4})$ . For comparison, an oracle algorithm that knows the true topological order of $\mathsf{G}$ in advance would still need to compute $d$ regression functions, and hence would have complexity $O(dT)$ . Thus, the extra complexity of learning the topological order is only $O(nr)=O(nd)$ , which is linear in the dimension and the number of samples. Furthermore, under additional assumptions on the sparsity and/or structure of the DAG, the time complexity can be reduced further, however, our analysis makes no such assumptions.

Algorithm 1 Population algorithm for learning nonparametric DAGs

Set $S_{0}=\emptyset$ and for $j=0,1,2,\ldots$ , let

\displaystyle k_{j}

\displaystyle=\operatorname*{arg\,min}_{\ell\notin S_{j}}\mathbb{E}\operatorname{var}(X_{\ell}\,|\,S_{j}),\qquad S_{j+1}=S_{j}\cup\{k_{j}\}.

2.

Return the DAG $\mathsf{G}$ that corresponds to the topological sort $(k_{1},\ldots,k_{d})$ .

Algorithm 2 NPVAR algorithm

Input: $X^{(1)},\ldots,X^{(n)}$ , $\eta>0$ .

1.

Set $\widehat{L}_{0}=\emptyset$ , $\widehat{\sigma}_{\ell 0}^{2}=\widehat{\operatorname{var}}(X_{\ell})$ , $k_{0}=\operatorname*{arg\,min}_{\ell}\widehat{\sigma}_{\ell 0}^{2}$ , $\widehat{\sigma}_{0}^{2}=\sigma_{k_{0}0}^{2}$ .
2.

Set $\widehat{L}_{1}\mathrel{\mathop{\mathchar 58\relax}}=\{\ell\mathrel{\mathop{\mathchar 58\relax}}|\widehat{\sigma}_{\ell 0}^{2}-\widehat{\sigma}_{0}^{2}|<\eta\}$ .

For $j=2,3,\ldots$ :

(a)

Randomly split the $n$ samples in half and let $\widehat{A}_{j}\mathrel{\mathop{\mathchar 58\relax}}=\cup_{m=1}^{j}\widehat{L}_{m}$ .
(b)

For each $\ell\notin\widehat{A}_{j}$ , use the first half of the sample to estimate $f_{\ell j}(X_{\widehat{A}_{j}})=\mathbb{E}[X_{\ell}\,|\,\widehat{A}_{j}]$ via a nonparametric estimator $\widehat{f}_{\ell j}$ .

(c)

For each $\ell\notin\widehat{A}_{j}$ , use the second half of the sample to estimate the residual variances via the plug-in estimator

\displaystyle\widehat{\sigma}_{\ell j}^{2}

\displaystyle=\frac{1}{n/2}\sum_{i=1}^{n/2}(X_{\ell}^{(i)})^{2}-\frac{1}{n/2}\sum_{i=1}^{n/2}\widehat{f}_{\ell j}(X_{\widehat{A}_{j}}^{(i)})^{2}.

(3)

(d)

Set $k_{j}=\operatorname*{arg\,min}_{\ell\notin\widehat{A}_{j}}\widehat{\sigma}_{\ell j}^{2}$ and $\widehat{L}_{j+1}=\{\ell\mathrel{\mathop{\mathchar 58\relax}}|\widehat{\sigma}_{\ell j}^{2}-\widehat{\sigma}_{k_{j}j}^{2}|<\eta,\,\ell\notin\widehat{A}_{j}\}$ .

4.

Return $\widehat{L}=(\widehat{L}_{1},\ldots,\widehat{L}_{\widehat{r}})$ .

3.3 Comparison to existing algorithms

Compared to existing algorithms based on order search and equal variances, NPVAR applies to more general models without parametric assumptions, independent noise, or additivity. It is also instructive to make comparisons with greedy score-based algorithms such as causal additive models (CAM, (Bühlmann et al., 2014)) and greedy DAG search (GDS, (Peters and Bühlmann, 2013)). We focus here on CAM since it is more recent and applies in nonparametric settings, however, similar claims apply to GDS as well.

CAM is based around greedily minimizing the log-likelihood score for additive models with Gaussian noise. In particular, it is not guaranteed to find a global minimizer, which is as expected since it is based on a nonconvex program. This is despite the global minimizer—if it can be found—having good statistical properties. The next example shows that, in fact, there are identifiable models for which CAM will find the wrong graph with high probability.

Example 5.

Consider the following three-node additive noise model with $z_{j}\sim\mathcal{N}(0,1)$ :

\displaystyle\begin{aligned} X_{1}&=z_{1},\\ X_{2}&=g(X_{1})+z_{2},\\ X_{3}&=g(X_{1})+g(X_{2})+z_{3}.\end{aligned}

(4)

In the supplement (Appendix D), we show the following: There exist infinitely many nonlinear functions $g$ for which the CAM algorithm returns an incorrect order under the model (4). This is illustrated empirically in Figure 1(a) for the nonlinearities $g(u)=\text{sgn}(u)|u|^{1.4}$ and $g(u)=\sin u$ . In each of these examples, the model satisfies the identifiability conditions for CAM as well as the conditions required in our work.

We stress that this example does not contradict the statistical results in Bühlmann et al. (2014): It only shows that the algorithm may not find a global minimizer and as a result, returns an incorrect variable ordering. Correcting this discrepancy between the algorithmic and statistical results is a key motivation behind our work. In the next section, we show that NPVAR provably learns the true ordering—and hence the true DAG—with high probability.

4 Sample complexity

Our main result analyzes the sample complexity of NPVAR (Algorithm 2). Recall the layer decomposition $L(\mathsf{G})$ from Section 2 and define $d_{j}\mathrel{\mathop{\mathchar 58\relax}}=|A_{j}|$ . Let $f_{\ell j}(X_{A_{j}})=\mathbb{E}[X_{\ell}\,|\,A_{j}]$ .

Condition 1 (Regularity).

For all $j$ and all $\ell\notin A_{j}$ , (a) $X_{j}\in[0,1]$ , (b) $f_{\ell j}\mathrel{\mathop{\mathchar 58\relax}}[0,1]^{d_{j}}\to[0,1]$ , (c) $f_{\ell j}\in L^{\infty}([0,1]^{d_{j}})$ , and (d) $\operatorname{var}(X_{\ell}\,|\,A_{j})\leq\zeta_{0}<\infty$ .

These are the standard regularity conditions from the literature on nonparametric statistics (Györfi et al., 2006; Tsybakov, 2009), and can be weakened (e.g. if the $X_{j}$ and $f_{\ell j}$ are unbounded, see (Kohler et al., 2009)). We impose these stronger assumptions in order to simplify the statements and focus on technical details pertinent to graphical modeling and structure learning. The next assumption is justified by Theorem 3.1, and as we have noted, can also be weakened.

Condition 2 (Identifiability).

$\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))\equiv\sigma^{2}$ does not depend on $j$ .

Our final condition imposes some basic finiteness and consistency requirements on the chosen nonparametric estimator $\widehat{f}$ , which we view as a function for estimating $\mathbb{E}[Y\,|\,Z]$ from an arbitrary distribution over the pair $(Y,Z)$ .

Condition 3 (Estimator).

The nonparametric estimator $\widehat{f}$ satisfies (a) $\mathbb{E}[Y\,|\,Z]\in L^{\infty}\implies\widehat{f}\in L^{\infty}$ and (b) $\mathbb{E}_{\widehat{f}}\|\widehat{f}(Z)-\mathbb{E}[Y\,|\,Z]\|_{2}^{2}\to 0$ .

This is a mild condition that is satisfied by most popular estimators including kernel smoothers, nearest neighbours, and splines, and in particular, Condition 3(a) is only used to simplify the theorem statement and can easily be relaxed.

Theorem 4.1.

Assume Conditions 1-3. Let $\Delta_{j}>0$ be such that $\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})>\sigma^{2}+\Delta_{j}$ for all $\ell\notin A_{j}$ and define $\Delta\mathrel{\mathop{\mathchar 58\relax}}=\inf_{j}\Delta_{j}$ . Let $\delta^{2}\mathrel{\mathop{\mathchar 58\relax}}=\sup_{\ell,j}\mathbb{E}_{\widehat{f}_{\ell j}}\|f_{\ell j}(X_{A_{j}})-\widehat{f}_{\ell j}(X_{A_{j}})\|_{2}^{2}$ . Then for any $\delta\sqrt{d}<\eta<\Delta/2$ ,

\displaystyle\mathbb{P}(\widehat{L}=L(\mathsf{G}))

\displaystyle\gtrsim 1-\frac{\delta^{2}}{\eta^{2}}rd

(5)

Once the layer decomposition $L(\mathsf{G})$ is known, the graph $\mathsf{G}$ can be learned via standard nonlinear variable selection methods (see Appendix A in the supplement).

A feature of this result is that it is agnostic to the choice of estimator $\widehat{f}$ , as long as it satisfies Condition 3. The dependence on $\widehat{f}$ is quantified through $\delta^{2}$ , which depends on the sample size $n$ and represents the rate of convergence of the chosen nonparametric estimator. Instead of choosing a specific estimator, Theorem 4.1 is stated so that it can be applied to general estimators. As an example, suppose each $f_{\ell j}$ is Lipschitz continuous and $\widehat{f}$ is a standard kernel smoother. Then

\displaystyle\mathbb{E}_{\widehat{f}_{\ell j}}\|f_{\ell j}(X_{L_{j}})-\widehat{f}_{\ell j}(X_{L_{j}})\|_{2}^{2}\leq\delta^{2}\lesssim n^{-\tfrac{2}{2+0pt}}.

Thus we have the following special case:

Corollary 4.2.

Assume each $f_{\ell j}$ is Lipschitz continuous. Then $\widehat{L}$ can be computed in $O(nd^{5})$ time and $\mathbb{P}(\widehat{L}=L(\mathsf{G}))\geq 1-\varepsilon$ as long as $n=\Omega((rd/(\eta^{2}\varepsilon))^{1+d/2})$ .

This is the best possible rate attainable by any algorithm without imposing stronger regularity conditions (see e.g. §5 in (Györfi et al., 2006)). Furthermore, $\delta^{2}$ can be replaced with the error of an arbitrary estimator of the residual variance itself (i.e. something besides the plug-in estimator (3)); see Proposition C.4 in Appendix C for details.

To illustrate these results, consider the problem of finding the direction of a Markov chain $X_{1}\to X_{2}\to\cdots\to X_{d}$ whose transition functions $\mathbb{E}[X_{j}\,|\,X_{j-1}]$ are each Lipschitz continuous. Then $r=d$ , so Corollary 4.2 implies that $n=\Omega((d^{2}/(\eta\sqrt{\varepsilon}))^{1+d/2})$ samples are sufficient to learn the order—and hence the graph as well as each transition function—with high probability. Since $r=d$ for any Markov chain, this particular example maximizes the dependence on $d$ ; at the opposite extreme a bipartite graph with $r=2$ would require only $n=\Omega((\sqrt{d}/(\eta\sqrt{\varepsilon}))^{1+d/2})$ . In these lower bounds, it is not necessary to know the type of graph (e.g. Markov chain, bipartite) or the depth $r$ .

Choice of $\eta$

The lower bound $\eta>\delta\sqrt{d}$ is not strictly necessary, and is only used to simplify the lower bound in (5). In general, taking $\eta$ sufficiently small works well in practice. The main tradeoff in choosing $\eta>0$ is computational: A smaller $\eta$ may lead to “splitting” one of the layers $L_{j}$ . In this case, NPVAR still recovers the structure correctly, but the splitting results in redundant estimation steps in Step 3 (i.e. instead of estimating $L_{j}$ in one iteration, it takes multiple iterations to estimate correctly). The upper bound, however, is important: If $\eta$ is too large, then we may include spurious nodes in the layer $L_{j}$ , which would cause problems in subsequent iterations.

Nonparametric rates

Theorem 4.1 and Corollary 4.2 make no assumptions on the sparsity of $\mathsf{G}$ or smoothness of the mean functions $\mathbb{E}[X_{\ell}\,|\,A_{j}]$ . For this reason, the best possible rate for a naïve plug-in estimator of $\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})$ is bounded by the minimax rate for estimating $\mathbb{E}[X_{\ell}\,|\,A_{j}]$ . For practical reasons, we have chosen to focus on an agnostic analysis that does not rely on any particular estimator. Under additional sparsity and smoothness assumptions, these rates can be improved, which we briefly discuss here.

For example, by using adaptive estimators such as RODEO (Lafferty and Wasserman, 2008) or GRID (Giordano et al., 2020), the sample complexity will depend only on the sparsity of $f_{\ell j}(X_{A_{j}})$ , i.e. $d^{*}=\max_{j}\max_{\ell\notin A_{j}}|\{k\in A_{j}\mathrel{\mathop{\mathchar 58\relax}}\partial_{k}f_{\ell j}\neq 0\}|$ , where $\partial_{k}$ is the $k$ th partial derivative. Another approach that does not require adaptive estimation is to assume $|L_{j}|\leq w$ and define $r^{*}\mathrel{\mathop{\mathchar 58\relax}}=\sup\{|i-j|\mathrel{\mathop{\mathchar 58\relax}}e=(e_{1},e_{2})\in E,e_{1}\in L_{i},e_{2}\in L_{j}\}$ . Then $\delta^{2}\asymp n^{-2/(2+wr^{*})}$ , and the resulting sample complexity depends on $wr^{*}$ instead of $d$ . For a Markov chain with $w=r^{*}=1$ this leads to a substantial improvement.

Instead of sparsity, we could impose stronger smoothness assumptions: Let $\beta_{*}$ denote the smallest Hölder exponent of any $f_{\ell j}$ . Then if $\beta_{*}\geq d/2$ , then one can use a one-step correction to the plug-in estimator (3) to obtain a root- $n$ consistent estimator of $\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})$ (Robins et al., 2008; Kandasamy et al., 2015). Another approach is to use undersmoothing (Doksum et al., 1995). In this case, the exponential sample complexity improves to polynomial sample complexity. For example, in Corollary 4.2, if we replace Lipschitz with the stronger condition that $\beta_{*}\geq d/2$ , then the sample complexity improves to $n=\Omega(rd/(\eta^{2}\varepsilon))$ .

5 Experiments

Finally, we perform a simulation study to compare the performance of NPVAR against state-of-the-art methods for learning nonparametric DAGs. The algorithms are: RESIT (Peters et al., 2014), CAM (Bühlmann et al., 2014), EqVar (Chen et al., 2019), NOTEARS (Zheng et al., 2020), GSGES (Huang et al., 2018), PC (Spirtes and Glymour, 1991), and GES (Chickering, 2003). In our implementation of NPVAR, we use generalized additive models (GAMs) for both estimating $\widehat{f}_{\ell j}$ and variable selection. One notable detail is our implementation of EqVar, which we adapted to the nonlinear setting by using GAMs instead of subset selection for variable selection (the order estimation step remains the same). Full details of the implementations used as well as additional experiments can be found in the supplement. Code implementing the NPVAR algorithm is publicly available at https://github.com/MingGao97/NPVAR.

We conducted a series of simulation on different graphs and models, comparing the performance in both order recovery and structure learning. Due to space limitations, only the results for structure learning in the three most difficult settings are highlighted in Figure 2. These experiments correspond to non-sparse graphs with non-additive dependence given by either a Gaussian process (GP) or a generalized linear model (GLM):

•

Graph types. We sampled three families of DAGs: Markov chains (MC), Erdös-Rényi graphs (ER), and scale-free graphs (SF). For MC graphs, there are exactly $d$ edges, whereas for ER and SF graphs, we sample graphs with $kd$ edges on average. This is denoted by ER4/SF4 for $k=4$ in Figure 2. Experiments on sparser DAGs can be found in the supplement.
•

Probability models. For the Markov chain models, we used two types of transition functions: An additive sine model with $\mathbb{P}(X_{j}\,|\,X_{j-1})=\mathcal{N}(\sin(X_{j-1}),\sigma^{2})$ and a discrete model (GLM) with $X_{j}\in\{0,1\}$ and $\mathbb{P}(X_{j}\,|\,X_{j-1})\in\{p,1-p\}$ . For the ER and SF graphs, we sampled $\mathbb{E}[X_{j}\,|\,\operatorname{pa}(j)]$ from both additive GPs (AGP) and non-additive GPs (NGP).

Full details as well as additional experiments on order recovery, additive models, sparse graphs, and misspecified models can be found in the supplement (Appendix E).

Structure learning

To evaluate overall performance, we computed the structural Hamming distance (SHD) between the learned DAG and the true DAG. SHD is a standard metric used for comparison of graphical models. According to this metric, the clear leaders are NPVAR, EqVar, and CAM. Consistent with existing results, existing methods tend to suffer as the edge density and dimension of the graphs increase, however, NPVAR is more robust in these settings. Surprisingly, the CAM algorithm remains quite competitive for non-additive models, although both EqVar and NPVAR clearly outperform CAM. On the GLM model, which illustrates a non-additive model with non-additive noise, EqVar and NPVAR performed the best, although PC showed good performance with $n=1000$ samples. Both CAM and RESIT terminated with numerical issues on the GLM model.

These experiments serve to corroborate our theoretical results and highlight the effectiveness of the NPVAR algorithm, but of course there are tradeoffs. For example, algorithms such as CAM which exploit sparse and additive structure perform very well in settings where sparsity and additivity can be exploited, and indeed outperform NPVAR in some cases. Hopefully, these experiments can help to shed some light on when various algorithms are more or less effective.

Misspecification and sensitivity analysis

We also considered two cases of misspecification: In Appendix B.1, we consider an example where Condition 2 fails, but NPVAR still successfully recovers the true ordering. This experiment corroborates our claims that this condition can be relaxed to handle unequal residual variances. We also evaluated the performance of NPVAR on linear models as in (1), and in all cases it was able to recover the correct ordering.

6 Discussion

In this paper, we analyzed the sample complexity of a polynomial-time algorithm for estimating nonparametric causal models represented by a DAG. Notably, our analysis avoids many of the common assumptions made in the literature. Instead, we assume that the residual variances are equal, similar to assuming homoskedastic noise in a standard nonparametric regression model. Our experiments confirm that the algorithm, called NPVAR, is effective at learning identifiable causal models and outperforms many existing methods, including several recent state-of-the-art methods. Nonetheless, existing algorithms such as CAM are quite competitive and apply in settings where NPVAR does not.

We conclude by discussing some limitations and directions for future work. Although we have relaxed many of the common assumptions made in the literature, these assumptions have been replaced by an assumption on the residual variances that may not hold in practice. An interesting question is whether or not there exist provably polynomial-time algorithms for nonparametric models in under less restrictive assumptions. Furthermore, although the proposed algorithm is polynomial-time, the worst-case $O(d^{5})$ dependence on the dimension is of course limiting. This can likely be reduced by developing more efficient estimators of the residual variance that do not first estimate the mean function. This idea is common in the statistics literature, however, we are not aware of such estimators specifically for the residual variance (or other nonlinear functionals of $\mathbb{P}(X)$ ). Furthermore, our general approach can be fruitfully applied to study various parametric models that go beyond linear models, for which both computation and sample efficiency would be expected to improve. These are interesting directions for future work.

Acknowledgements

We thank the anonymous reviewers for valuable feedback, as well as Y. Samuel Wang and Edward H. Kennedy for helpful discussions. B.A. acknowledges the support of the NSF via IIS-1956330 and the Robert H. Topel Faculty Research Fund. Y.D.’s work has been partially supported by the NSF (CCF-1439156, CCF-1823032, CNS-1764039).

Appendix A Reduction to order search

The fact that DAG learning can be reduced to learning a topological sort is well-known. For example, this fact is the basis of exact algorithms for score-based learning based on dynamic programming Silander and Myllymaki (2006); Ott et al. (2004); Perrier et al. (2008); Singh and Moore (2005) as well as recent algorithms for linear models (Chen et al., 2019; Ghoshal and Honorio, 2017, 2018). See also (Shojaie and Michailidis, 2010). This fact has also been exploited in the overdispersion scoring model developed by Park and Raskutti (2017) as well as for nonlinear additive models (Bühlmann et al., 2014). In fact, more can be said: Any ordering defines a minimal I-map of $\mathbb{P}(X)$ via a simple iterative algorithm (see §3.4.1, Algorithm 3.2 in (Koller and Friedman, 2009)), and this minimal I-map is unique as long as $\mathbb{P}(X)$ satisfies the intersection property. This is guaranteed, for example, if $\mathbb{P}(X)$ has a positive density, but holds under weaker conditions (see (Peters, 2015) for necessary and sufficient conditions assuming $\mathbb{P}(X)$ has a density and (Dawid, 1980), Theorem 7.1, for the general case). This same algorithm can then be used to reconstruct the true DAG $\mathsf{G}$ from the true ordering $\prec$ . As noted in Section 2, a further reduction can be obtained by considering the layer decomposition $L(\mathsf{G})$ , from which all topological orders $\prec$ of $\mathsf{G}$ can be deduced.

Once the ordering is known, existing nonlinear variable selection methods (Lafferty and Wasserman, 2008; Rosasco et al., 2013; Giordano et al., 2020; Miller and Hall, 2010; Bertin and Lecué, 2008; Comminges and Dalalyan, 2011) suffice to learn the parent sets $\operatorname{pa}(j)$ and hence the graph $\mathsf{G}$ . More specifically, given an order $\prec$ , to identify $\operatorname{pa}(j)$ , let $f(S_{j})\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}[X_{j}\,|\,S_{j}]$ , where $S_{j}\mathrel{\mathop{\mathchar 58\relax}}=\{X_{k}\mathrel{\mathop{\mathchar 58\relax}}X_{k}\prec X_{j}\}$ . The parent set of $X_{j}$ is given by the active variables in this conditional expectation, i.e. $\operatorname{pa}(j)=\{k\mathrel{\mathop{\mathchar 58\relax}}\partial_{k}f\neq 0\}$ , where $\partial_{k}$ is the partial derivative of $f$ with respect to the $k$ th argument.

In our experiments, we use exactly this procedure to learn $\mathsf{G}$ from the order $\prec$ , based on the data. Specifically, we use generalized additive models, similar to the pruning step in (Bühlmann et al., 2014). See Appendix E for more details.

Appendix B Proof of Theorem 3.1

The key lemma is the following, which is easy to prove for additive noise models via (2), and which we show holds more generally in non-additive models:

Lemma B.1.

Let $A\subset V$ be an ancestral set in $\mathsf{G}$ . If $\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))\equiv\sigma^{2}$ does not depend on $j$ , then for any $j\notin A$ ,

	$\displaystyle\mathbb{E}\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle=\sigma^{2}\quad\text{if $\operatorname{pa}(j)\subset A$,}$
	$\displaystyle\mathbb{E}\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle>\sigma^{2}\quad\text{otherwise.}$

Proof of Lemma B.1.

Let $B_{j}=\operatorname{pa}(X_{j})$ , $\overline{B_{j}}\mathrel{\mathop{\mathchar 58\relax}}=B_{j}-A$ , and $\overline{\mathsf{G}}$ be the subgraph of $\mathsf{G}$ formed by removing the nodes in the ancestral set $A$ . Then

	$\displaystyle\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle=\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{A},X_{B_{j}})\,\|\,X_{A}]+\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{A},X_{B_{j}}]\,\|\,X_{A}]$
		$\displaystyle=\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{A},X_{\overline{B_{j}}})\,\|\,X_{A}]+\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{A},X_{\overline{B_{j}}}]\,\|\,X_{A}].$

There are two cases: (i) $\overline{B_{j}}=\emptyset$ , and (ii) $\overline{B_{j}}\neq\emptyset$ . In case (i), it follows that $B_{j}\subset A$ and hence $\operatorname{pa}(j)\subset A$ . Since $X_{j}$ is conditionally independent of its nondescendants (e.g. ancestors) given its parents, it follows that $\operatorname{var}(X_{j}\,|\,X_{A})=\operatorname{var}(X_{j}\,|\,X_{B_{j}})$ and hence

\displaystyle\mathbb{E}\operatorname{var}(X_{j}\,|\,X_{A})=\mathbb{E}\operatorname{var}(X_{j}\,|\,X_{B_{j}})=\sigma^{2}.

In case (ii), it follows that

	$\displaystyle\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle=\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{A},X_{\overline{B_{j}}})\,\|\,X_{A}]+\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{A},X_{\overline{B_{j}}}]\,\|\,X_{A}]$
		$\displaystyle=\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{B_{j}})\,\|\,X_{A}]+\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{B_{j}}]\,\|\,X_{A}],$

where again we used that $X_{j}$ is conditionally independent of its nondescendants (e.g. ancestors) given its parents to replace conditioning on $(X_{A},X_{\overline{B_{j}}})=X_{A\cup B_{j}}$ with conditioning on $B_{j}$ .

Now suppose $X_{k}$ is in case (i) and $X_{j}$ is in case (ii). We wish to show that $\mathbb{E}\operatorname{var}(X_{j}\,|\,X_{A})>\mathbb{E}\operatorname{var}(X_{k}\,|\,X_{A})=\sigma^{2}$ . Then

	$\displaystyle\mathbb{E}\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle=\mathbb{E}\big{[}\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{B_{j}})\,\|\,X_{A}]\big{]}+\mathbb{E}\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{B_{j}}]\,\|\,X_{A}]$
		$\displaystyle>\mathbb{E}\big{[}\mathbb{E}[\operatorname{var}(X_{j}\,\|\,\operatorname{pa}(j))\,\|\,X_{A}]\big{]}$
		$\displaystyle=\mathbb{E}\operatorname{var}(X_{j}\,\|\,\operatorname{pa}(j))$
		$\displaystyle=\mathbb{E}\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))$
		$\displaystyle=\sigma^{2},$

where we have invoked the assumption that $\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))$ does not depend on $j$ to conclude $\mathbb{E}[\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))\,|\,X_{A}]=\mathbb{E}[\operatorname{var}(X_{k}\,|\,\operatorname{pa}(k))\,|\,X_{A}]$ . This completes the proof. ∎

Theorem 3.1 is an immediate corollary of Lemma B.1. For completeness, we include a proof below.

Proof of Theorem 3.1.

Let $S(\mathsf{G})$ denote the set of sources in $\mathsf{G}$ and note that Lemma B.1 implies that if $X_{s}\in S(\mathsf{G})$ , then $\operatorname{var}(X_{s})<\operatorname{var}(X_{j})$ for any $s\neq j$ . Thus, $S(\mathsf{G})$ is identifiable. Let $\mathsf{G}_{1}$ denote the subgraph of $\mathsf{G}$ formed by removing the nodes in $S(\mathsf{G})$ . Since $S(\mathsf{G})=L_{1}=A_{1}$ , $S(\mathsf{G})$ is an ancestral set in $\mathsf{G}$ . After conditioning on $A_{1}$ , we can thus apply Lemma B.1 once again to identify the sources in $\mathsf{G}_{1}$ , i.e. $S(\mathsf{G}_{1})=L_{2}$ . By repeating this procedure, we can recursively identify $L_{1},\ldots,L_{r}$ , and hence any topological sort of $\mathsf{G}$ . ∎

B.1 Generalization to unequal variances

In this appendix, we illustrate how Theorem 3.1 can be extended to the case where residual variances are different, i.e. $\sigma_{j}^{2}=\mathbb{E}\operatorname{var}(X_{j}|\operatorname{pa}(j))$ is not independent of $j$ . Let $\operatorname{de}(i)$ be the descendant of node $i$ and $[a\mathrel{\mathop{\mathchar 58\relax}}b]=\{a,a+1,\ldots,b-1,b\}$ . Note also that for any nodes $X_{u}$ and $X_{v}$ in the same layer $L_{m}$ of the graph, if we interchange the position of $u$ and $v$ in some true ordering $\pi$ consistent with the graph to get $\pi_{u}$ and $\pi_{v}$ , both $\pi_{u}$ and $\pi_{v}$ are correct orderings.

The following result is similar to existing results on unequal variances (Ghoshal and Honorio, 2018; Park and Kim, 2020; Park, 2020), with the exception that it applies to general DAGs without linearity, additivity, or independent noise.

Theorem B.2.

Suppose there exists an ordering $\pi$ such that for all $j\in[1\mathrel{\mathop{\mathchar 58\relax}}d]$ and $k\in\pi_{[j+1\mathrel{\mathop{\mathchar 58\relax}}d]}$ , the following conditions holds:

If $i=\pi_{j}$ and $k$ are not in the same layer $L_{m}$ , then

\displaystyle\sigma_{i}^{2}

\displaystyle<\sigma_{k}^{2}+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,|\,\operatorname{pa}(k))|X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}}).

(6)

2.

If $i$ and $k$ are in the same layer $L_{m}$ , then either $\sigma^{2}_{i}=\sigma^{2}_{k}$ or (6) holds.

Then the order $\pi$ is identifiable.

In this condition not only do we need to control the descendants of a node, but also the other non-descendants that have not been identified.

Before proving this result, we illustrate it with an example.

Example 6.

Let’s consider a very simple case: A Markov chain with three nodes $X_{1}\rightarrow X_{2}\rightarrow X_{3}$ such that

\displaystyle\begin{split}&X_{1}=z_{1}\sim N(0,1)\\ &X_{2}=\frac{1}{2}X_{1}^{2}+z_{2},\ \ \ \ z_{2}\sim N(0,\tfrac{2}{3})\\ &X_{3}=\frac{1}{3}X_{2}^{2}+z_{3},\ \ \ \ z_{3}\sim N(0,\tfrac{1}{2}).\end{split}

Here we have unequal residual variances. We now check that this model satisfies the conditions in Theorem B.2. Let $f(u)=u^{2}$ and note that the true ordering is $X_{1}\prec X_{2}\prec X_{3}$ . Starting from the first source node $X_{1}$ , we have

\displaystyle\begin{split}\sigma^{2}_{1}&=1\\ \sigma_{2}^{2}+\mathbb{E}\operatorname{var}(f(X_{\pi(2)}))&=2/3+\operatorname{var}(X_{1}^{2}/2)=2/3+1/2>1=\sigma^{2}_{1}\\ \sigma^{2}_{3}+\mathbb{E}\operatorname{var}(f(X_{\pi(3)}))&=1/2+\operatorname{var}(X_{2}^{2}/3)\\ &=1/2+\frac{1}{9}\Big{(}\operatorname{var}(X_{1}^{4})/16+\operatorname{var}(z_{2}^{2})+\operatorname{var}(X_{1}^{2}z_{2})\\ &+2\operatorname{cov}(X_{1}^{4}/4,X_{1}^{2}z_{2})+2\operatorname{cov}(z_{2}^{2},X_{1}^{2}z_{2})\Big{)}\\ &=1/2+8/9+8/81>1=\sigma^{2}_{1}\end{split}

Then for the second source node $X_{2}$ ,

\displaystyle\begin{split}\sigma^{2}_{2}&=2/3\\ \sigma^{2}_{3}+\mathbb{E}\operatorname{var}(f(X_{\pi(3)})\,|\,X_{1})&=1/2+\frac{1}{9}(\operatorname{var}(z^{2}_{2})+\mathbb{E}\operatorname{var}(X^{2}_{1}z_{2}\,|\,X_{1}))\\ &=1/2+1/3-1/81>2/3\end{split}

Thus the condition is satisfied.

If instead we have $\sigma^{2}_{3}=\operatorname{var}(z_{3})=1/3$ , the condition would be violated. It is easy to check that nothing changes for $X_{1}$ . For the second source node $X_{2}$ , things are different:

\displaystyle\sigma^{2}_{3}+\mathbb{E}\operatorname{var}(f(X_{\pi(3)})\,|\,X_{1})=1/3+1/3-1/81<2/3

Thus, the order of $X_{2}$ and $X_{3}$ would be flipped for this model.

This example is easily confirmed in practice. We let $n$ range from 50 to 1000, and check if the estimated order is correct for the two models ( $\sigma_{3}^{2}=1/2$ and $\sigma_{3}^{2}=1/3$ ). We simulated this 50 times and report the averages in Figure 3.

Proof of Theorem B.2.

We first consider the case where every node in the same layer has a different residual variance, i.e. $X_{u},X_{v}\in L_{m}\implies\sigma^{2}_{u}\neq\sigma^{2}_{v}$ . We proceed with induction on the element $j$ of the ordering $\pi$ . Let $i=\pi_{j}$ ,

When $j=1$ and $i=\pi_{1}$ , $X_{i}$ must be a source node and we have for all $k\in\pi_{[2\mathrel{\mathop{\mathchar 58\relax}}d]}$ ,

	$\displaystyle\operatorname{var}(X_{i})=\sigma_{i}^{2}<$	$\displaystyle\sigma_{k}^{2}+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k)))$
	$\displaystyle=$	$\displaystyle\mathbb{E}\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))+\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k)))$
	$\displaystyle=$	$\displaystyle\operatorname{var}(X_{k}).$

Thus the first node to be identified must be $i=\pi_{1}$ , as desired.

Now suppose the the first $j-1$ nodes in the ordering $\pi$ are correctly identified. The parent of node $i=\pi_{j}$ must have been identified in $\pi_{[j-1]}$ or it is a source node. Then we have for all $k\in\{\pi_{j+1},\ldots,\pi_{d}\}$ ,

	$\displaystyle\mathbb{E}\operatorname{var}(X_{i}\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})=\sigma_{i}^{2}<$	$\displaystyle\sigma_{k}^{2}+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$
	$\displaystyle=$	$\displaystyle\mathbb{E}\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$
	$\displaystyle=$	$\displaystyle\mathbb{E}[\mathbb{E}(\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})]+\mathbb{E}[\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\operatorname{var}(X_{k}\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$

Then the $j^{\text{th}}$ node to be identified must be $i=\pi_{j}$ . The induction is completed.

Finally, if $X_{u},X_{v}\in L_{m}$ are in the same layer and $\sigma_{u}^{2}=\sigma_{v}^{2}$ , then this procedure may choose either $X_{u}$ or $X_{v}$ first. For example, if $X_{u}$ is chosen first, then $X_{v}$ will be incorporated into the same layer as $X_{u}$ . Since both these nodes are in the same layer, swapping $X_{u}$ and $X_{v}$ in any ordering still produces a valid ordering of the DAG. The proof is complete. ∎

Appendix C Proof of Theorem 4.1

The proof of Theorem 4.1 will be broken down into several steps. First, we derive an upper bound on the error of the plug-in estimator used in Algorithm 2 (Appendix C.1), and then we derive a uniform upper bound (Appendix C.2). Based on this upper bound, we prove Theorem 4.1 via Proposition C.4 (Appendix C.3). Appendix C.4 collects various technical lemmas that are used throughout.

C.1 A plug-in estimate

Let $(X,Y)\in\mathbb{R}^{m}\times\mathbb{R}$ be a pair of random variables and $\widehat{f}$ be a data-dependent estimator of the conditional expectation $\mathbb{E}[Y\,|\,X]\mathrel{\mathop{\mathchar 58\relax}}=f(X)$ . Assume we have split the sample into two groups, which we denote by $(U^{(1)},V^{(1)}),\ldots,(U^{(n_{1})},V^{(n_{1})})\sim\mathbb{P}(X,Y)$ and $(X^{(1)},Y^{(1)}),\ldots,(X^{(n_{2})},Y^{(n_{2})})\sim\mathbb{P}(X,Y)$ for clarity. Given these samples, define an estimator of $\sigma^{2}_{\operatorname{RV}}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}\operatorname{var}(Y\,|\,X)$ by

\displaystyle\widehat{\sigma}_{\operatorname{RV}}^{2}\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}(Y^{(i)})^{2}-\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\widehat{f}(X^{(i)})^{2}.

(7)

Note here that $\widehat{f}$ depends on $(U^{(i)},V^{(i)})$ , and is independent of the second sample $(X^{(i)},Y^{(i)})$ . We wish to bound the deviation $\mathbb{P}(|\widehat{\sigma}^{2}_{\operatorname{RV}}-\sigma^{2}_{\operatorname{RV}}|\geq t)$ .

Define the target $\theta^{*}=\mathbb{E}f^{2}(X)$ and its plug-in estimator

\displaystyle\theta(g;q)=\int g(x)^{2}\operatorname{d\!}q(x).

Letting $\mathbb{P}_{X}$ denote the true marginal distribution with respect to $x$ , we have $\theta^{*}=\theta(f;\mathbb{P}_{X})$ and $\widehat{\theta}\mathrel{\mathop{\mathchar 58\relax}}=\theta(\widehat{f};\widehat{\mathbb{P}})$ , where $\widehat{\mathbb{P}}=n_{2}^{-1}\sum_{i}\delta_{X^{(i)}}$ is the empirical distribution. We will also make use of more general targets $\theta(g;\mathbb{P}_{X})=\mathbb{E}g^{2}(X)$ for general functions $g$ .

Finally, as a matter of notation, we adopt the following convention: For a random variable $Z$ , $\|Z\|_{p}\mathrel{\mathop{\mathchar 58\relax}}=(\mathbb{E}_{Z}|Z|^{p})^{1/p}$ is the usual $L^{p}$ -norm of $Z$ as a random variable, and for a (nonrandom) function $f$ , $\|f\|_{p}\mathrel{\mathop{\mathchar 58\relax}}=(\int|f|^{p}\,\operatorname{d\!}\zeta)^{1/p}$ , where $\zeta$ is a fixed base measure such as Lebesgue measure. In particular, $\|f(X)\|_{p}\neq\|f\|_{p}$ . The difference of course lies in which measure integrals are taken with respect to. Moreover, we shall always explicitly specify with respect to which variables probabilities and expectations are taken, e.g. to disambiguate $\mathbb{E}_{\widehat{f},X}=\mathbb{E}_{\widehat{f}}\mathbb{E}_{X}$ , $\mathbb{E}_{\widehat{f}}$ , and $\mathbb{E}_{X}$ .

We first prove the following result:

Proposition C.1.

Assume $\|f(X)\|_{\infty},\|\widehat{f}(X)\|_{\infty}\leq B_{\infty}$ . Then

\displaystyle\mathbb{P}(|\widehat{\theta}-\theta^{*}|\geq 2t)\lesssim\frac{\mathbb{E}_{\widehat{f}}\|f(X)-\widehat{f}(X)\|_{2}^{2}}{t^{2}}+\frac{\|f(X)\|_{4}^{4}+\mathbb{E}_{\widehat{f}}\|\widehat{f}(X)-f(X)\|_{4}}{n_{2}t^{2}}.

(8)

Proof.

We have

\displaystyle\mathbb{P}_{\widehat{f},X}(|\widehat{\theta}-\theta^{*}|\geq 2t)\leq\mathbb{P}_{\widehat{f},X}(|\theta^{*}-\theta(\widehat{f};\mathbb{P}_{X})|\geq t)+\mathbb{P}_{\widehat{f},X}(|\widehat{\theta}-\theta(\widehat{f};\mathbb{P}_{X})|\geq t).

The second term is easy to dispense with since

\displaystyle\mathbb{P}_{\widehat{f},X}(|\widehat{\theta}-\theta(\widehat{f};\mathbb{P}_{X})|\geq t)

\displaystyle\leq\frac{\mathbb{E}_{\widehat{f}}\mathbb{E}_{X}(\widehat{\theta}-\theta(\widehat{f};\mathbb{P}_{X}))^{2}}{t^{2}}\leq\frac{\mathbb{E}_{\widehat{f}}\operatorname{var}_{X}(\widehat{f}^{2}(X))}{n_{2}t^{2}}.

(9)

It follows from Lemma C.5 with $p=4$ that

	$\displaystyle\|\mathbb{E}_{X}\widehat{f}(X)^{4}-\mathbb{E}_{X}f(X)^{4}\|$	$\displaystyle\leq\\|\widehat{f}(X)-f(X)\\|_{4}\sum_{k=0}^{3}\\|\widehat{f}(X)\\|_{4}^{k}\\|f(X)\\|_{4}^{3-k}$
		$\displaystyle\leq 4(2B_{\infty})^{3}\\|\widehat{f}(X)-f(X)\\|_{4}$

and hence

\displaystyle\operatorname{var}_{X}(\widehat{f}^{2}(X))\leq\mathbb{E}_{X}\widehat{f}^{4}(X)\lesssim\|f(X)\|_{4}^{4}+\|\widehat{f}(X)-f(X)\|_{4}.

(10)

Combined with (9), we finally have

\displaystyle\mathbb{P}(|\widehat{\theta}-\theta(\widehat{f};\mathbb{P}_{X})|\geq t)\lesssim\frac{\|f(X)\|_{4}^{4}+\mathbb{E}_{\widehat{f}}\|\widehat{f}(X)-f(X)\|_{4}}{n_{2}t^{2}}.

(11)

For the first term, since $f(X),\widehat{f}(X)\in L^{\infty}$ , Lemma C.7 implies

	$\displaystyle\mathbb{E}_{\widehat{f}}(\theta^{*}-\theta(\widehat{f};\mathbb{P}_{X}))^{2}$	$\displaystyle=\mathbb{E}_{\widehat{f}}\mathbb{E}_{X}(f^{2}(X)-\widehat{f}^{2}(X))^{2}$
		$\displaystyle\lesssim\mathbb{E}_{\widehat{f}}\\|f(X)-\widehat{f}(X)\\|_{2}^{2},$

and thus

\displaystyle\mathbb{P}_{\widehat{f}}(|\theta^{*}-\theta(\widehat{f};\mathbb{P}_{X})|\geq t)

\displaystyle\leq\frac{\mathbb{E}_{\widehat{f}}(\theta^{*}-\theta(\widehat{f};\mathbb{P}_{X}))^{2}}{t^{2}}\lesssim\frac{\mathbb{E}_{\widehat{f}}\|f(X)-\widehat{f}(X)\|_{2}^{2}}{t^{2}}.

Therefore

	$\displaystyle\mathbb{P}_{\widehat{f},X}(\|\widehat{\theta}-\theta^{*}\|\geq 2t)$	$\displaystyle\leq\mathbb{P}_{\widehat{f},X}(\|\theta^{*}-\theta(\widehat{f};\mathbb{P}_{X})\|\geq t)+\mathbb{P}_{\widehat{f},X}(\|\widehat{\theta}-\theta(\widehat{f};\mathbb{P}_{X})\|\geq t)$
		$\displaystyle\lesssim\frac{\mathbb{E}_{\widehat{f}}\\|f(X)-\widehat{f}(X)\\|_{2}^{2}}{t^{2}}+\frac{\\|f(X)\\|_{4}^{4}+\mathbb{E}_{\widehat{f}}\\|\widehat{f}(X)-f(X)\\|_{4}}{n_{2}t^{2}}.\qed$

Finally, we conclude the following:

Corollary C.2.

If $\|f(X)\|_{\infty},\|\widehat{f}(X)\|_{\infty}\leq B_{\infty}$ , then

\displaystyle\mathbb{P}(|\widehat{\sigma}_{\operatorname{RV}}^{2}-\sigma_{\operatorname{RV}}^{2}|\geq t)

\displaystyle\lesssim\frac{4}{t^{2}}\Big{(}\mathbb{E}_{\widehat{f}}\|f(X)-\widehat{f}(X)\|_{2}^{2}+\frac{\operatorname{var}(Y)+\|f(X)\|_{4}^{4}+\mathbb{E}_{\widehat{f}}\|\widehat{f}(X)-f(X)\|_{4}}{n_{2}}\Big{)}.

C.2 A uniform bound

For any $j=1,\ldots,r$ and $\ell\notin A_{j}$ , define $\sigma_{\ell j}^{2}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})$ and $\widehat{\sigma}_{\ell j}^{2}$ the corresponding plug-in estimator from (7). By Proposition C.2, we have for $f_{\ell j}(X_{A_{j}})\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{E}[X_{\ell}\,|\,X_{A_{j}}]$ ,

\displaystyle\begin{aligned} \mathbb{P}(|\widehat{\sigma}_{\ell j}^{2}-\sigma_{\ell j}^{2}|\geq t)\lesssim\frac{4}{t^{2}}\Big{(}\mathbb{E}_{\widehat{f}}\|f_{\ell j}(X_{A_{j}})&-\widehat{f}_{\ell j}(X_{A_{j}})\|_{2}^{2}\\ &+\frac{\operatorname{var}(X_{\ell})+\|f_{\ell j}(X_{A_{j}})\|_{4}^{4}+\mathbb{E}_{\widehat{f}}\|\widehat{f}_{\ell j}(X_{A_{j}})-f_{\ell j}(X_{A_{j}})\|_{4}}{n_{2}}\Big{)}\end{aligned}

(12)

Thus we have the following result:

Proposition C.3.

Assume for all $j$ and $\ell\notin A_{j}$ :

1.

$\|f_{\ell j}(X_{A_{j}})\|_{\infty},\|\widehat{f}_{\ell j}(X_{A_{j}})\|_{\infty}\leq B_{\infty}$ ;
2.

$\|f_{\ell j}(X_{A_{j}})\|_{4}^{4}+\operatorname{var}(X_{\ell})\leq B_{\infty}$ ;
3.

$\mathbb{E}_{\widehat{f}}\|f_{\ell j}(X_{A_{j}})-\widehat{f}_{\ell j}(X_{A_{j}})\|_{2}^{2}\to 0$ .

Then

\displaystyle\sup_{\ell,j}\mathbb{P}(|\widehat{\sigma}_{\ell j}^{2}-\sigma_{\ell j}^{2}|>t)\lesssim\frac{4}{t^{2}}\Big{(}\mathbb{E}_{\widehat{f}}\|f_{\ell j}(X_{A_{j}})-\widehat{f}_{\ell j}(X_{A_{j}})\|_{2}^{2}+\frac{1}{n_{2}}\Big{)}.

(13)

For example, under Conditions 1-3, we have $B_{\infty}\mathrel{\mathop{\mathchar 58\relax}}=\sup_{\ell,j}2\|f_{\ell j}(X_{A_{j}})\|_{\infty}^{4}+\zeta_{0}$ .

C.3 Proof of Theorem 4.1

Recall $\sigma_{\ell j}^{2}=\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})$ and $\widehat{\sigma}_{\ell j}^{2}$ is the plug-in estimator defined by (7). Let $\xi>0$ be such that

\displaystyle\sup_{\ell,j}\mathbb{P}(|\widehat{\sigma}_{\ell j}^{2}-\sigma_{\ell j}^{2}|>t)\leq\frac{\xi^{2}}{t^{2}}.

(14)

For example, Proposition C.3 implies $\xi^{2}\asymp\delta^{2}+n_{2}^{-1}$ . Recall also $\Delta\mathrel{\mathop{\mathchar 58\relax}}=\inf_{j}\Delta_{j}$ , where $\Delta_{j}>0$ is the smallest number such that $\mathbb{E}\operatorname{var}(X_{\ell}\,|\,A_{j})>\sigma^{2}+\Delta_{j}$ for all $\ell\notin A_{j}$ .

Theorem 4.1 follows immediately from Proposition C.4 below, combined with Proposition C.3 to bound $\xi$ by $\delta$ .

Proposition C.4.

Define $\xi>0$ as in (14). Then for any threshold $\xi\sqrt{d}<\eta<\Delta/2$ , we have

\displaystyle\mathbb{P}(\widehat{L}=L(\mathsf{G}))

\displaystyle\geq 1-\frac{\xi^{2}}{\eta^{2}}rd.

Proof.

Define $\mathcal{E}_{j-1}\mathrel{\mathop{\mathchar 58\relax}}=\{\widehat{L}_{1}=L_{1},\ldots,\widehat{L}_{j-1}=L_{j-1}\}$ . It follows that

\displaystyle\mathbb{P}(\widehat{L}=L(\mathsf{G}))

\displaystyle=\mathbb{P}(\widehat{L}_{1}=L_{1},\ldots,\widehat{L}_{r}=L_{r})=\prod_{j=1}^{r}\mathbb{P}(\widehat{L}_{j}=L_{j}\,|\,\mathcal{E}_{j-1}).

Clearly, if $\widehat{L}_{1}=L_{1},\ldots,\widehat{L}_{r}=L_{r}$ then $\widehat{r}=r$ . By definition, we have $\sigma_{\ell j}^{2}>\sigma^{2}+\Delta$ , i.e. $\Delta$ is the smallest “gap” between any source in a subgraph $\mathsf{G}[V-A_{j-1}]$ and the rest of the nodes.

Now let $\widehat{\sigma}^{2}\mathrel{\mathop{\mathchar 58\relax}}=\min_{\ell}\widehat{\sigma}_{\ell 0}^{2}$ and consider $L_{1}$ :

\displaystyle\mathbb{P}(\widehat{L}_{1}=L_{1})=\mathbb{P}(|\widehat{\sigma}_{\ell 0}^{2}-\widehat{\sigma}^{2}|\leq\eta\,\forall\ell\in A_{1},\,|\widehat{\sigma}_{\ell 0}^{2}-\widehat{\sigma}^{2}|>\eta\,\forall\ell\notin A_{1})

Now for any $k,\ell\in L_{j}$ ,

\displaystyle|\widehat{\sigma}_{\ell j}^{2}-\widehat{\sigma}_{kj}^{2}|

\displaystyle\leq|\widehat{\sigma}_{\ell j}^{2}-\sigma^{2}|+|\widehat{\sigma}_{kj}^{2}-\sigma^{2}|,

and for any $\ell\notin L_{j}$ and $k\in L_{j}$ ,

\displaystyle|\widehat{\sigma}_{\ell j}^{2}-\widehat{\sigma}_{kj}^{2}|

\displaystyle>\Delta_{j}-|\sigma_{\ell j}^{2}-\widehat{\sigma}_{\ell j}^{2}|-|\widehat{\sigma}_{kj}^{2}-\sigma_{kj}^{2}|.

Thus, with probability $1-d\xi^{2}/t^{2}$ , we have

	$\displaystyle\|\widehat{\sigma}_{\ell j}^{2}-\widehat{\sigma}_{kj}^{2}\|$	$\displaystyle\leq 2t\quad\text{if $k,\ell\in L_{j}$, and}$
	$\displaystyle\|\widehat{\sigma}_{\ell j}^{2}-\widehat{\sigma}_{kj}^{2}\|$	$\displaystyle>\Delta_{j}-2t\quad\text{if $\ell\notin L_{j}$ and $k\in L_{j}$}.$

Now, as long as $t<\Delta/4$ , we have $\Delta_{j}-2t>2t$ , which implies that $\eta\mathrel{\mathop{\mathchar 58\relax}}=2t<\Delta/2$ .

Finally, we have

\displaystyle\mathbb{P}(\widehat{L}_{1}\neq L_{1})

\displaystyle\leq\frac{\xi^{2}}{\eta^{2}}d\implies\mathbb{P}(\widehat{L}_{1}=L_{1})\geq 1-\frac{\xi^{2}}{\eta^{2}}d.

Recall that $d_{j}\mathrel{\mathop{\mathchar 58\relax}}=|L_{j}|$ . Then by a similar argument

\displaystyle\mathbb{P}(\widehat{L}_{2}=L_{2}\,|\,\widehat{L}_{1}=L_{1})

\displaystyle\geq 1-\frac{\xi^{2}}{\eta^{2}}(d-d_{1}).

Recalling $\mathcal{E}_{j-1}\mathrel{\mathop{\mathchar 58\relax}}=\{\widehat{L}_{1}=L_{1},\ldots,\widehat{L}_{j-1}=L_{j-1}\}$ , we have just proved that $\mathbb{P}(\widehat{L}_{2}=L_{2}\,|\,\mathcal{E}_{1})\geq 1-(d-d_{1})(\xi^{2}/\eta^{2})$ . A similar argument proves that $\mathbb{P}(\widehat{L}_{j}=L_{j}\,|\,\mathcal{E}_{j-1})\geq 1-(\xi^{2}/\eta^{2})(d-d_{j-1})$ . Since $\eta>\xi\sqrt{d}$ , the inequality $\prod_{j}(1-x_{j})\geq 1-\sum_{j}x_{j}$ implies

	$\displaystyle\mathbb{P}(\widehat{L}=L(\mathsf{G}))$	$\displaystyle=\prod_{j=1}^{r}\mathbb{P}(\widehat{L}_{j}=L_{j}\,\|\,\mathcal{E}_{j-1})$
		$\displaystyle=\prod_{j=1}^{r}\Big{(}1-\frac{\xi^{2}}{\eta^{2}}(d-d_{j-1})\Big{)}$
		$\displaystyle\geq 1-\sum_{j=1}^{r}\frac{\xi^{2}}{\eta^{2}}(d-d_{j-1})$
		$\displaystyle\geq 1-\frac{\xi^{2}}{\eta^{2}}rd$

as desired. ∎

C.4 Technical lemmas

Lemma C.5.

\displaystyle|\mathbb{E}X^{p}-\mathbb{E}Y^{p}|\leq\|X-Y\|_{p}\sum_{k=0}^{p-1}\|X\|_{p}^{k}\|Y\|_{p}^{p-1-k}

Proof.

Write $\mathbb{E}X^{p}-\mathbb{E}Y^{p}$ as a telescoping sum:

	$\displaystyle\|\mathbb{E}X^{p}-\mathbb{E}Y^{p}\|=\|\\|X\\|_{p}^{p}-\\|Y\\|_{p}^{p}\|$	$\displaystyle=\|\\|X\\|_{p}-\\|Y\\|_{p}\|\cdot\sum_{k=0}^{p-1}\\|X\\|_{p}^{k}\\|Y\\|_{p}^{p-k-1}$
		$\displaystyle\leq\\|X-Y\\|_{p}\cdot\sum_{k=0}^{p-1}\\|X\\|_{p}^{k}\\|Y\\|_{p}^{p-k-1}.\qed$

Lemma C.6.

Fix $p>2$ and $\delta>0$ and suppose $\|f\|_{p+\delta},\|g\|_{p+\delta}\leq B_{p+\delta}$ . Then

\displaystyle\|f-g\|_{p}\leq C_{p,\delta}\cdot\|f-g\|_{2}^{\gamma_{p,\delta}},\quad C_{p,\delta}=(2B_{p+\delta})^{\tfrac{(p-2)(p+\delta)}{p(p+\delta-2)}},\quad\gamma_{p,\delta}=\frac{2\delta}{p(p+\delta-2)}.

The exponent $\gamma_{p,\delta}$ satisfies $\gamma_{p,\delta}\leq 2/p<1$ and $\gamma_{p,\delta}\to 2/p$ as $\delta\to\infty$ , and the constant $C_{p,\delta}\to 1-\tfrac{2}{p}$ as $\delta\to\infty$ . Thus, if $f-g\in L^{\infty}$ , then

\displaystyle\|f-g\|_{p}\lesssim\|f-g\|_{2}^{2/p}.

Proof.

Use log-convexity of $L^{p}$ -norms with $2=q<p<r=p+\delta$ . ∎

Lemma C.7.

Assume $\|f\|_{\infty}\leq B_{\infty}<\infty$ and $g\in L^{4}$ . Then

\displaystyle\begin{aligned} \mathbb{E}_{X}(f(X)^{2}-g(X)^{2})^{2}\leq\|f(X)&-g(X))\|_{4}^{4}+\\ &4B_{\infty}\|f(X)-g(X)\|_{3}^{3}+4B_{\infty}^{2}\|f(X)-g(X)\|_{2}^{2}.\end{aligned}

(15)

If additionally $g(X)\in L^{\infty}$ , then

\displaystyle\mathbb{E}_{X}(f(X)^{2}-g(X)^{2})^{2}

\displaystyle\lesssim\|f(X)-g(X)\|_{2}^{2}.

(16)

Proof.

Note that

	$\displaystyle(f(X)^{2}-g(X)^{2})^{2}$	$\displaystyle=(g(X)-f(X))^{4}+4f(X)(g(X)-f(X))^{3}+4f(X)^{2}(g(X)-f(X))^{2}$
		$\displaystyle\leq(g(X)-f(X))^{4}+4\|f(X)\|\|g(X)-f(X)\|^{3}+4f(X)^{2}(g(X)-f(X))^{2}.$

Thus

	$\displaystyle\mathbb{E}(f(X)^{2}-g(X)^{2})^{2}$	$\displaystyle\leq\mathbb{E}(g(X)-f(X))^{4}+$
		$\displaystyle\qquad\qquad 4\mathbb{E}\big{[}\|f(X)\|\|g(X)-f(X)\|^{3}\big{]}+4\mathbb{E}\big{[}f(X)^{2}(g(X)-f(X))^{2}\big{]}$
		$\displaystyle\leq\mathbb{E}(g(X)-f(X))^{4}+$
		$\displaystyle\qquad\qquad 4B_{\infty}\mathbb{E}\|g(X)-f(X)\|^{3}+4B_{\infty}^{2}\mathbb{E}(g(X)-f(X))^{2}.$

This proves the first inequality. The second follows from taking $\delta\to\infty$ in Lemma C.6 and using $\|f-g\|_{p}^{p}\lesssim\|f-g\|_{2}^{2}$ for $p=3,4$ . ∎

Appendix D Comparison to CAM algorithm

In this section, we justify the claim in Example 5 that there exist infinitely many nonlinear functions $g$ for which the CAM algorithm returns an incorrect graph under the model (4). To show this, we first construct a linear model on which CAM returns an incorrect ordering. Since CAM focuses on nonlinear models, we then show that this extends to any sufficiently small nonlinear perturbation of this model.

The linear model is

\displaystyle\left\{\begin{aligned} &X_{1}\sim\mathcal{N}(0,1)\\ &X_{2}=X_{1}+z_{2}\ \ \ \ z_{2}\sim\mathcal{N}(0,1)\\ &X_{3}=X_{1}+X_{2}+z_{3}\ \ \ \ z_{3}\sim\mathcal{N}(0,1).\end{aligned}\right.

(17)

The graph is

\displaystyle\begin{split}X_{1}\rightarrow&X_{2}\\ \searrow&\downarrow\\ &X_{3}\end{split}

which corresponds to the adjacency matrix

\displaystyle\begin{pmatrix}0&1&1\\ 0&0&1\\ 0&0&0\\ \end{pmatrix}.

The IncEdge step of the CAM algorithm (§5.2, (Bühlmann et al., 2014)) is based on the following score function:

\displaystyle\sum_{j=1}^{d}\log\Big{(}\mathbb{E}\operatorname{var}(X_{j}\,|\,\operatorname{pa}(j))\Big{)}.

The algorithm starts with the empty DAG (i.e. $\operatorname{pa}(j)=\emptyset$ for all $j$ ) and proceeds by greedily adding edges that decrease this score the most in each step. For example, in the first step, CAM searches for the pair $(X_{i},X_{j})$ that maximizes $\log\operatorname{var}(X_{j})-\log\mathbb{E}\operatorname{var}(X_{j}\,|\,X_{i})$ , and adds the edge $X_{i}\to X_{j}$ to the estimated DAG. The second proceeds similarly until the estimated order is determined. Thus, it suffices to study the log-differences $\omega(j,i,S)\mathrel{\mathop{\mathchar 58\relax}}=\log\operatorname{var}(X_{j}\,|\,X_{S})-\log\mathbb{E}\operatorname{var}(X_{j}\,|\,X_{S\cup i})$ .

The following are straightforward to compute for the model (17):

\displaystyle\begin{split}&\operatorname{var}(X_{1})=1\ \ \ \ \mathbb{E}\operatorname{var}(X_{2}|X_{1})=1\ \ \ \ \mathbb{E}\operatorname{var}(X_{1}|X_{2})=\frac{1}{2}\\ &\operatorname{var}(X_{2})=2\ \ \ \ \mathbb{E}\operatorname{var}(X_{3}|X_{1})=2\ \ \ \ \mathbb{E}\operatorname{var}(X_{1}|X_{3})=\frac{3}{2}\\ &\operatorname{var}(X_{3})=6\ \ \ \ \mathbb{E}\operatorname{var}(X_{3}|X_{2})=\frac{1}{3}\ \ \ \ \!\mathbb{E}\operatorname{var}(X_{2}|X_{3})=\frac{1}{2}\\ \end{split}

Then

\displaystyle\begin{split}\log\Big{(}\frac{\operatorname{var}(X_{2})}{\mathbb{E}\operatorname{var}(X_{2}|X_{1})}\Big{)}=\log\Big{(}\frac{2}{1}\Big{)}=\log 2\ \ \ \ &\log\Big{(}\frac{\operatorname{var}(X_{1})}{\mathbb{E}\operatorname{var}(X_{1}|X_{2})}\Big{)}=\log\Big{(}\frac{1}{1/2}\Big{)}=\log 2\\ \log\Big{(}\frac{\operatorname{var}(X_{3})}{\mathbb{E}\operatorname{var}(X_{3}|X_{1})}\Big{)}=\log\Big{(}\frac{6}{2}\Big{)}=\log 3\ \ \ \ &\log\Big{(}\frac{\operatorname{var}(X_{1})}{\mathbb{E}\operatorname{var}(X_{1}|X_{3})}\Big{)}=\log\Big{(}\frac{1}{1/3}\Big{)}=\log 3\\ \log\Big{(}\frac{\operatorname{var}(X_{3})}{\mathbb{E}\operatorname{var}(X_{3}|X_{2})}\Big{)}=\log\Big{(}\frac{6}{3/2}\Big{)}=\log 4\ \ \ \ &\log\Big{(}\frac{\operatorname{var}(X_{2})}{\mathbb{E}\operatorname{var}(X_{2}|X_{3})}\Big{)}=\log\Big{(}\frac{2}{1/2}\Big{)}=\log 4\\ \end{split}

Now, if $X_{3}\to X_{2}$ is chosen first, then the order is incorrect and we are done. Thus suppose CAM instead chooses $X_{2}\to X_{3}$ , then in the next step it would update the score for $X_{1}\to X_{3}$ to be

\displaystyle\log\Big{(}\frac{\mathbb{E}\operatorname{var}(X_{3}|X_{2})}{\mathbb{E}\operatorname{var}(X_{3}|X_{1},X_{2})}\Big{)}=\log\Big{(}\frac{3/2}{1}\Big{)}=\log\frac{3}{2}<\log\Big{(}\frac{\mathbb{E}\operatorname{var}(X_{1}|X_{2})}{\mathbb{E}\operatorname{var}(X_{1}|X_{3},X_{2})}\Big{)}=\log 3

Therefore, for the next edge, CAM would choose $X_{3}\to X_{1}$ , which also leads to the wrong order. Thus regardless of which edge is selected first, CAM will return the wrong order.

Thus, when CAM is applied to data generated from (17), it is guaranteed to return an incorrect ordering. Although the model (17) is identifiable, it does not satisfy the identifiability condition for CAM (Lemma 1, (Bühlmann et al., 2014)), namely that the structural equation model is a nonlinear additive model. Thus, we need to extend this example to an identifiable, nonlinear additive model.

Since this depends only on the scores $\omega(j,i,S)$ , it suffices to construct a nonlinear model with similar scores. For this, we consider a simple nonlinear extension of (17): Let $g$ be an arbitrary bounded, nonlinear function, and define $g_{\delta}(u)\mathrel{\mathop{\mathchar 58\relax}}=u+\delta g(u)$ . The nonlinear model is given by

\displaystyle\left\{\begin{aligned} &X_{1}\sim\mathcal{N}(0,1)\\ &X_{2}=g_{\delta}(X_{1})+z_{2}\ \ \ \ z_{2}\sim\mathcal{N}(0,1)\\ &X_{3}=g_{\delta}(X_{1})+g_{\delta}(X_{2})+z_{3}\ \ \ \ z_{3}\sim\mathcal{N}(0,1).\end{aligned}\right.

(18)

This model satisfies both our identifiability condition (Condition 2) and the identifiability condition for CAM (Lemma 1, (Bühlmann et al., 2014)).

We claim that for sufficiently small $\delta$ , the CAM algorithm will return the wrong ordering (see Proposition D.1 below for a formal statement). It follows that the scores $\omega(j,i,S;\delta)$ corresponding to the model (18) can be made arbitrarily close to $\omega(j,i,S)=\omega(j,i,S;0)$ , which implies that CAM will return the wrong ordering for sufficiently small $\delta>0$ .

In Figure 4, we illustrate this empirically. In addition to the model (18), we also simulated from the following model, which shows that this phenomenon is not peculiar to the construction above:

\displaystyle\left\{\begin{aligned} &X_{1}\sim\mathcal{N}(0,1)\\ &X_{2}=X_{1}^{2}+z_{2}\ \ \ \ z_{2}\sim\mathcal{N}(0,1)\\ &X_{3}=4X_{1}^{2}+h(X_{2})+z_{3}\ \ \ \ z_{3}\sim\mathcal{N}(0,1).\end{aligned}\right.

(19)

In all eight examples, NPVAR perfectly recovers the ordering while CAM is guaranteed to return an inconsistent order for sufficiently large $n$ (i.e. once the scores are consistently estimated).

Proposition D.1.

Let $\mathbb{E}_{\delta}$ and $\operatorname{var}_{\delta}$ be taken with respect to model (18). Then for all $i,j\in\{1,2,3\}$ , as $\delta\to 0$ ,

\displaystyle\begin{split}|\operatorname{var}_{\delta}(X_{i})-\operatorname{var}_{0}(X_{i})|&=o(1),\\ |\mathbb{E}_{\delta}\operatorname{var}_{\delta}(X_{i}|X_{j})-\mathbb{E}_{0}\operatorname{var}_{0}(X_{i}|X_{j})|&=o(1).\end{split}

and $|\mathbb{E}_{\delta}\operatorname{var}_{\delta}(X_{3}|X_{1},X_{2})-\mathbb{E}_{0}\operatorname{var}_{0}(X_{3}|X_{1},X_{2})|=o(1)$ .

Proof sketch of Proposition D.1.

The proof is consequence of the fact that the differences are continuous functions of $\delta$ . We sketch the proof for $i=2$ ; the remaining cases are similar.

We have

\begin{split}&\operatorname{var}_{\delta}(X_{2})=\operatorname{var}_{\delta}(g(X_{1}))+\operatorname{var}_{\delta}(\epsilon)\\ &\operatorname{var}_{0}(X_{2})=\operatorname{var}_{0}(X_{1})+\operatorname{var}_{0}(\epsilon)\\ \end{split}

Let $\varphi(t)$ be the standard normal density. We only need to analyze and compare $\operatorname{var}_{0}(X_{1})-\operatorname{var}_{\delta}(g(X_{1}))=\operatorname{var}_{0}(X_{1})-\operatorname{var}_{0}(g(X_{1}))$ in two parts:

\begin{split}&\int(X_{1}^{2}-g(X_{1})^{2})\varphi(X_{1})dX_{1}\\ &(\int X_{1}\varphi(X_{1})dX_{1})^{2}-(\int g(X_{1})\varphi(X_{1})dX_{1})^{2}.\end{split}

Since $|X_{1}-g(X_{1})|\leq\delta$ ,

\begin{split}&|\int(X_{1}^{2}-g(X_{1})^{2})\varphi(X_{1})dX_{1}|\leq\delta\int|X_{1}+g(X_{1})|\varphi(X_{1})dX_{1}\leq\delta\int 2|X_{1}|\varphi(X_{1})dX_{1}+\delta^{2}=C\delta+\delta^{2}\\ &|\int g(X_{1})\varphi(X_{1})dX_{1}|=|\int(X_{1}-g(X_{1}))\varphi(X_{1})dX_{1}|\leq\delta\int\varphi(X_{1})dX_{1}=\delta\\ &|\int(X_{1}+g(X_{1})\varphi(X_{1}))dX_{1}|=|\mathbb{E}g(X_{1})|\leq\delta.\end{split}

Thus

|(\int X_{1}\varphi(X_{1})dX_{1})^{2}-(\int g(X_{1})\varphi(X_{1})dX_{1})^{2}|\leq\delta^{2},

so that

|\operatorname{var}_{\delta}(X_{2})-\operatorname{var}_{0}(X_{2})|=o(1)

as claimed. ∎

Appendix E Experiment details

In this appendix we outline the details of our experiments, as well as additional simulations.

E.1 Experiment settings

For graphs, we used

•

Markov chain (MC). Graph where there is one edge from $X_{i-1}$ to $X_{i}$ for all nodes $i=2,\ldots,d$ .
•

Erdős Rényi (ER). Random graphs whose edges are added independently with specified expected number of edges.
•

Scale-free networks (SF). Networks simulated according to the Barabasi-Albert model.

For models, we refer to the nonlinear functions in SEM. We specify the nonlinear functions in $X_{j}=f_{j}(X_{\text{pa(j)}})+z_{j}$ for all $j=1,2,\ldots,d$ , where $z_{j}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2})$ with variance $\sigma^{2}\in\{0.2,0.5,0.8\}$

•

Additive sine model (SIN): $f_{j}(X_{\text{pa}(j)})=\sum_{k\in\text{pa}(j)}f_{jk}(X_{k})$ where $f_{jk}(X_{k})=\sin(X_{k})$ .
•

Additive Gaussian process (AGP): $f_{j}(X_{\text{pa}(j)})=\sum_{k\in\text{pa}(j)}f_{jk}(X_{k})$ where $f_{jk}$ is a draw from Gaussian process with RBF kernel with length-scale one.
•

Non-additive Gaussian process (NGP): $f_{j}$ is a draw from Gaussian process with RBF kernel with length-scale one.
•

Generalized Linear Model (GLM): This is a special case with non-additive model and non-additive noise. We specify the model $\mathbb{P}(X_{j}=1)=f_{j}(X_{\text{pa}(j)})$ by a parameter $p\in\{0.1,0.3\}$ . Given $p$ , if $j$ is a source node, $\mathbb{P}(X_{j}=1)=p$ . For the Markov chain model we define:

$\displaystyle f_{j}(X_{\text{pa}(j)})=\left\{\begin{aligned} p&\ \ \ \ &X_{k}=1\\ 1-p&\ \ \ \ &X_{k}=0\\ \end{aligned}\right.$

or vice versa.

We generated graphs from each of the above models with $\{d,4d\}$ edges each. These are denoted by the shorthand XX-YYY- $k$ , where XX denotes the graph type, YYY denotes the model, and $k$ indicates the graphs have $kd$ edges on average. For example, ER-SIN-1 indicates an ER graph with $d$ (expected) edges under the additive sine model. SF-NGP-4 indicates an SF graph with $4d$ (expected) edges under a non-additive Gaussian process. Note that there is no difference between additive or non-additive GP for Markov chains, so the results for MC-NGP- $k$ are omitted from the Figures in Appendix E.5.

Based on these models, we generated random datasets with $n$ samples. For each simulation run, we generated $n\in\{100,200,500,750,1000\}$ samples for graphs with $d\in\{5,10,20,40,50,60,70\}$ nodes.

E.2 Implementation and baselines

Code implementing NPVAR can be found at https://github.com/MingGao97/NPVAR. We implemented NPVAR (Algorithm 2) with generalized additive models (GAMs) as the nonparametric estimator for $\widehat{f}_{\ell j}$ . We used the gam function in the R package mgcv with P-splines bs=’ps’ and the default smoothing parameter sp=0.6. In all our implementations, including our own for Algorithm 2, we used default parameters in order to avoid skewing the results in favour of any particular algorithm as a result of hyperparameter tuning.

We compared our method with following approaches as baselines:

•

Regression with subsequent independence test (RESIT) identifies and disregards a sink node at each step via independence testing (Peters et al., 2014). The implementation is available at http://people.tuebingen.mpg.de/jpeters/onlineCodeANM.zip. Uses HSIC test for independence testing with alpha=0.05 and gam for nonparametric regression.
•

Causal additive models (CAM) estimates the topological order by greedy search over edges after a preliminary neighborhood selection (Bühlmann et al., 2014). The implementation is available at https://cran.r-project.org/src/contrib/Archive/CAM/. By default, CAM applies an extra pre-processing step called preliminary neighborhood search (PNS). Uses gam to compute the scores and for pruning, and mboost for preliminary neighborhood search. The R implementation of CAM does not use the default parameters for gam or mboost, and instead optimizes these parameters at runtime.
•

NOTEARS uses an algebraic characterization of DAGs for score-based structure learning of nonparametric models via partial derivatives (zheng2018notears; Zheng et al., 2020). The implementation is available at https://github.com/xunzheng/notears. We used neural networks for the nonlinearities with a single hidden layer of 10 neurons. The training parameters are lambda1=0.01, lambda2=0.01 and the threshold for adjacency matrix is w_threshold=0.3.
•

Greedy equivalence search with generalized scores (GSGES) uses generalized scores for greedy search without assuming model class (Huang et al., 2018). The implementation is available at https://github.com/Biwei-Huang/Generalized-Score-Functions-for-Causal-Discovery/. We used cross-validation parameters parameters.kfold = 10 and parameters.lambda = 0.01.
•

Equal variance (EqVar) algorithm identifies source node by minimizing conditional variance in linear SEM (Chen et al., 2019). The implementation is available at https://github.com/WY-Chen/EqVarDAG. The original EqVar algorithm estimates the error variances in a linear SEM via the covariance matrix $\Sigma=\mathbb{E}XX^{T}$ , and then uses linear regression (e.g. best subset selection) to learn the structure of the DAG. We adapted this algorithm to the nonlinear setting in the obivous way by using GAMs (instead of subset selection) for variable selection. The use of the covariance matrix to estimate the order remains the same.
•

PC algorithm and greedy equivalence search (GES) are standard baselines for structure learning (Spirtes and Glymour, 1991; Chickering, 2003). The implementation is available from R package pcalg. For PC algorithm, use correlation matrix as sufficient statistic. Independence test is implemented by gaussCItest with significance level alpha=0.01. For GES, set the score to be GaussL0penObsScore with lambda as $\log n/2$ , which corresponds to the BIC score.

The experiments were conducted on an Intel E5-2680v4 2.4GHz CPU with 64 GB memory.

E.3 Metrics

We evaluated the performance of each algorithm with the following two metrics:

•

$\mathbb{P}(\text{correct order})$ : The percentage of runs in which the algorithm gives a correct topological ordering, over $N$ runs. This metric is only sensible for algorithms that first estimate an ordering or return an adjacency matrix which does not contain undirected edges, including RESIT, CAM, EqVar and NOTEARS.
•

Structural Hamming distance (SHD): A standard benchmark in the structure learning literature that counts the total number of edge additions, deletions, and reversals needed to convert the estimated graph into the true graph.

Since there may be multiple topological orderings of a DAG, in our evaluations of order recovery, we check whether or not the order returned is any of the possible valid orderings. For PC, GES, and GSGES, they all return a CPDAG that may contain undirected edges, in which case we evaluate them favourably by assuming correct orientation for undirected edges whenever possible. Since CAM, RESIT, EqVar and NPVAR each first estimate a topological ordering then estimate a DAG. To estimate a DAG from an ordering, we apply the same pruning step to each algorithm for a fair comparison, which is adapted from (Bühlmann et al., 2014). Specifically, given an estimated ordering, we run a gam regression for each node on its ancestors, then determine the parents of the node by the $p$ -values with significance level cutoff=0.001 for estimating the DAG.

E.4 Timing

For completeness, runtime comparisons are reported in Tables 1 and 2. Algorithms based on linear models such as EqVar, PC, and GES are by far the fastest. These algorithms are also the most highly optimized. The slowest algorithms are GSGES and RESIT. Timing comparisons against CAM are are difficult to interpret since by default, CAM first performs preliminary neighbourhood search, which can easily be applied to any of the other algorithms tested. The dramatic difference with and without this pre-processing step by comparing Tables 1 and 2: For $d=40$ , with this extra step CAM takes just over 90 seconds, whereas without it, on the same data, it takes over 8.5 hours. For comparison, NPVAR takes around two minutes (i.e. without pre-processing or neighbourhood search).

Algorithm	$d$	$n$	Runtime (s)
EqVar	20	1000	$0.0017\pm 0.0003$
PC	20	1000	$0.056\pm 0.016$
GES	20	1000	$0.060\pm 0.034$
NPVAR	20	1000	$10.76\pm 0.23$
NOTEARS	20	1000	$31.46\pm 8.79$
CAM (w/ PNS)	20	1000	$40.56\pm 1.29$
CAM (w/o PNS)	20	1000	$559.01\pm 9.49$
RESIT	20	1000	$652.15\pm 7.26$
GSGES	20	1000	$3216.00\pm 95.0$

Table 1: Runtime comparisons for

d=20

. Timing for CAM is presented with and without preliminary neighbourhood selection (PNS), which is a pre-processing step that can be applied to any algorithm. In our experiments, only CAM used PNS.

Algorithm	$d$	$n$	Runtime (s)
EqVar	40	1000	$0.0043\pm 0.0003$
GES	40	1000	$0.12\pm 0.0052$
PC	40	1000	$0.019\pm 0.030$
NOTEARS	40	1000	$76.05\pm 19.16$
CAM (w/ PNS)	40	1000	$95.59\pm 6.33$
NPVAR	40	1000	$118.33\pm 2.25$
CAM (w/o PNS)	40	1000	$31644.56\pm 1329.31$

Table 2: Runtime comparisons for

d=40

. Timing for CAM is presented with and without preliminary neighbourhood selection (PNS), which is a pre-processing step that can be applied to any algorithm. In our experiments, only CAM used PNS.

E.5 Additional experiments

Here we collect the results of our additional experiments. Since the settings MC-AGP- $k$ and MC-NGP- $k$ are equivalent (i.e. since there is only parent for each node), the plots for MC-NGP- $k$ are omitted. Some algorithms might be skipped due to high computational cost or numerical issue.

•

Figure 5: SHD vs. $n$ with $d=5$ fixed, across all graphs and models tested.
•

Figure 6: SHD vs. $n$ with $d=10$ fixed, across all graphs and models tested.
•

Figure 7: SHD vs. $n$ with $d=20$ fixed, across all graphs and models tested.
•

Figure 8: SHD vs. $n$ with $d=40$ fixed, across all graphs and models tested with GSGES and RESIT skipped (due to high computational cost).
•

Figure 9: SHD vs. $n$ with $d=50$ fixed, across all graphs and models tested with GSGES, RESIT and NOTEARS skipped (due to high computational cost).
•

Figure 10: SHD vs. $n$ with $d=60$ fixed, across all graphs and models tested with GSGES, RESIT and NOTEARS skipped (due to high computational cost).
•

Figure 11: SHD vs. $n$ with $d=70$ fixed, across all graphs and models tested with GSGES, RESIT and NOTEARS skipped (due to high computational cost).
•

Figure 12: SHD vs. $n$ with $d$ ranging from $5$ to $70$ on GLM with CAM and RESIT skipped (due to numerical issues).
•

Figure 13: SHD vs. $d$ with $n=500$ fixed, across all graphs and models tested.
•

Figure 14: SHD vs. $d$ with $n=1000$ fixed, across all graphs and models tested.
•

Figure 15: Ordering recovery vs. $n$ with $d=5$ fixed, across all graphs and models tested.
•

Figure 16: Ordering recovery vs. $n$ with $d=10$ fixed, across all graphs and models tested.
•

Figure 17: Ordering recovery vs. $n$ with $d=20$ fixed, across all graphs and models tested.
•

Figure 18: Ordering recovery vs. $n$ with $d=40$ fixed, across all graphs and models tested with RESIT skipped (due to high computational cost).
•

Figure 19: Ordering recovery vs. $n$ with $d=50$ fixed, across all graphs and models tested with RESIT and NOTEARS skipped (due to high computational cost).
•

Figure 20: Ordering recovery vs. $n$ with $d=60$ fixed, across all graphs and models tested with RESIT and NOTEARS skipped (due to high computational cost).
•

Figure 21: Ordering recovery vs. $n$ with $d=70$ fixed, across all graphs and models tested with RESIT and NOTEARS skipped (due to high computational cost).

References

Aragam and Zhou (2015) B. Aragam and Q. Zhou. Concave Penalized Estimation of Sparse Gaussian Bayesian Networks. Journal of Machine Learning Research, 16(69):2273–2328, 2015.
Aragam et al. (2015) B. Aragam, A. A. Amini, and Q. Zhou. Learning directed acyclic graphs with penalized neighbourhood regression. arXiv:1511.08963, 2015.
Aragam et al. (2019) B. Aragam, A. Amini, and Q. Zhou. Globally optimal score-based learning of directed acyclic graphs in high-dimensions. In Advances in Neural Information Processing Systems 32, pages 4450–4462. 2019.
Bertin and Lecué (2008) K. Bertin and G. Lecué. Selection of variables and dimension reduction in high-dimensional non-parametric regression. Electronic Journal of Statistics, 2:1224–1241, 2008.
Bottou (2015) L. Bottou. Two big challenges in machine learning. ICML 2015 keynote lecture, 2015.
Bühlmann (2018) P. Bühlmann. Invariance, causality and robustness. arXiv preprint arXiv:1812.08233, 2018.
Bühlmann et al. (2014) P. Bühlmann, J. Peters, and J. Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. Annals of Statistics, 42(6):2526–2556, 2014.
Chen et al. (2019) W. Chen, M. Drton, and Y. S. Wang. On causal discovery with an equal-variance assumption. Biometrika, 106(4):973–980, 09 2019. ISSN 0006-3444. doi: 10.1093/biomet/asz049.
Chicharro et al. (2019) D. Chicharro, S. Panzeri, and I. Shpitser. Conditionally-additive-noise models for structure learning. arXiv preprint arXiv:1905.08360, 2019.
Chickering (2003) D. M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507–554, 2003.
Comminges and Dalalyan (2011) L. Comminges and A. S. Dalalyan. Tight conditions for consistent variable selection in high dimensional nonparametric regression. In Proceedings of the 24th Annual Conference on Learning Theory, pages 187–206, 2011.
Dawid (1980) A. P. Dawid. Conditional independence for statistical operations. Annals of Statistics, pages 598–617, 1980.
Doksum et al. (1995) K. Doksum, A. Samarov, et al. Nonparametric estimation of global functionals and a measure of the explanatory power of covariates in regression. The Annals of Statistics, 23(5):1443–1473, 1995.
Ghoshal and Honorio (2017) A. Ghoshal and J. Honorio. Learning identifiable gaussian bayesian networks in polynomial time and sample complexity. In Advances in Neural Information Processing Systems 30, pages 6457–6466. 2017.
Ghoshal and Honorio (2018) A. Ghoshal and J. Honorio. Learning linear structural equation models in polynomial time and sample complexity. In A. Storkey and F. Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1466–1475, Playa Blanca, Lanzarote, Canary Islands, 09–11 Apr 2018. PMLR.
Giordano et al. (2020) F. Giordano, S. N. Lahiri, and M. L. Parrella. GRID: A variable selection and structure discovery method for high dimensional nonparametric regression. Annals of Statistics, to appear, 2020.
Györfi et al. (2006) L. Györfi, M. Kohler, A. Krzyzak, and H. Walk. A distribution-free theory of nonparametric regression. Springer Science & Business Media, 2006.
Hastie and Tibshirani (1990) T. J. Hastie and R. J. Tibshirani. Generalized additive models, volume 43. CRC press, 1990.
Hoyer et al. (2009) P. O. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in neural information processing systems, pages 689–696, 2009.
Huang et al. (2018) B. Huang, K. Zhang, Y. Lin, B. Schölkopf, and C. Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1551–1560. ACM, 2018.
Kandasamy et al. (2015) K. Kandasamy, A. Krishnamurthy, B. Poczos, L. Wasserman, et al. Nonparametric von mises estimators for entropies, divergences and mutual informations. In Advances in Neural Information Processing Systems, pages 397–405, 2015.
Kohler et al. (2009) M. Kohler, A. Krzyżak, and H. Walk. Optimal global rates of convergence for nonparametric regression with unbounded data. Journal of Statistical Planning and Inference, 139(4):1286–1296, 2009.
Koller and Friedman (2009) D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
Lachapelle et al. (2019) S. Lachapelle, P. Brouillard, T. Deleu, and S. Lacoste-Julien. Gradient-based neural dag learning. arXiv preprint arXiv:1906.02226, 2019.
Lafferty and Wasserman (2008) J. Lafferty and L. Wasserman. Rodeo: sparse, greedy nonparametric regression. Annals of Statistics, 36(1):28–63, 2008.
Lauritzen (1996) S. L. Lauritzen. Graphical models. Oxford University Press, 1996.
Loh and Bühlmann (2014) P.-L. Loh and P. Bühlmann. High-dimensional learning of linear causal networks via inverse covariance estimation. Journal of Machine Learning Research, 15:3065–3105, 2014.
Miller and Hall (2010) H. Miller and P. Hall. Local polynomial regression and variable selection. In Borrowing Strength: Theory Powering Applications–A Festschrift for Lawrence D. Brown, pages 216–233. Institute of Mathematical Statistics, 2010.
Mitrovic et al. (2018) J. Mitrovic, D. Sejdinovic, and Y. W. Teh. Causal inference via kernel deviance measures. In Advances in Neural Information Processing Systems, pages 6986–6994, 2018.
Monti et al. (2019) R. P. Monti, K. Zhang, and A. Hyvarinen. Causal discovery with general non-linear relationships using non-linear ica. arXiv preprint arXiv:1904.09096, 2019.
Mooij et al. (2016) J. M. Mooij, J. Peters, D. Janzing, J. Zscheischler, and B. Schölkopf. Distinguishing cause from effect using observational data: methods and benchmarks. The Journal of Machine Learning Research, 17(1):1103–1204, 2016.
Ng et al. (2019) I. Ng, Z. Fang, S. Zhu, and Z. Chen. Masked gradient-based causal structure learning. arXiv preprint arXiv:1910.08527, 2019.
Nowzohour and Bühlmann (2016) C. Nowzohour and P. Bühlmann. Score-based causal learning in additive noise models. Statistics, 50(3):471–485, 2016.
Ott et al. (2004) S. Ott, S. Imoto, and S. Miyano. Finding optimal models for small gene networks. In Pacific symposium on biocomputing, volume 9, pages 557–567. Citeseer, 2004.
Park (2018) G. Park. Learning generalized hypergeometric distribution (ghd) dag models. arXiv preprint arXiv:1805.02848, 2018.
Park (2020) G. Park. Identifiability of additive noise models using conditional variances. Journal of Machine Learning Research, 21(75):1–34, 2020.
Park and Kim (2020) G. Park and Y. Kim. Identifiability of gaussian linear structural equation models with homogeneous and heterogeneous error variances. Journal of the Korean Statistical Society, pages 1–17, 2020.
Park and Raskutti (2017) G. Park and G. Raskutti. Learning quadratic variance function (QVF) dag models via overdispersion scoring (ODS). The Journal of Machine Learning Research, 18(1):8300–8342, 2017.
Pearl (2009) J. Pearl. Causality. Cambridge university press, 2009.
Perrier et al. (2008) E. Perrier, S. Imoto, and S. Miyano. Finding optimal bayesian network given a super-structure. Journal of Machine Learning Research, 9(Oct):2251–2286, 2008.
Peters (2015) J. Peters. On the intersection property of conditional independence and its application to causal discovery. Journal of Causal Inference, 3(1):97–108, 2015.
Peters and Bühlmann (2013) J. Peters and P. Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1):219–228, 2013.
Peters et al. (2014) J. Peters, J. M. Mooij, D. Janzing, and B. Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15(1):2009–2053, 2014.
Peters et al. (2017) J. Peters, D. Janzing, and B. Schölkopf. Elements of causal inference: foundations and learning algorithms. MIT press, 2017.
Robins et al. (2008) J. Robins, L. Li, E. Tchetgen, A. van der Vaart, et al. Higher order influence functions and minimax estimation of nonlinear functionals. In Probability and statistics: essays in honor of David A. Freedman, pages 335–421. Institute of Mathematical Statistics, 2008.
Rosasco et al. (2013) L. Rosasco, S. Villa, S. Mosci, M. Santoro, and A. Verri. Nonparametric sparsity and regularization. The Journal of Machine Learning Research, 14(1):1665–1714, 2013.
Rothenhäusler et al. (2018) D. Rothenhäusler, J. Ernest, P. Bühlmann, et al. Causal inference in partially linear structural equation models. The Annals of Statistics, 46(6A):2904–2938, 2018.
Schölkopf (2019) B. Schölkopf. Causality for machine learning. arXiv preprint arXiv:1911.10500, 2019.
Shojaie and Michailidis (2010) A. Shojaie and G. Michailidis. Penalized likelihood methods for estimation of sparse high-dimensional directed acyclic graphs. Biometrika, 97(3):519–538, 2010.
Silander and Myllymaki (2006) T. Silander and P. Myllymaki. A simple approach for finding the globally optimal bayesian network structure. In Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence, 2006.
Singh and Moore (2005) A. P. Singh and A. W. Moore. Finding optimal bayesian networks by dynamic programming. 2005.
Spirtes and Glymour (1991) P. Spirtes and C. Glymour. An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review, 9(1):62–72, 1991.
Spirtes et al. (2000) P. Spirtes, C. Glymour, and R. Scheines. Causation, prediction, and search, volume 81. The MIT Press, 2000.
Tagasovska et al. (2018) N. Tagasovska, T. Vatter, and V. Chavez-Demoulin. Nonparametric quantile-based causal discovery. arXiv preprint arXiv:1801.10579, 2018.
Tsybakov (2009) A. Tsybakov. Introduction to nonparametric estimation. Springer Series in Statistics, New York, page 214, 2009. cited By 1.
van de Geer and Bühlmann (2013) S. van de Geer and P. Bühlmann. $\ell_{0}$ -penalized maximum likelihood for sparse directed acyclic graphs. Annals of Statistics, 41(2):536–567, 2013.
Wang and Drton (2020) Y. S. Wang and M. Drton. High-dimensional causal discovery under non-gaussianity. Biometrika, 107(1):41–59, 2020.
Wu and Fukumizu (2020) P. Wu and K. Fukumizu. Causal mosaic: Cause-effect inference via nonlinear ica and ensemble method. arXiv preprint arXiv:2001.01894, 2020.
Yu et al. (2019) Y. Yu, J. Chen, T. Gao, and M. Yu. Dag-gnn: Dag structure learning with graph neural networks. arXiv preprint arXiv:1904.10098, 2019.
Zhang and Hyvärinen (2009) K. Zhang and A. Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proceedings of the twenty-fifth conference on uncertainty in artificial intelligence, pages 647–655. AUAI Press, 2009.
Zheng et al. (2020) X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing. Learning Sparse Nonparametric DAGs. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 3414–3425, 2020.
Zhu and Chen (2019) S. Zhu and Z. Chen. Causal discovery with reinforcement learning. arXiv preprint arXiv:1906.04477, 2019.

	$\displaystyle\mathbb{E}\operatorname{var}(X_{j}\,\|\,X_{A})$	$\displaystyle=\mathbb{E}\big{[}\mathbb{E}[\operatorname{var}(X_{j}\,\|\,X_{B_{j}})\,\|\,X_{A}]\big{]}+\mathbb{E}\operatorname{var}[\mathbb{E}[X_{j}\,\|\,X_{B_{j}}]\,\|\,X_{A}]$
		$\displaystyle>\mathbb{E}\big{[}\mathbb{E}[\operatorname{var}(X_{j}\,\|\,\operatorname{pa}(j))\,\|\,X_{A}]\big{]}$
		$\displaystyle=\mathbb{E}\operatorname{var}(X_{j}\,\|\,\operatorname{pa}(j))$
		$\displaystyle=\mathbb{E}\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))$
		$\displaystyle=\sigma^{2},$

	$\displaystyle\mathbb{E}\operatorname{var}(X_{i}\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})=\sigma_{i}^{2}<$	$\displaystyle\sigma_{k}^{2}+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$
	$\displaystyle=$	$\displaystyle\mathbb{E}\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))+\mathbb{E}\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$
	$\displaystyle=$	$\displaystyle\mathbb{E}[\mathbb{E}(\operatorname{var}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})]+\mathbb{E}[\operatorname{var}(\mathbb{E}(X_{k}\,\|\,\operatorname{pa}(k))\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\operatorname{var}(X_{k}\,\|\,X_{\pi_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}})$

A polynomial-time algorithm for learning nonparametric causal graphs

Abstract

1 Introduction

Contributions

Related work

Notation

2 Background

Equal variances

Layer decomposition of a DAG

3 Identifiability and algorithmic consequences

3.1 Nonparametric identifiability

Theorem 3.1.

Example 1 (Causal pairs, (Mooij et al., 2016)).

Example 2 (Binomial models, (Park and Raskutti, 2017)).

Example 3 (Generalized linear models).

Example 4 (Additive noise models, (Peters et al., 2014)).

Unequal variances

3.2 A polynomial-time algorithm

3.3 Comparison to existing algorithms

Example 5.

4 Sample complexity

Condition 1 (Regularity).

Condition 2 (Identifiability).

Condition 3 (Estimator).

Theorem 4.1.

Corollary 4.2.

Choice of η\eta

Nonparametric rates

5 Experiments

Structure learning

Misspecification and sensitivity analysis

6 Discussion

Acknowledgements

Appendix A Reduction to order search

Appendix B Proof of Theorem 3.1

Lemma B.1.

Proof of Lemma B.1.

Proof of Theorem 3.1.

B.1 Generalization to unequal variances

Theorem B.2.

Example 6.

Proof of Theorem B.2.

Appendix C Proof of Theorem 4.1

C.1 A plug-in estimate

Proposition C.1.

Proof.

Corollary C.2.

C.2 A uniform bound

Proposition C.3.

C.3 Proof of Theorem 4.1

Proposition C.4.

Proof.

C.4 Technical lemmas

Lemma C.5.

Proof.

Lemma C.6.

Proof.

Lemma C.7.

Proof.

Appendix D Comparison to CAM algorithm

Proposition D.1.

Proof sketch of Proposition D.1.

Appendix E Experiment details

E.1 Experiment settings

E.2 Implementation and baselines

E.3 Metrics

E.4 Timing

E.5 Additional experiments

References

Choice of $\eta$