Optimal estimation of Gaussian DAG models

Introduction

A significant open question in the literature on structure learning is the optimal sample complexity of learning a directed acyclic graphical model. The problem of deriving upper bounds on the sample complexity for this problem goes back decades (Zuk et al., 2006; Friedman and Yakhini, 1996), and in recent years there has been significant progress (Ghoshal and Honorio, 2017a, 2018; Chen et al., 2019; Park and Raskutti, 2017; Park, 2018; Park and Park, 2019; Park, 2020; Wang and Drton, 2020; Gao et al., 2020; Gao and Aragam, 2021). Nonetheless, despite these upper bounds, a tight characterization of the optimal sample complexity is missing. This is to be contrasted with the situation for learning undirected graphs (UGs), also known as Markov random fields (MRFs), for which optimal rates were established approximately ten years ago (Santhanam and Wainwright, 2012; Wang et al., 2010), alongside similar results for support recovery in linear models (Wainwright, 2009a, b). In fact, this is unsurprising given the connection between these two problems via neighbourhood regression. Unfortunately, learning a directed acyclic graph (DAG) does not reduce to neighbourhood regression as it involves a more difficult order recovery step.

In this paper, we resolve this question for the special case of linear Gaussian DAG models with equal error variances. The identifiability of these models was established in Peters and Bühlmann (2013), and eventually led to the development of several polynomial-time algorithms under the equal variance assumption (Ghoshal and Honorio, 2017a, 2018; Chen et al., 2019; Gao et al., 2020). Nonetheless, it was not known whether or not any of these algorithms were optimal for this precise statistical setting. We will show that a variant of the EQVAR algorithm from Chen et al. (2019) is indeed optimal. This involves the derivation of new lower bounds and a novel analysis of the EQVAR algorithm that sharpens the existing sample complexity upper bound from $O(q^{2}\log d)$ to $O(q\log(d/q))$, where $q$ is the maximum number of parents in the DAG and $d$ is the number of nodes. This upper bound is optimal up to constants, and allows for the high-dimensional regime with $d\gg n$, where as usual $n$ denotes the sample size. Moreover, in Section 4, we extend this result to the case of general Gaussian models with known ordering. Our results also extend to more general identification assumptions (e.g. allowing for unequal error variances) as well as subgaussian error terms; see Remark 1.

As a problem of independent interest, we further compare the complexity of learning Gaussian graphical models (GGMs) and Gaussian DAG models under the equal variance assumption. Given the additional complexity of the order recovery problem in DAG learning, the folklore has generally been that learning DAGs is harder than learning UGs. Despite this folklore, few results are available to rigorously characterize the hardness of these problems on an equal footing (besides known NP-hardness results for both problems, see Srebro, 2003; Chickering, 1996; Chickering et al., 2004). The equal variance assumption gives us the opportunity to make an apples-to-apples comparison under the same assumptions. As we will show, the optimal sample complexity for both problems scales as $O(q\log(d/q))$. In other words, learning a DAG is statistically no harder than learning a GGM under the equal variance assumption. It is worth emphasizing that this comparison is purely statistical: The computational complexity of the algorithm we analyze is exponential in $q$ whereas learning GGMs can be done efficiently; see also Remark 2.

To the best of our knowledge, these are the first results giving a tight characterization of the optimal sample complexity for learning DAG models from observational data.

The rest of this paper is organized as follows: In the remainder of Section 1, we discuss related work and the problem setting. In Sections 2 and 3 we present our main results for learning equal variance DAGs. Then in Section 4 we consider the special case of known ordering, and in Section 5 make further comparisons with learning undirected GGMs. An illustrative simulation study is presented in Section 6 before concluding with some open questions in Section 7.

Given a directed graph $G=(V,E)$ with $|V|=d$ nodes, we make the following standard definitions:

The descendants $\operatorname{de}(k)$ to which $k$ has at least one directed path;

The nondescendents $\operatorname{nd}(k)=V\setminus\operatorname{de}(k)$;

The ancestors $\operatorname{an}(k)$ any of which has at least one directed path to $k$.

Given a random vector $X=(X_{1},\ldots,X_{d})$, we say that $G$ is a Bayesian network for $X$ (or more precisely, its joint distribution $P$), if the following factorization holds:

In this case, we abuse notation by identifying the random vector $X$ with the vertex set $V$, i.e. $V=X=[d]=\{1,2,\ldots,d\}$. We denote the class of all DAGs with $d$ nodes and at most $q$ parents per node (i.e. in-degree $\leq q$) by $\mathcal{G}_{d,q}$.

1 Related work

To provide context, we begin by reviewing the related problem of learning the structure of an undirected graph (e.g. MRF, GGM, etc.) from data. Early work establishing consistency and rates of convergence includes Meinshausen and Bühlmann (2006); Banerjee et al. (2008); Ravikumar et al. (2010), with information-theoretic lower bounds following in Santhanam and Wainwright (2012); Wang et al. (2010). More recently, sample optimal and computationally efficient algorithms have been proposed (Vuffray et al., 2016; Misra et al., 2020). Part of the reason for the early success of MRFs is owed to the identifiability and convexity of the underlying problems. By contrast, DAG learning is notably nonidentifiable and nonconvex. This has led to a line of work to better understand identifiability (e.g. Hoyer et al., 2009; Zhang and Hyvärinen, 2009; Peters et al., 2014; Peters and Bühlmann, 2013; Park and Raskutti, 2017) as well as efficient algorithms that circumvent the nonconvexity of the score-based problem (Ghoshal and Honorio, 2017a, 2018; Chen et al., 2019; Gao et al., 2020; Gao and Aragam, 2021). The latter class of algorithms begins by finding a topological sort of the DAG; once this is known the problem reduces to a variable selection problem. Our paper builds upon this line of work.

Other approaches include score-based learning, for which various consistency results are known (van de Geer and Bühlmann, 2013; Bühlmann et al., 2014; Loh and Bühlmann, 2014; Aragam et al., 2015; Nowzohour and Bühlmann, 2016; Nandy et al., 2018; Rothenhäusler et al., 2018; Aragam et al., 2019), but for which optimality results are missing. It is interesting to note that recent work has explicitly connected the equal variance assumption we use here to score-based learning via a greedy search algorithm (Rajendran et al., 2021). We also note here important early work on the constraint-based PC algorithm, which also establishes finite-sample rates under the strong faithfulness assumption (Kalisch and Bühlmann, 2007).

2 Problem setting

Although our results extend to more general settings, we focus on the special case of linear Gaussian Bayesian networks under equal variances. See Remark 1 for a discussion of generalizations. Specifically, let $X=(X_{1},\ldots,X_{d})$,

that is, each node is a linear combination of its parents with independent Gaussian noise. The variance $\sigma^{2}$ of each noise term is assumed to be the same; this is the key identifiability assumption that is imposed on the model.

More compactly, let $B=(\beta_{jk})$ denote the coefficient matrix such that $\beta_{jk}\neq 0$ is equivalent to the existence of the edge $j\to k$. Then letting $\epsilon=(\epsilon_{1},\ldots,\epsilon_{d})$ we have

The matrix $B$ defines a graph $G=G(B)$ by its nonzero entries, i.e.

Whenever $G$ is acyclic, it is easy to check that (1) holds, and hence $G$ is a Bayesian network for $X$. In the sequel we assume that $G$ is acyclic.

The following quantities are important in the sequel: The largest in-degree of any node is denoted by $q$, i.e.

The absolute values of the coefficients are lower bounded by $\beta_{\min}$, i.e.

Let the class of distributions satisfying the above conditions (3), (4), and (5) be denoted by $\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$. For any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ we have

This follows directly from (2) and $\operatorname{cov}(\varepsilon)=\sigma^{2}I$. Since the DAG is identifiable from the observational distribution, we denote $G(F)$ to be the DAG associated with the distribution $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$. Finally, we introduce the variance gap:

The proof of this lemma is a straightforward calculation; see Appendix E for details.

Both our upper and lower bounds can be generalized as follows: Although we assume Gaussianity for simplicity, everything extends to subgaussian families without modification. This is because the upper bound analysis relies only on subgaussian concentration, and the lower bounds easily extend to subgaussian models (i.e. since subgaussian also contains Gaussian as a subclass). Furthermore, the equal variance assumption can be relaxed to more general settings as long as $B$ can be identified by Algorithm 1. Examples include (a) the “unequal variance” condition from Ghoshal and Honorio (2018) (see Assumption 1 therein) and (b) if noise variances are known up to some ratio as in Loh and Bühlmann (2014). Moreover, both of these identifiability conditions include the naive equal variance condition as a special case, hence the lower bounds still apply. This implies more general optimality results for a wider class of Bayesian networks.

Algorithm and upper bound

We begin with stating the sufficient conditions on the sample size for DAG recovery under the equal variance assumption. Namely, we present an algorithm (Algorithm 1) that takes samples from a distribution $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ as an input and returns the DAG $G(F)$ with high probability. We first state an upper bound for the number of samples required in Algorithm 1 in Theorem 2.1.

For any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, let $\widehat{G}$ be the DAG return by Algorithm 1 with $\gamma=\Delta/2$. If

then $P(\widehat{G}=G(F))\gtrsim 1-\delta$.

The proof of this result can be found in Appendix A. The obtained sample complexity depends on the variance gap $\Delta$, which serves as signal strength, and covariance matrix norm $M$, which shows up when estimating conditional variances. Treating these parameters as fixed, the sample complexity scales with $q\log(d/q)$. The order of this complexity arises mainly from counting all possible conditioning sets. The proof follows the correctness of Algorithm 1, which consists of two main steps: Learning ordering and Learning parents.

Algorithmically, the first step is the same as Chen et al. (2019), however, our analysis is sharper: We separately analyze the estimation of each conditional variance directly rather than indirectly via the inverse covariance matrix. This leads to the improved sample complexity in Theorem 2.1. This step is where we exploit the equal variance assumption: The conditional variance $\operatorname{var}(X_{k}\,|\,C)$ of each random variable $X_{k}$ is a constant $\sigma^{2}$ if and only if $\operatorname{pa}(k)\subseteq C$ for any nondescendant set $C$. This implies that the variance of any non-source node in the corresponding subgraph would be larger than $\sigma^{2}$. Therefore, when all conditional variances $v_{kC}$ are correctly estimated with error within some small factor of the signal $\Delta$ (see Lemma A.1), identifying the node with the smallest $\sigma_{k}$ yields a source node in the underlying subgraph. Recall that $\sigma_{k}$ is the minimum variance estimation that node $k$ can achieve conditioned on at most $q$ nondescendants. Finally, recursively applying the above step leads to a valid topological sort.

In the second step, given the correct ordering, we use Best Subset Selection (BSS) along with a backward phase to learn the parents for each node. Note that BSS is already applied in the step 1.(c).i. of Algorithm 1 and the candidate set $C_{j}$ can be stored for each $\widehat{\tau}_{j}$, thus there is no additional computational cost. Again, when all conditional variances are well approximated by their sample counterpart $v_{kC}$, $C_{j}$ would be a superset of the true parents of current node $\widehat{\tau}_{j}$, otherwise the minimum would not be achieved. Meanwhile, removal of any true parent $i\in\operatorname{pa}(\widehat{\tau}_{j})$ from $C_{j}$ would induce a significant change in conditional variances, which is quantified by $\Delta$ as well. This is used to design a tuning parameter $\gamma$ in the backward phase for pruning $C_{j}$. Finally, we show the tail probability of conditional variance estimation error is well bounded to get the desired sample complexity in Lemma A.2.

When the true variance gap $\Delta$ is unknown, we can select the tuning parameter $\gamma$ according to the following theorem:

For any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, let $\widehat{G}$ be the DAG return by Algorithm 1 with tuning parameter

then $P(\widehat{G}=G(F))\gtrsim 1-\exp(-q\log(d/q))$.

The proof of this result can be found in Appendix A.3.

Lower bound

We will now present the necessary conditions on the sample size for DAG recovery under the equal variance assumption. Namely, we present a subclass of $\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ such that any estimator that successfully recovers the underlying DAG in this subclass with high probability requires a prescribed minimum sample size. For this, we rely on Fano’s inequality, which is a standard technique for establishing necessary conditions for graph recovery. See Corollary B.2 for the exact variant we use.

In Theorem 3.1, we state two sample complexity lower bounds: $q\log(d/q)/(M^{2}-1)$ and $\log d/\beta_{\min}^{2}$. Though the first one dominates when fixing other parameters as constants, the second one reveals the dependency on the signal strength $\beta_{\min}$. This can also be seen from the upper bound in Theorem 2.1 by replacing $\Delta=\beta_{\min}^{2}\sigma^{2}$ (cf. Lemma 1.1). We will present two ensembles for each bound. The first one is the whole set of sparse DAGs $\mathcal{G}_{d,q}$, and the second is the set of DAGs with only one edge, which is constructed to study the dependency on the coefficient $\beta_{\min}$.

For the first ensemble, we borrow the ideas from Santhanam and Wainwright (2012) to count the number of DAGs in $\mathcal{G}_{d,q}$. The only difference is we consider DAGs instead of undirected graphs. Also, it is easy to bound the KL divergence between any two distributions in this ensemble due to Gaussianity, which would lead to the bound $q\log(d/q)/(M^{2}-1)$. For the second ensemble, it is easy to count the size of this ensemble since we consider the DAGs with only one edge. Then all possibilities of any different pair of edges are analyzed to bound the KL divergence. This ensemble gives us the bound $\log d/\beta_{\min}^{2}$. The detailed proof can be found in Appendix B.

For comparison, Ghoshal and Honorio (2017b) previously established a lower bound for general Gaussian DAGs (i.e. without equal variances) of

which is a comparable lower bound. This is interesting since by restricting to simpler equal variance models (i.e. a smaller family), the problem should become easier, however, our analysis shows this is not the case. In particular, our lower bounds do not follow from previous work, and require a slightly different analysis as outlined in Appendix B.

Reconstructing a DAG from its ordering

The second step of Algorithm 1 may be of interest in its own right: Abstracted away, this step seeks to reconstruct a DAG from knowledge of its topological sort. We claim that the second step of Algorithm 1 is in fact sample optimal for learning the parents of each node (and hence all of $G$) given the true ordering of $G$ under more general assumptions.

Dropping the equal variance condition from $\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, define $\sigma^{2}_{k}\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{var}(\epsilon_{k})$ and let $\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma_{\max}^{2},M)$ denote the class of Gaussian distributions such that (3), (4), and (5) hold and

i.e. $\sigma^{2}_{k}$ is allowed to depend on $k$. Note that $\mathcal{F}_{d,q}(\beta_{\min},\sigma_{\max}^{2},M)\subset\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma_{\max}^{2},M)$. Furthermore, we modify the definition of the variance gap for $\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma_{\max}^{2},M)$ as follows:

Finally, given a known topological sort $\tau$ of $G$, let $\widehat{G}(\tau)$ be the DAG returned by the second step of Algorithm 1 with $\gamma=\overline{\Delta}/2$.

As with the rest of our results, these result extend to subgaussian models without issue. See Remark 1.

Using the second part of Lemma A.1, Lemma A.2 and following the proof in Appendix A.2, we have an upper bound on the sample complexity for recovering $G$ from its ordering:

For any $F\in\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma_{\max}^{2},M)$, given a valid topological sort $\tau$ of $G$, let $\widehat{G}(\tau)$ be the DAG returned by the second step of Algorithm 1 with $\gamma=\overline{\Delta}/2$. If

then $P(\widehat{G}(\tau)=G(F)\,|\,\tau)\gtrsim 1-\delta$.

The “given $\tau$” in the probability is to emphasize that the estimator has the access to the true ordering $\tau$.

Unsurprisingly, this approach of using best subset selection with a backwards phase is indeed optimal: We have a matching lower bound (up to constants).

then given the knowledge of true ordering $\tau$ of DAG $G$, for any estimator $\widehat{G}$,

The proof uses known lower bounds from the sparse support recovery literature (Wainwright, 2009b); see Appendix C for details.

This more general optimality result for the second step shows that it is only in the first step (learning parents) that the equal variance assumption is operational. Moreover, although it may be possible to improve the sample complexity of the second step for the smaller class $\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, since the sample complexity upper bound of the first step of Algorithm 1 matches the lower bound for recovering the whole graph, such improvements would not change the optimal sample complexity for learning $G$.

Comparison with undirected graphs

Our results on DAG learning under equal variances raise an interesting question: Is learning an equal variance DAG statistically more difficult than learning its corresponding Gaussian graphical model (i.e. inverse covariance matrix)? This is especially intriguing given the folklore intuition that learning a DAG is more difficult than learning an undirected graph (UG). In fact, it is common to learn an undirected graph first as a pre-processing step in order to reduce the search space and sample complexity for DAG learning (Perrier et al., 2008; Loh and Bühlmann, 2014; Bühlmann et al., 2014; Aragam et al., 2019). In this section we explore this question and show that in fact, at least in the special case of equal variance Gaussian models, the sample complexity of both problems is the same.

First, let us recall some basics about undirected graphical models, also known as Markov random fields (MRFs). When $X\sim\mathcal{N}(0,\Sigma)$ as in this paper, an MRF can be read off from the inverse covariance matrix $\Gamma=(\gamma_{jk})\mathrel{\mathop{\mathchar 58\relax}}=\Sigma^{-1}$. More precisely, the zero pattern of $\Gamma$ defines an undirected graph $U=U(\Gamma)$ that is automatically an MRF for $X$:

Wang et al. (2010) showed that the optimal sample complexity for learning a GGM is $n\asymp s\log d$, where

is the degree of $U$ or maximum neighborhood size, and Misra et al. (2020) developed an efficient algorithm that matches this information-theoretic lower bound. Given $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, let $U(F)$ be the undirected graph induced by the covariance matrix of $F$ (cf. 6). It follows that for any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ we can learn the structure $U(F)$ with $\Theta(s\log d)$ samples. Note that this sample complexity scales with $s$ instead of $q$.

Consider the DAGs $G_{1}$ and $G_{2}$ in Figure 1. In $G_{1}$, we have $q=O(d)$ since $T$ has $d$ parents, whereas in $G_{2}$ we have $q=O(1)$ since each $S_{k}$ has only one parent. Thus, we expect that learning $G_{1}$ will require $\Theta(d)$ samples and learning $G_{2}$ will require $\Theta(\log d)$ samples. By comparison, the UGs associated with each model, given by $U_{1}$ and $U_{2}$ have $s=d$, and hence if we use the previous approaches to learn each $U_{k}$ we will need $\Theta(d\log d)$ samples each. Of course, this is to be expected: One should expect that a specialized estimator that exploits the structure of the family $\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ to perform better.

2 Optimal estimation of equal variance GGMs

Example 1 shows that there is a gap between existing “universal” algorithms for learning GGMs (i.e. algorithms that do not exploit the equal variance assumption) and the sample complexity for learning equal variance DAGs. A natural question then is: What is the optimal sample complexity for learning the structure of $U(F)$ for any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$?

We begin by establishing a lower bound that matches the lower bound in Theorem 3.1 (up to constants):

To derive an upper bound for this problem, we use the well-known trick of moralization; see Lauritzen (1996) for details. Since $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, we can first learn the DAG $G=G(F)$ via Algorithm 1. Given the output $\widehat{G}$, we then form the moralized graph and define $\widehat{U}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{M}(\widehat{G})$. See Algorithm 2.

This approach is justified by a result due to Loh and Bühlmann (2014). First, we need the following condition:

Let precision matrix $\Gamma=\Sigma^{-1}=[\operatorname{cov}(F)]^{-1}$, $\Gamma_{ij}=0\Leftrightarrow\beta_{ij}=0\text{ and }\beta_{ik}\beta_{jk}=0$ for all $k\neq i,j$.

When we sample nonzero entries of $B$ from some continuous distribution independently, Condition 1 is satisfied except on a set of Lebesgue measure zero. For example, it is easy to check that the examples in Figure 1 satisfy this condition. Under this condition, moralization is guaranteed to return $U$:

If Condition 1 holds, then $U=\mathcal{M}(G(F))$.

Under Condition 1, we have the following upper bound, which matches the lower bound in Theorem 5.1:

Assuming Condition 1, for any $F\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$ let $\widehat{U}$ be the UG returned by Algorithm 2 with $\gamma=\Delta/2$. If

then $P(\widehat{U}=U(F))\gtrsim 1-\delta$.

The proof is straightforward by Theorem 2.1 and Lemma 5.2. This answers the question proposed at the beginning of this section: Under the equal variance assumption, learning a DAG is no harder than learning its corresponding UG.

It is an interesting question whether or not similar results hold without Condition 1; i.e. is there a direct estimator of $U$—not based on moralizing a DAG—that matches the sample complexity of learning an equal variance DAG?

To compare Theorem 5.1 with previous work, Wang et al. (2010) showed the optimal sample complexity is $\Theta(s\log d/\lambda^{2})$ where $\lambda^{2}=\min_{(s,t)\in E}\Gamma_{st}^{2}/(\Gamma_{ss}\Gamma_{tt})$, which is dominated by $\beta_{\min}^{2}$ when $\beta_{\min}$ is small. For equal variance GGMs, our result gives $\Omega(q\log(d/q)/M^{2}+\log d/\beta^{2}_{\min})$. Under Condition 1, $s$ is always greater than $q$, so our new lower bound is strictly smaller, along with a matching upper bound that shows the $s$ dependence for general GGMs is suboptimal for equal variance GGMs. Another related work (Cai et al., 2016) derives lower bounds on precision matrix estimation under certain matrix norms, which is distinct from the graph recovery problem we consider in this work. For comparison, put in our setting, their lower bound becomes $\Omega(s^{2}\log d)$, which again depends on $s$ instead of $q$.

Experiments

To illustrate the effectiveness of Algorithm 1, we report the results of a simulation study. We note that existing variants of Algorithm 1 have been compared against other approaches such as greedy DAG search (GDS, Peters and Bühlmann, 2013), see Chen et al. (2019) for details. In our experiments, we controlled the number of parents when generating the DAG, fix all noise variances to be the same, and sample nonzero entries of $\beta_{k}$’s uniformly from given intervals.

To generate random DAGs, we first randomly permute $[d]$ to obtain an ordering $\tau$. Then for each $j\in[d]$, we randomly draw a set $S$ of $q$ nodes from $\tau_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}$ and set

We then generate random nonzero coefficients according to $\beta_{j}\sim\operatorname{Rad}\times\operatorname{Unif}(0.5,1)$, where $\operatorname{Rad}$ is a Rademacher random variable. Finally, we generate data by the resulting Gaussian linear model:

We consider graphs with $d\in\{20,30,\ldots,90\}$ nodes and in-degree $q\in\{2,3,4\}$. For each setting, the total number of replications is $N=100$. For each replication, we generate a random graph and a dataset with sample size $n\in\{80,160,\ldots,560\}$. Finally, we report $\#\{\widehat{G}=G\}/N$ to approximate $P(\widehat{G}=G)$.

2 Implementation

We implement the learning ordering phase of Algorithm 1 using code from Chen et al. (2019)The code can be found at https://github.com/WY-Chen/EqVarDAG/blob/master/R/EqVarDAG_HD_TD.R., which inputs the oracle in-degree $q$ and outputs a topological sort. For the learning parents phase of Algorihm 1, we use the R package leaps (Lumley and Lumley, 2013) for Best Subset Selection with BIC.

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

3 Results

The results are shown in Figure 2. As expected, the probability of successfully recovering true DAG goes to one quickly across different settings of the number of nodes $d$ and maximum in-degree $q$. Since Best Subset Selection is computationally expensive, we are not able to examine higher dimensions systematically. Nonetheless, to test higher dimensional cases, we checked several (random) cases for $d=200,300,500$, and the results shows with 90% chance the DAG is successfully recovered for moderate sample sizes $n=480,560,900$.

Conclusion

In this paper, we derived the optimal sample complexity for learning Gaussian linear DAG models under an equal variance condition that has been extensively studied in the literature. These results extend to subgaussian errors under similar assumptions as well as more general models with unequal variances as long as the DAG remains identifiable by the proposed algorithm, which is easy to implement and simulations corroborate our theoretical findings. We also investigated the sub-problem of learning a linear DAG from its ordering and made comparisons with the classical problem of learning GGMs, showing the sample complexity of both problems is the same.

Finally, it would be of interest to generalize the results in Section 5 to more general families, i.e. beyond equal variances and its generalizations (Ghoshal and Honorio, 2018). This would require the derivation of new identifiability conditions, as in Gao and Aragam (2021) and Rajendran et al. (2021).

References

Appendix A Proof of upper bound

We first show that if all conditional variances are estimated sufficiently well, then Algorithm 1 is able to identify the true DAG.

If for all $k\in V$ and $C\subset V\setminus\{k\}$, $|C|\leq q$,

We start by showing that $\widehat{\tau}$ is a valid ordering for $G$, which is equivalent to saying $\widehat{\tau}_{j}$ is a source node of the subgraph $G[V\setminus\widehat{\tau}_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}]$ for all $j$. We proceed by induction. For $j=1$, it reduces to compare marginal variances.

Now we look at the second step of Algorithm 1, this step is to remove false parents from candidate set returned by Best Subset Selection. For any $j$, let $\widehat{\tau}_{j}=j$ for ease of notation. Given that $\widehat{\tau}$ is a valid ordering, $\operatorname{pa}(j)\subseteq\widehat{\tau}_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}$. We first conclude $\operatorname{pa}(j)\subseteq C_{j}$, otherwise there exists $C_{j}^{\prime}\subseteq\widehat{\tau}_{[1\mathrel{\mathop{\mathchar 58\relax}}j-1]}$ with $\operatorname{pa}(j)\subseteq C_{j}^{\prime}$ such that

Next we bound the estimation error tail probability:

For all $k\in V$ and $C\subset V\setminus\{k\}$, $|C|\leq q$,

Denote the covariance between $k$ and set of nodes $C$ at

Note that $\rVert\Sigma_{kC}\lVert,\rVert\Sigma_{CC}\lVert,\rVert\Sigma^{-1}_{CC}\lVert\leq\rVert\Sigma_{t}\lVert\leq\lVert\Sigma\rVert\leq M$.

Then the estimation error for conditional variance

The first inequality is by the triangular inequality, and the second simply bounds $\lVert\Sigma_{kC}\rVert$ by $M$. The third inequality introduces the estimation error of $\widehat{\Sigma}_{kC}$ and the final inequality replaces this with the estimation error of the full covariance matrix $\widehat{\Sigma}_{t}$. To set the RHS to be smaller than $\epsilon>0$, we consider three estimation errors. The first two can be controlled via standard sub-exponential concentration, whereas the third can be controlled via Theorem 6.5 from Wainwright (2019):

for some constants $A_{3},A_{4}$. The largest error is from

This is just another Gaussian covariance matrix estimation error, i.e.

for some constant $A_{5}$. Now we require all the errors to be bounded by $\zeta=\epsilon/M^{2}$ such that the conditional variance estimation error is within $\epsilon$. Thus

A.2 Proof of Theorem 2.1

Combine Lemma A.1 and A.2, we can have success probability:

The last equality is by $(\frac{d}{q})^{q}=(\frac{d}{q})^{q-1}\times(\frac{d}{q})\geq 2^{q-1}\times(\frac{d}{q})\geq q\times\frac{d}{q}=d$, thus $q\log(d/q)\gtrsim\log d$. In the end, replace $\epsilon$ by $\Delta/4$, solve for the sample size $n$ such that

we can have the desired sample complexity. ∎

A.3 Proof of Theorem 2.2

In the proof of Lemma A.1, denote the estimation error to be upper bounded by $\epsilon$, i.e. for all $k\in V$ and $C\subset V\setminus\{k\}$, $|C|\leq q$,

for the correctness of second phase to proceed. Therefore, let $\gamma=3\epsilon$ and require $\epsilon<\Delta/10$. Finally, set

Then we have failure probability bounded:

And to satisfy the requirement $\epsilon<\Delta/10$, we need sample size

Appendix B Proof of lower bound

Let’s start with recalling Fano’s inequality and its corollary under the structure learning setting. Let $\theta(F)$ be a parameter associated to some observational distribution $F$.

For a class of distributions $\mathcal{F}$ and its subclass $\mathcal{F}^{\prime}=\{F_{1},\ldots,F_{N}\}\subseteq\mathcal{F}$,

Set $\theta(F)=G(F)$, $\mathbf{dist}(\cdot,\cdot)=\mathbf{1}\{\cdot\neq\cdot\}$. One consequence of Lemma B.1 is as follows:

Consider some subclass $\mathcal{G}^{\prime}=(G_{1},\ldots,G_{N})\subseteq\mathcal{G}_{d,q}$, and let $\mathcal{F}^{\prime}=\{F_{1},\ldots,F_{N}\}$, each of whose elements is generated by one distinct $G\in\mathcal{G}^{\prime}$. If the sample size is bounded as

then the any estimator for $G$ is $\delta$-unreliable:

Thus the strategy for building lower bound is to find a subclass of original problem such that

Pairwise KL divergence between any two distributions is small.

Now we do some counting for the number of DAGs with $d$ nodes and in-degree bounded by $q$.

For $q\leq d/2$, the number of DAGs with $d$ nodes and in-degree bounded by $q$ scales as $\Theta(dq\log(d/q))$.

The proof construction is similar to Santhanam and Wainwright (2012). We can upper bound by number of directed graphs (DG), and lower bound by one particular subclass of DAGs.

For lower bound, we look at one subclass of DAGs. Suppose $d/(q+1)$ is an integer, otherwise discard remaining nodes. First partition $d$ nodes into $q+1$ groups with equal size $d/(q+1)$. Then for the first group, build directed edges from nodes in group $2,3,\ldots,q+1$ to group $1$, which requires $q$ permutations on $d/(q+1)$ nodes within one particular group. Then the nodes in group $1$ has exactly degree $q$. Similarly, for group $2$, build directed edges from group $3,4,\ldots,q+1$, which requires $q-1$ permutations on $d/(q+1)$ nodes. Therefore, for the subclass of DAGs generated in this way of partition, we have

many DAGs, any of which is valid DAG and has degree bounded by $q$. Then the cardinality

Thus the total number of DAGs scales as $\Theta(dq\log\frac{d}{q})$. ∎

B.2 Proof of Theorem 3.1

In this ensemble, we consider all possible DAGs with bounded in-degree. Note that by Lemma B.3, we know $N\asymp dq\log(d/q)$, it remains to provide an upper bound for the KL divergence between any two distributions within the class. For any two $F_{j},F_{k}\in\mathcal{F}_{d,q}(\beta_{\min},\sigma^{2},M)$, denote their covariance matrices to be $\Sigma_{j},\Sigma_{k}$. Due to Gaussianity, It is easy to see that

Therefore, we can establish the first lower bound that

Ensemble B

For this ensemble, we consider the DAGs with exactly one edge $u\to v$ and coefficient $\beta_{\min}$, denoted as $G^{uv}$. There are 2 directions and $d(d-1)/2$ many edges, so the cardinality of this ensemble would be $N=d(d-1)\asymp d^{2}$. Then denote the distribution defined according to $G^{uv}$ as $F^{uv}$, the log likelihood

and the difference between any two cases is

Then take expectation over $F^{uv}$ we get the KL divergence:

For fixed edge $(u,v)$, any other edges $(j,k)$ has relationship and corresponding KL divergence below:

$j\neq u,k\neq v$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}$

$j=u,k\neq v$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}$

$j\neq u,k=v$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}$

$j=v,k=u$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}+\beta_{\min}^{4}/2-\beta_{\min}$

$j=v,k\neq u$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}+\beta_{\min}^{4}/2$

$j\neq v,k=u$, $\mathbf{KL}(F^{uv}||F^{jk})=\beta_{\min}^{2}$

Among them, the largest KL between $F^{uv},F^{jk}$ is $\beta_{\min}^{2}+\beta_{\min}^{4}/2$. Therefore, we can conclude a lower bound

Appendix C Proof of Proposition 4.2

For simplicity we consider DAGs with $d+1$ nodes. We first recall a known lower bound for sparsity recovery: Consider the linear model $Y=\beta^{\top}X+\epsilon$ with $X\sim\mathcal{N}(0,\Sigma)$ and $\epsilon\sim\mathcal{N}(0,\sigma^{2})$. The support of $\beta$ is $S\subset[d]$ and $|S|=q$. Let $\beta_{\min}\mathrel{\mathop{\mathchar 58\relax}}=\min_{j\mathrel{\mathop{\mathchar 58\relax}}\beta_{j}\neq 0}|\beta_{j}|$. Then informally,

then with $q$ known, for any instance from the linear model and any estimator $\widehat{S}$ for $S$,

Now we adapt this result to our setting, $\omega_{bu}(\Sigma)\beta_{\min}^{2}$ is the variance explained by $X$ under regression model, thus upper bounded by $M\beta^{2}_{\min}$ when regarding $X$ as parents in DAG. Additionally, since every Gaussian with positive definite $\Sigma$ has a minimal I-map and given ordering, the parents can be read off through regression, i.e. the model class in Lemma C.1 is equivalent to the one generated in $\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma^{2}_{\max},M)$. For any estimator $\widehat{G}=\widehat{G}(\tau)$, denote $\widehat{\operatorname{pa}}(k)\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{pa}_{\widehat{G}}(k)$ for any node $k$. Then if

The first inequality is by relaxing the problem to simply finding the parents of the last node from all preceding nodes. The second inequality is because we can restrict at a sub-ensemble of $\overline{\mathcal{F}}_{d,q}(\beta_{\min},\sigma^{2}_{\max},M)$ whose last node of ordering has $q$ parents and maximum noise $\operatorname{var}(\epsilon_{\tau_{d}})=\sigma^{2}_{\max}$. The third inequality is because knowing the number of parents only makes the problem easier. The final inequality is by noticing the equivalence to sparsity recovery problem and applying Lemma C.1. ∎

Appendix D Proof of lower bound of GGM (Theorem 5.1)

We introduce two useful lemmas from Wang et al. (2010):

Consider a restricted ensemble $\widetilde{\mathcal{U}}\subseteq\mathcal{U}$ consisting of $N=|\widetilde{\mathcal{U}}|$ models, and let model index $\theta$ be chosen uniformly at random from $\{1,...,N\}$. Given the observations $X$, the error probability for any estimator $\widehat{U}$

The mutual information is upper bounded by $I(\theta;X)\leq\frac{n}{2}R(\widetilde{\mathcal{U}})$, where

with $a,b\to 0$ and $pb\to 0$, the determinant $\log\det A\approx pa$.

Finally, let’s consider three ensembles of UGs generated by DAGs. We describe the ensembles by showing how the DAGs generate the UGs.

In this first Ensemble, we consider an empty DAG, then add one edge from node $S$ to $T$ with linear coefficient $\beta_{\min}$. Specifically,

It remains to figure out the structure of covariance matrix and find out the corresponding determinants. Without loss of generality, let the first two nodes to be $S,T$, then the covariance matrix of any model (denoted as $j$th) is

It is easy to see that $\log\det\Sigma_{j}=\log(1+\beta_{\min}^{2}-\beta_{\min}\times\beta_{\min})=0$ for all models in this subclass. To compute the average $\bar{\Sigma}$, by symmetry, all diagonal and off-diagonal entries are the same respectively. For entries on diagonal, there are two situations: whether it corresponds to node $T$ or not. For off-diagonal entries, there are two situations: corresponds to edge $S-T$ or not. Different situations behave differently Table 1 with total counts $N$:

Thus we conclude the entries in $\bar{\Sigma}$:

Using Lemma D.3, we conclude $\log\det\bar{\Sigma}\asymp\beta^{2}_{\min}$, and invoking Lemma D.1 and D.2, we obtain a lower bound as

Ensemble B

Here we can adopt the same construction as the first ensemble for DAG in Appendix A, which applies analogously through Lemma B.1. Since the joint distribution remains to be the same, we have KL divergence upper bounded by $(M^{2}-1)d$.

For number of models inside this class, firstly we know that for a UG with degree bounded by $s$, there are $\Theta(ds\log{d/s})$ many UGs (Lemma 1(b) of Santhanam and Wainwright (2012)). By Lemma 5.2, $U=\mathcal{M}(G)$, so $q\leq s$, thus the number of UGs would be greater than $\Theta(dq\log(d/q))$, which leads to the same lower bound:

Appendix E Proof of Lemma 1.1

Immediate from the law of total variance: