Learning Gaussian  DAG Models without Condition Number Bounds
Authors: Constantinos Costis Daskalakis
ANTHIMOS VARDIS KANDIROS
Rui Yao
Venue: ICML 2025 poster (2025)
Keywords: Bayesian network
graphical model
Gaussian DAG models
condition number
Primary area: deep_learning->robustness
TL;DR: We provide algorithms for learning Gaussian DAG Models with sample complexity independent of the condition number.

Abstract: We study the problem of learning the topology of a directed Gaussian Graphical Model under the equal-variance assumption, where the graph has $n$ nodes and maximum in-degree $d$. Prior work has established that $O(d \log n)$ samples are sufficient for this task. However, an important factor that is often overlooked in these analyses is the dependence on the condition number of the covariance matrix of the model. Indeed, all algorithms from prior work require a number of samples that grows polynomially with this condition number. In many cases this is unsatisfactory, since the condition number could grow polynomially with $n$, rendering these prior approaches impractical in high-dimensional settings. In this work, we provide an algorithm that recovers the underlying graph and prove that the number of samples required is independent of the condition number. Furthermore, we establish lower bounds that nearly match the upper bound up to a $d$-factor, thus providing an almost tight characterization of the true sample complexity of the problem. Moreover, under a further assumption that all the variances of the variables are bounded, we design a polynomial-time algorithm that recovers the underlying graph, at the cost of an additional polynomial dependence of the sample complexity on $d$. 
We complement our theoretical findings with simulations on synthetic datasets that confirm our predictions.

Learning Gaussian DAG Models without Condition Number Bounds
Constantinos Daskalakis
Vardis Kandiros 1
Rui Yao
Abstract
of observations (Duncan, 2014). In particular, this is one of the fundamental goals of causal inference, where the empha- sis is discovering cause-and-effect relationships between variables. Perhaps the most popular way of approaching this problem is through the formalization of Directed Acyelic Graph (DAG) models, or Bayesian Networks, where vari- ables are sampled sequentially, with the order induced by a DAG G(Lauritzen, 1996; Wainwright et al., 2008). These models have been widely adopted in the social sciences, starting from their introduction in the seminal work (Pearl et al., 2000). They have also found applicability in fields as diverse as human geography (Malioutov et al., 2006; Ao & Chang, 2020), the environment (Grace et al., 2007). and ecology (Byrnes et al., 201 1). The problem of determining G from observations of the variables is referred in the literature causal structure learning.
We study the problem of learning the topology of a directed Gaussian Graphical Model under the equal-variance assumption, where the graph has n nodes and maximum in-degree d. Prior work has established that Od samples are sufficient for this task. However, an important factor that is often overlooked in these analyses is the dependence on the condition number of the covariance matrix of the model. Indeed, all al- gorithms from prior work require a number of samples that grows polynomially with this condi- tion number. In many cases this unsatisfactory, since the condition number could grow polyno- mially with n, rendering these prior approaches impractical in nigh-dimensional settings. In this work, we provide an algorithm that recovers the underlying graph and prove that the number of samples required is independent of the condition number. Furthermore, we establish lower bounds that nearly match the upper bound up a d d-factor, thus providing an almost tight characterization of the true sample complexity of the problem. More- over, under a further assumption that all the vari- ances of the variables are bounded, we design a polynomial-time algorithm that recovers the un- derlying graph, at the cost of an additional poly- nomial dependence of the sample complexity on d. We complement our theoretical findings with simulations on synthetic datasets that confirm our predictions.
In the general case, some popular approaches for solving this problem include the PC algorithm (Spirtes et al., 2001) and Greedy Equivalence Search (GES) (Chickering, 2002b). However, these often require strict assumptions, such as strong faithfulness (Uhler et al., 2013), or encounter compu- tational issues. Indeed, general, the problem of recovering the true Bayesian Network is NP-hard (Chickering, 1996).
The problem becomes more structured using the formal- ization of Structural Equation Models (SEMs), where each variable is equal to a function of its parents with some added independent noise. Identifiability for these models has been established when the functions are sufficiently non-linear or when the noise is non-Gaussian (Peters et al., 2014; Shimizu et al., 2006). For linear Gaussian SEMs, which is the fo- cus of this work, finding the graph in general is impossi- ble. Identifiability was established in (Peters & Bühlmann, 2014), under the additional assumption that all e; have the same variance (equal variances assumption). A subsequent flurry of works (Ghoshal & Honorio, 2017b;a; 2018; Chen et al., 2019: Gao al., 2020: 2022) established upper and lower bounds the number of samples needed (sample complexity) for inferring the graph in the high dimensional sparse setting where the maximum in-degree of any node is bounded by d and the number of nodes n >> d. Often, the goal in these works is to obtain sample complexity that scales as (d no...
1. Introduction and Background
common problem that arises in many scientific disciplines is inferring the underlying dependency structure of a vector
'Equal contribution 'Department of EECS, Massachusetts Institute of Technology, Cambridge, Massachusetts. United Constantinos Daskalakis <costis mrostis@csail.nne misail.mit.edu2,, Vardis Kandiros <kandiros mit.edu>, Rui Yao <rayyao mii..duu>.
States. Correspondence
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
However, an important factor that is often overlooked in

---

Learning Gaussian DAG Models without Condition Number Bounds
prior work is the role of the condition number x of the covariance matrix of observations. Indeed, all previous anal- yses typically assume some type of condition that implies bound on K, independent of n. However, in many interest- ing examples and simple examples, K scales polynomially with the number of nodes n. This would imply that all al- gorithms in prior work require poly(n) samples, which is problematic in nigh-dimensional settings. For example, in the directed path graph the tree graph, grows polynomi- ally with n, but the structure can clearly be identified easily using (log n) samples. Thus, a natural question that arises is whether it is even information theoretically possible to design an algorithm that is truly high-dimensional, with a sample complexity independent of the condition number.
more detailed technical comparison of our results with prior work in Section 2.4. Finally, we give a proof sketch for the main results in Section 3 and simulation results in Section 4.
1.1. Prior work
There is a long line of work that focuses on recovering the underlying directed acyclic graph (DAG) in structural equation models. A popular approach is the PC algorithm (Spirtes et al., 2001), which runs in polynomial time if the graph is sparse. (Robins et al., 2003) shows that it is point- wise but not uniformly consistent under the classical faith- fulness assumption on the underlying distribution. (Zhang & Spirtes, 2012; Kalisch & Bühhlann 2007) show uniform consistency under the strong faithfulness assumption, which can be sometimes restrictive (Uhler et al., 2013). Other ap- proaches include score-based model Spiieeto et al., 1993; Heckerman et al., 1995; Chickering, 2002b;a), where searching for a graph can be cast as an optimization problem.
In this work, we answer this question in the affirmative. In particular, we design and analyze algorithm that achieves graph recovery with sample complexity independent of the condition number while at the same time retaining the opti- mal (d log n) scaling, We identify a key quantity TGG) of the graph (see (2)) that determines this sample complexity, which is always upper bounded by the condition number k. Furthermore, we provide lower bound for the sample com- plexity of any algorithm, which also contains r and matches our upper bound up to single factor of d (this d-factor dis- crepancy of the upper and lower bound is a known limitation of all prior lower bound approaches for Markov Random Fields (Wang et al., 2010; Misra et al., 2020; Kelner et al., 2020)). Thus, we obtain an almost tight characterization of the optimal sample complexity of the problem. Another design an efficient al-
Closer to our work is line of research that considers learn- ing structural equation models (SEMs). Various works show identifiability when the functions are sufficiently non- linear (Peters et al., 2014) or when they are linear and the noise is non-Gaussian (Shimizu et al., 2006). In general, when the noise is Gaussian, it impossible to determine the direction of the edges, which makes the problem non- identifiable. (Peters & Bühlmann. 2014) proved identifia- bility under the equal variances assumption, followed by an investigation on the optimal sample complexity (Ghoshal & Honorio, 2017a;b; Chen et al., 2019; Gao et al., 2020; Gao & Aragam, 2021: Gao et al., 2022). In particular, in (Gao et al., 2022) they provide an algorithm that requires (d log samples. As we mentioned earlier, all these prior analyses assume that the condition number is bounded, and the sample complexity is implicitly dependent on that with a polynomial rate. Another notable work is (Gao et al., 2023). This paper focuses on support recovery in a linear regres- sion setting where the underlying design satisfies additional structural constraints. The authors provide an alternative ap- proach Best Subset Selection and show it performs better in a setting where there is no generalized path cancellation.
natural question is whether gorithm that solves the problem without condition number assumptions. To this end, we provide such an algorithm, which works under the natural assumption that the max- imum variance of a node is bounded by a constant. We also show that this is a strictly weaker assumption than the boundedness of the condition number.
In terms of techniques, the problem of learning DAG is information-theorrticll different in nature from that of learning undirected graphical models due to the combinato- rial nature involved in finding the right topological ordering of the nodes. Despite that, our approach brings together a variety of tools from the literature on learning directed as well as undirected graphical models. We hope this serves as an inspiration for further fruitful investigations that would uncover connections between the two problems.
Finally. another popular branch of the literature consid- ers learning undirected Markov Random Fields (MRFs) (Bresler, 2015; Vuffray et al., 2020; 2016; Lokhov et al., 2018; Klivans & Meka, 2017: Wu et al., 2019; Misra et al., 2020). In general, this problem possesses an inherent con- vexity and can be solved by neighborhood regression, which is why it is considered easier than the directed case. In the special case where the underlying distribution is Gaus- sian, the so-called undirected Gaussian Graphical Model (GGM), a variety of algorithms have been proposed when
The rest of the paper is organized as follows. Section 1.1 discusses prior work related to the problem, and Section 1.2 contains some preliminary definitions and facts about the model. The main Algorithm that achieves information- theoretic discovery is given in Section 1. with the matching lower bound following in Section 2.2. An efficient Algo- rithm is presented in 2.3. For completeness, we provide
2

---

Learning Gaussian DAG Models without Condition Number Bounds
the underlying graph is sparse, based on estimating the co- variance matrix(Meinshausen & Bühlmann, 2006; Yuan &
following relation holds
Lin, 2007; Cai et al., 2011; Anandkumar et al., 2012: Cai et al., 2016; Johnson et al., 2012; Friedman et al., 2008). Under minimal assumptions, (Wang et al., 2010) proved a sample complexity lower bound of ( (miin) log ri), where 3min is the minimum normalized edge strength of the model. (Misra et al., 2020) provides an algorithm that suc- ceeds with ((d//3minn (miu) log n) samples, and (Kelner al., 2020) designs computationally efficient alternative with similar sample complexity under the additional assumption of valk-summability (Malioutov et al., 2006). An interest- ing folklore question is to compare the complexity of the problems of identifying the undirected and directed graphs of Gaussian models(Gao al., 2022). Our work can thus be seen as establishing a separation between these two tasks, which are govemned by different quantities, Aman in the case of undirected and in the case of directed models. Finally, a related line of work studies Gaussian models with the ad- ditional assumption that the topology is a tree, see e.g. et al., 2010; 2011).
B
where X is the vector of random variables indexed by the same topological order, and is a standard n-dimensional Gaussian vector. This has the straightforward implication that the covariance E and inverse covariance O of vector X can be written as
It is clear that a simply controls the scale of all variables in the problem and doesn't affect structure recovery. Therefore, for simplicity in the remainder of this work, we assume that 1. If was unknown, our results extend straightfor- wardly. Algorithm 1 has the same sample complexity but with an added first step to estimate the smallest variance. The number of samples in Algorithm would scale by o -2, More details can found in Supplementary E.3. We make the following assumption throughout this work. Assumption 1.1.
1.2. Preliminaries and Problem Settings
min
bbyl buin 0
SEV SEV.EEaa()
Let G = (V, be a Directed Acyelic Graph (DAG) on n nodes, where the maximum in-degree of a node is bounded by d. For j E V. we sometimes write -> j if (i..) E F directed edge, As discussed in Section 1, G can be thought of as the fixed ground truth graph that we seek to identify. We denote by pa(i)= = [j E V jj - i} the parent- set of node i. By assumption, |pa(t) < d for all i E V. A topological order of the nodes in G is any permutation of the nodes T: : V + [n]] such that i + j implics (ii) < (f). Since G does not have directed cycles, this ordering always exists but not necessarily unique. With each node in G, we associate random variable X; E TR. A vector of random variables (X1, X2........ Xn) E IR", is generated according to DAG G if the following structural equation holds
Assumption 1.1 places a lower bound on the minimum strength of an edge between two nodes. This lower bound
is essential information thaattiinnttoorticlllly.... if we want to detect the presence of each edge. In this paper, we assume bmin and the sparsity d are known to the learner. We discuss ex- tending our result to the case when bmin and d are unknown in the Appendix E.3. More specifically, we can perform the full Algorithm if mmin is unknown but the guarantee alters. We can also perform the first phase of our Algorithm if d is unknown.
We also define for a graph G the following quantity 7G)) which is the maximum sum of squares of the coefficients for the out-edges of a node.
x,= Σ b t jEM(F)
=GG = 1 max Σ
6 1-j-
3
We will see in the sequel how r(G) arises naturally in both the upper and lower bounds on the sample complexity of the problem.
where E; ~ o3) is i.i.d Gaussian noise of equal variance. In particular, if a node i has no parents, X, is just (O0, 1). The above model is the so-called equal variances model. We note that when this assumption is not satisfied, might be impossible identify the DAG, even up to the Markov Equivalence Class. For details see 1. For a set T, we denote X, to be the vector of (X,) for i E 7.
2. Results
In Section 2.1 we present algorithm based neigh- borhood selection and analyze its sample complexity. In Section 2.2 we present an information theoretic lower bound that matches this sample complexity up to a factor of d. Fi- nally, in Section 2.3 we provide a polynomial time algorithm with slightly worse sample complexity.
We can collect all coefficients by: in a matrix B E R_ so that B, = dij is non-zero only if i + j. If we index the rows and columns of B according to a topological or- dering, we can think of B as being upper triangular, and the
3

---

Learning Gaussian DAG Models without Condition Number Bounds
2.1. Information Theoretic Recovery
error for these matrices, which leads sample complexity that depends on condition number assumption. Instead, in our analysis we take advantage of the relation between this conditional variance and the mean squared error of regress- ing the node with the previous nodes in the ordering, which results in significant gains in sample complexity.
The algorithm is split into two phases, as shown below.
Algorithm 1 Estimating the DAG Information Theoretically
input Samples (X _x""),...,x(ny" X(m)
output Topology G
Thus, for Phase I, we still calculate the conditional vari- ances as in (Gao et al., 2022) and (Chen et al., 2019). We observe that for any node i and any set J, if J contains all the parents of i, then Var[X; [X,] = 1. Otherwise, the conditional variance is larger than 1 + b Thus, there is a gap between these two cases (see Proposition A.1) which helps us distinguish between them.
1: Phase I: Finding the Topological Ordering
2: T=
3: Find all X; with Var(X,) S b1+ /2) and put i E T
4: while [7" < 71 do
5: . = 0
6: for E [T] do
7:
for J S T, [J| = d do
8:
33 + Regress X; on X,
Phase 1 of the algorithm outputs an a list T, a valid topolog- ical order of all the nodes, In line 2, we initialize T as an empty list. We know that the first node in the topologiaal order has no parent, thus it has variance of 1, the smallest variance among all random variables. Therefore, in line 3, we choose such that X, has a small empirical variance to be the first in the list.
MSE + 1/2 E, x ()) - 2
9:
10:
2 then if MSE - 1+ 6 /2 T - T.append(i)
11:
end if 13: end for
12:
14: end for
Let now discuss lines 4 to 15 of Algorithm 1. Suppose, inductively, that we have identified the correct T for each iteration. In that case, all the nodes in T appear earlier in the topological order than V\T. This means that in the next iteration, the node that comes right after T in the topological order has all its parents in T. We call that node i. Therefore, the conditional variance, Var [X, [Xr], is (see Proposition A. and we can choose the one with smallest conditional variance accordingly.
15: end while
16: Phase 2: Finding the parents
17: for E do
18:
S 5 nodes in T before
19:
for S. J| = min([S),d) do
20:
PASSED + True
21:
for K E ()) S min[[||] [J], dd do
22:
BroK < Regress X; on XJuk
23:
if E K : 18 2 bmin/2 then
24:
PASSED t False, break.
However, we cannot directly calculate the conditional vari- ance on Xy or even estimate it with bounded accuracy by vanilla OLS. If we do so, we need O(n) samples, which is larger than our target O(lg n) samples, to esti- mate the conditional variance of X; given Xr. Neverthe- less, we can take advantage of the sparsity of the graph, i.e. [pa(i)] I d. By conditional independence, Var[X, (araX))Xpato] = Var[X,
25:
end end for
26:
27:
ir PASSED = True then
28:
-- regress X; on X,
29:
pa(i) + (i E J: 18,1 2 2n) } bbreak.
30:
end if
have 1. Therefore, instead of calculating the conditional variance Var[X,|X7] we can calculate Var [ar[X,]X]]] for all < d-element subsets J of T' Hence, lines 7 to 13 do the following: For all i, we loop through J C |,|J| S d to find the subset that makes the empirical conditional variance smaller than 1 + b /2. If we find such i, we append to topological order. We repeat this procedure until there is node outside T.
31:
end for
32:
end for
In Phase 1, we aim to find the topological order. A fun- damental property of the equal-variances model is that the variance increases along directed paths. Thus, the strategy is to start with the node of smallest variance. Intuitively, we want to keep adding the node with the smallest conditional variance, conditioned on the nodes we have already added. One approach for computing this variance is by using the conditional variance formula for Gaussians, as is done in prior work (Chen et al., 2019; Gao et al., 2022). One impor- tant challenge in our setting though is that the expression for this conditional variance involves covariance matrices with arbitrarily bad condition numbers. Prior analyses in (Chen al., 2019; Gao et al., 2022) bound the estimation
After finding topological ordering in the first phase, the second phase aims to recover the parents of each node, in order to completely determine the topology. Here, we draw inspiration from a technique developed in (Misra et al., 2020) for undirected Graphical Models. For any vertex i, we would like to figure out which of the candidate subsets of size at most d that come before it is the correct parent set.
4

---

Learning Gaussian DAG Models without Condition Number Bounds
The only guarantee we have in this sample-starved setting is that linear regression with O(d) covariates has small error. Thus, for any such candidate parent set A, we perform series of regressions of i with covariate set A U B, where B ranges over all subsets of d nodes before i again. If for at least one of these sets B, some regression coefficient of a variable in X; E B comes out sufficiently large, then we know that A does not contain all parents and hence reject it. The reason is that conditioned on the parent set, X; should be independent of all other variables that come before it in the ordering. On the other hand, if for all sets B, the only large coefficients come from variables in set A, then we know that A contains all parents of i. We then pick the parents of i as the variables of A with large coefficients. A detailed description of this procedure given in Algorithm 1. We remark that, for the second phase of finding the parents, it is possible to use as an alternative the Best Subset Selection algorithm. Even though the authors in (Gao et al., 2023) do not consider the equal variances model, it is possible to use their analysis about support recovery together with some additional steps to obtain similar sample complexity guarantees for the second phase of finding the parents.
relevant quantities for finding the topology are n, d, bmin with dependence K. An additional relevant quan- tity for the sample complexity of determining the parent set is rG).. We note here that this is always better than a condition number bound, as the following Lemma suggests.
Lemma 2.2. For any probabilistic model of the form (17) with DAG G. we have (G) S K.
The proof of this lemma can be seen in Appendix D. Thus, our bound stronger than the ones obtained using condition number assumptions (Chen et al., 2019; Gao et al., 2022). Furthermore, there are also cases where T is constant and k grows exponentially with n. For example, if we have the DAG X1 + X2 + + X,, and 6,,4+1 = k > 1, then the condition number is at least k=-1 but T(G) is constant. A related construction involving trees with similar guarantees given in the proof of Lemma 2.5.
For the bound rG),, in Section 2.2, we show that this de- pendence 7 is necessary. For more detailed discussion and comparisons of this result with that of prior work, we refer the reader to Section 2.4.
2.2. Information Theoretic Lower Bound
The next result provides an upper bound on the number of samples required by Algorithm 1 to find the correct topol- ogy with high probability (proof: A.) Note that Algorithm requires knowledge of the parameters a, d, bmin- In Supple- mentary E.3, we present an adaptive version of this algo- rithm with similar sample complexity guarantees.
In Section 2.1 we presented an algorithm that suc- ceeds with high probability as long as it is given O (max(bb= min =bm2n )d ldloog//d)) samples. We now present lower bound on the sample complexity of any algorithm, which matches the upper bound up to a factor of d. Thus, we obtain fine-grained characterization of the optimal sample complexity for this problem.
Theorem 2.1. Suppose we run Algorithm/ I using in- dependent samples generated from DAG G according to (17). If Assumption 1.1 holds, then there exists an absolute constant > 0, such that the following guarantees hold.
Our general strategy will be to define family of distri- butions as in (17), each indexed by a different topology and then apply Fano's 's inequality (Yu, 1997). Each family will instantiate some part of the difficulty of the problem, which corresponds to some term in the sample complexity. Concretely, we can prove (in Appendix B) the following.
Phase of Algorithm 1 succeeds in finding a correct topological ordering of the nodes with probability at least 8, provided that
Theorem 2.3. We can construct family of DAGs g satis- fying Assumption 1.1 with the following property: suppose one DAG G E G is picked uniformly at random from G and m independent samples (1), (1),,..,'""".... X (mm) are generated according to G. Then, any estimator e Rmm + G that takes as input these samples satisfies
m 2 C JTT
(dlog(n/d) + log(L/o))
Provided that Phase 1 of Algorithm I succeeds in find- ing a correct ordering, Phase finds the correct parent set for each node with probability at least 1 8. pro- vided that
TP (xiL) x(mm)
1
C. (C) (dlog(n/d) log(1/1)) m
12
bmin
as long as m & C mmax(b (iii Tb log(n). where C' >>> 0
Notice that we obtain separate characterizations for the sam- ples required to find the correct topological ordering and the correct parent sets. We remark that for the task of finding the ordering, it suffices to assume that the noise is sub- gaussian (for details see Supplementary E.4). The only
is an absolute constant.
Thus, Algorithm 1 is minimax optimal up to a factor of d. Each of the two terms in the max reflects the difficulty of finding the correct topology and finding the parent set,
5

---

Learning Gaussian DAG Models without Condition Number Bounds
respectively. Hence, essentially characterize the com- plexity of these two problems separately up to a d-factor. This is explained in more detail in the Supplementary Mate- rial. For a more detailed comparison of this result with prior work see Section 2.4.
error of the Lasso to determine the next node in the ordering. Afterwards, to find the parents of each node, we regress on nodes that appear before it in the topological ordering. The details appear in Algorithm 2, where À is the regularization parameter.
2.3. Efficient Algorithm with Bounded Variance
Algorithm 2 Efficient algorithm for recovering the topology
input Samples (X ( ,...,, x (m) output Topology G
The runtime of Algorithm 1 is O(n2++2), which even for small values of the maximum degree becomes impractical. We are thus in search of a computationally efficient alter- native, which will also enjoy favorable sample complexity guarantees. One important observation about structure learn- ing is its similarity with Sparse Linear Regression (SLR), where one is tasked with estimating the regression vector in linear regression problem, assuming at most d of its n entries are non-zero. Information theoretically. (d log n) samples are sufficient to estimate the vector. However, it is conjectured that no polynomial time algorithm can succeed under general designs using o(n) samples (Kelner et al., 2024). The connection with our problem comes from the fact that finding the parents of each node basically involves solving a SLR problem where the non-zero coefficients correspond to the parents of the node. Given this close con- nection, it seems unlikely that a polynomial time algorithm succeeds in finding the DAG using o(n) samples, at least in the general case.
1: T 1: Var S (1 + 6 /2
while [T] < n do
.M =0
4:
for E [m]] ( ( do
5:
Choose À O v Logg////)))1/// vm)
8 - regression of X; on Xr
MSE + L, [ (3) x 2
7:
8:
if MSE S 1 1 mmnn/ / then
9:
M - M {+}
10:
pa(i) < ( J: 18,1 2 Umin\ 22
11: end if
12: end for
13: T +- T.append(M)
14: end while
Below we present an upper bound on the sample complexity of Algorithm (proof is in Appendix C.)
We thus need pose extra assumptions in order to make the problem tractable. One common assumption, which is also present in prior work (Ghoshal & Honorio, 2017b), is that all individual variances are bounded by a constant.
Theorem 2.6. Under Assumptions 1.1 and 2.4, there is an absolute constant C such that Algorithm 2 run with m samples recovers the correct DAG with probability at least 1, as long as
Assumption 2.4. For all i, Var (Xx) < R
m 2 C(R2 log(n/8) -3.73.66.. -mmin)
This can be seen mild assumption, as it rules out patho- logical cases where the variance of some nodes grows with the size of the DAG, a situation that could be considered unrealistic in most applications. It is also worth noting that this assumption is strictly weaker than a condition number bound. In particular, we can prove that a condition number bound implies bound on the maximum variance, but the converse is not necessarily true (proof in Appendix D.2).
Moreover, Algorithm 2 runs in time poly(n, (, m. T, R).
A comparison with prior work follows in the next Section.
2.4. Discussion and Comparison
For completeness of presentation, we now offer a more detailed technical comparison of our results with the most relevant prior work.
Lemma 2.5. (i) Suppose we have variables X1,..., Xn
generated according to (17) and let K be the condition number of their covariance matrix. Then Var (Xi) <
Regarding upper bounds on the sample complexity, (Ghoshal & Honorio, 2017b) obtain a sample complexity of O 1o44 log n)) where à is the size of the maximum Markov blanket, which could be much larger than the maximum in-degree. Furthermore, they need assume that the in- finity norm of θ is bounded, which is form of condition number assumption that we do not require. In (Chen et al., 2019) they obtain an improved (d log n)) dependence, but they also need an upper bound on the maximum variance of each node in the graph. Finally, (Gao et al., 2022) intro- duce an algorithm that needs (M log n) samples, which
k. f.for all i.
(ii) For every a E (0, 1), there exists sequence of DAGs {Gn)r In indexed by the number of nodes, with maximum variance Rn, and condition number Kny such that R., = O(1) but K., = ).
Our strategy for finding the DAG is to use Lasso to discover the neighborhood of each graph. We use the mean squared
6

---

Learning Gaussian DAG Models without Condition Number Bounds
has the optimal dependence on d and n, but also depends on M, which is upper bound on the square root of the condition number of K of Σ. In contrast, Algorithm I needs ( min.. Tb 2 )d log(n/d)) d)) samples, which is also op- timal in the dependence on n and d, but only depends on T. which is smaller than A by Lemma 2.5.
to the PC algorithm (Spirtes et al., 2001). Although PC algorithm does not require the structural equation as Equa- tion (17), the faithfulness assumption is not always satisfied because of the path cancellation. We also performed ex- perimental simulation for the PC algorithm, and the result shows that the PC algorithm behaves worse than Algorithms 1. 2, and (Gao et al., 2022) even for the case of small n, d. Further comparison and discussion can be found in Section E.2 in the supplementary material.
The quantity т is information theoretically necessary, as Theorem 2.3 implies an (( min) log(n)) lower bound. The only lower bound in prior work that is directly comparable to ours is the one in (Gao et al., 2022). Using our notation, Theorem 3.1 from their work establishes the d log(n/d)
Finally, Phase 1 of Algorithm I can be extended to the case when E; are only centered sub-Gaussian with equal variance (see Section E.4).
lower bound 52 ( M--11 Dumin n)). It is clear that our bound has the improved min factor compared to b of (Gao et al., 2022). At first glance, it seems that their bound 22( (dlog(n/d)) is stronger by a factor. However, we M=-1
3. Proof Sketch of the Analysis
can show (see Appendix D.6.) that 2 -1 1 2 db? min always holds, which means that the bound in (Gao et al., 2022) re- duces to log which is strictly worse than the lower bound b2 in Theorem 2.3. In fact, this d-factor discrepancy between upper and lower bounds is a known limitation in the lit- erature on learning undirected GGMs (Wang et al., 2010: Misra et al., 2020; Kelner et al., 2020). Thus, improving the lower bound for our problem could have significant impli- cations for other related questions as well. Also, (Ghoshal & Honorio, 2017a) establishes general lower bounds for Bayesian Networks where the conditional distributions are exponential families. However, when instantiated in the equal variances model, these results scale as & &logn where Wmnxx is the maximum L2-norm of the coefficients vector of the parents of a node. Thus, this bound is of the order of log n and fails to capture the dependence on T.
Before presenting the proof sketch for Theorem 1, we would like to highlight the technical novelty that leads to an al- gorithm with sample complexity that is independent of the condition number. The general technique is to use the gap between conditional variances if some parents are missing from the conditioning set. For Algorithm 1, we conduct a more detailed analysis of the conditional variance to achieve better sample complexity. Specifically, we use concentra- tion properties of the OLS error for both Phase I and Phase 2 in Algorithm 1. For Phase 1, the analysis of the conditional variance is sufficient to make the claim in Theorem 1, which is a bound already free of conditional number. However, for Phase 2, the variance gap only leads to Oddoog(n) IIA)" Pumi
sample complexity, where the ratio T(G)²
is not informa- tion optimal. Therefore, we shift the method to additionally examine the OLSS coefficients (Misra et al., 2020) to distin- guish the correct parent set. Therefore, we can then sharpen the sample complexity to (log(n)) /dd = by using this method. For Algorithm 2, we employ LASSO regression and also utilize the technique of examining the MSE and the coefficients. For the MSE part (similar to Phase 1 in Algorithm 1). we can reuse the proof in Theorem 2.1, and we replace the analysis for Phase in Algorithm with LASSO results for the coefficients in (Wainwright, 2019).
Our result also gives us the opportunity to compare the problems of learning directed and undirected graphs in rel- atively equal footing. In (Misra et al., 2020), they propose an algorithm that uses (d/ (d/3min) log n samples to find the undirected graph, where Suni is the minimum normalized edge strength. Thus, we can view 7 as the analog of Bmin for directed graphs.
Regarding the efficient algorithm, the most relevant prior work to compare with is (Ghoshal & Honorio, 2018). They perform L Li-regularized Gaussian MLE and their analysis assumes that the inverse covariance matrix of the variables is upper bounded in L norm. This assumption is qualitatively different from the ones we make, which makes the two results incomparable.
We now provide proof sketch for Theorem 2.1 (for Algo- rithm 1) and Theorem 2.6 (for Algorithm 2), which can be analyzed using the same basic principles. Suppose we have a variable X, and a subset C [m] \ {i) with //| < d. Since our algorithms operate inductively, we will always assume that all nodes in come before in the topological ordering. The basic observation comes from analyzing the regression
In addition, Algorithm 1 needs knowledge of d, bmin to perform. We give modification of the algorithm when o², d or bmin are unknown to the learner, and we refer the reader to Appendix E.3, where similar guarantee (for Phase 2) is obtained with similar order of sample complexity.
m
(x((1) 33 xxx1))² 2
B = argmin
m jml
) (xx(0) - -BT xx(1)) 21
8 = argmin F. B
We also note that Algorithm 1 is not directly comparable
7

---

Learning Gaussian DAG Models without Condition Number Bounds
In particular, we care about the empirical and population Mean Squared Error (MSE) of this regression, which is
- i, then for every / C T. we have Var(X;|Xx) > Var(bjic + (x(x,) = b2 + 1 2 1 b where we used the fact that j comes after I in the ordering. Thus, we can see that the variance difference between the two cases is at least be.... The empirical variance behaves like a chi- squared distribution, which is subexponential. Thus, in order to distinguish between these two cases, we need (bmin) samples. This how the first bound of Theorem 2.1 arises.
MS X1) = m j=1
m
BT -x x x
-3
MSE(X,, = EE (80) x
Since all the variables are Gaussian, we can also write
3.2. Phase 2: Finding the parents
MSE(X,.X)) = Var
After we find the ordering, we turn to locate the parents for each node Let B, be the nodes before i in the or- dering. Suppose subset 7 C B. contains all the parents of i. Then Var(X;|Xx)
In order to estimate the conditional variance Var (X)( previous analyses (CChen et al., 2019; Gao et al., 2022)) use the well known identity
1. Now suppose there exists j E pa(i) with j & I. As we observed earlier, in that case Var(X;|XX, 2 1+ Var(b);,X(IX)) = 1 + 6j, Var(X,IIXy). The situation now is different because could also con- tain some descendants of j. In general, we establish the following Lemma.
Var (X, Di Ei M12
where we use the indices to denote subsets of the rows and columns. They proceed by estimating each of these subblocks of the covariance matrix, which requires condi- tion number assumptions. Instead, we use the observation that this conditional variance is equal to the MSE of the regression of X, with X₁. Using simple properties of our model, we can establish that M MS,,X, Xr) will be close to M SS Xy) when using O(dlog n) samples, with no dependence on the condition number.
Lemma 3.1. For any i E [n] and any I S [m] \ (t). we have Var(X;|X,) 2 1 + b2 /
Thus, T arises exactly as a lower bound in the conditional variance of a node if we also condition on some of its descen- dants. From Lemma 3.1 it becomes clear that the variance gap is now lower bounded by / Thus, by imple- menting this strategy of finding the parents, we would need O /b ) samples to distinguish the true parent set, but it does not give the optimal dependence on or bmin-
Continuing with the proof, we distinguish between two possible cases:
If pa(i) I then clearly Var (X,
= Var (Xa(n) = Var(e,) = 1, by conditional in- dependence of X, from all previous nodes, given the parents.
To obtain the optimal dependence, we instead use a strategy that was employed in (Misra al., 2020) to learn undirected GGMs. Specifically, instead of using the MSE to distinguish the correct parent set, use the coefficients computed by the regression. The idea is that if subset 7 contains all the parents, then if we regress X; with XIuJ for any other J C B,, the coefficients of the regression vector in J will be small. On the other hand, if there exists j E pa(i) with j E I, then there exists set with j E I, such that I TUJ includes all parents. Thus, if we regress X, with X the coefficient for X3 will be large. So, this becomes our distinguishing criterion for the correct neighborhood. It turns out that if we have O(T/b2, samples, the accuracy with which we can compute the coefficients in the regression is sufficient to make that distinction. Hence, we obtain the optimal sample complexity.
If pa(t) \ # 0. meaning that some parents are omitted from 1, then Var (X,|Xpa(s)) = Var(G) +
Var E Ejepu(o) bycX,1x,) 1.
Thus, our ability to distinguish between these two cases depends on how big the variance gap between them is. A smaller gap would make the task harder. We now explain why the tasks of finding the ordering and finding the parents have different variance gaps, which will result in different factors appearing in the sample complexity of each task.
3.1. Phase 1: Finding the topological ordering
4. Simulation Results
Algorithm 1 works inductively. Suppose we have found the first nodes in the ordering, which is the subset T,. For every i & T,, i should be the next node in the ordering only all of its parents lie in T,.. If it does, then there exists a subset I C T, with [1] = d such that Var(X, = On the other hand, if there exists j & Tr such that
To validate our findings, we run Algorithm 1, Algorithm 2, and the algorithm in (Gao et al., 2022) on synthetic datasets generated by DAGs) and compare their performance. The samples are generated i..d. from a randomly generated DAG G, which constructed as follows: (1) Draw random
8

---

10, 15, 20, 25, Learning Gaussian DAG Models without Condition Number Bounds
permutation of {0, n 1} as (0)), ............................... a(n 1). (2) For each node o(j), choose a (O),...., o(j 1) as parent of ((j) i.i.d. with probability min(
and (Gao et al., 2020) have a trend of decreasing to zero, and Algorithms 1, 2 converge faster then (Gao et al., 2020) when the number of samples exceeds 600. More simulations and code, including comparisons with the PC algorithm, can be found in Supplementary Materials E.2, F, and G.
742 2) Then,
truncate the parent set to at most d nodes.
(3) For each edge i + J. we assign bi, i..dd uniformly sampled from [-bmax; bmin]] U [bming bmax] where 0 <
de4, n=10
de4, nw
buinn bmax are tunable parameters.
'22
6oo Aigrittmm Aggattim
Algoritim
We note that this is similar to the topology that was chosen in (Gao et al., 2022) to simulate their algorithm. After sampling the graph, we generate m independent samples according to the model specified by equation (17).
Algortim
G
10-
0.-
6.0
.000
Condition number
bujG)
Nax
17.5
n=22
1030
300
15.0
Algarthim
25-
330
Algarthm
Agaattim
20
5.0
2.5
0.0
Figure 1. Growth of condition number, and maximum variance as n grows, averaged over 1000 random graphs with same d = 4, bmin 0.5, bmax
1000
Figure 3. False positive results for three algorithms. We still use blue lines to represent the algorithm in (Gao et al.. 2022), the orange line is Alg. I, and the green line is Alg. 2. The x-axis represents the number of samples m. and the y-axis represents the average false positive edges. Also, the shade represents the 95% confidence region (centered with mean and 22r error region) of each case.
d-4, n-10
d=4, n5
3.8
0.-
3.8
3.7
The first thing to notice is that our algorithms seem to be less accurate than (Gao et al., 2022) for a small number of samples (up to around 400 samples). On the other hand, our methods clearly outperform (Gao et al., 2022) once the number of samples becomes large enough, with the latter even failing to be consistent. One possible explanation for this behavior is that the sample complexity of Algorithms I and 2 have a worse constant than (Gao et al., 2022). For Algorithm 1, in order to find the parents, we iterate over all d-sized subsets both to find the candidate neighborhood and to test it with respect to all other neighborhoods, while (Gao et al., 2022) simply choose the subset with the smallest conditional variance, which could result in better constant. But if the number of samples becomes even slightly larger than 400, especially for larger n, the effect of the condi- tion number begins to show, as (Gao et al., 2022) fails to converge to the truth. Indeed, in Figure 1, we have plotted the growth of the average condition number K, TGG), and maximum variance R as n increases. is clear that k grows superlinearly, 7 grows sublinearly, and R grows linearly. This could explain the deteriorating performance of (Gao et al., 2022) when x increases, compared to our methods, which remain consistent.
Goo 22
Gas et '22
3.6
Algenthin
Alprrtiin Algorthm2 2
Algenthin
1000
850
1.0
3.0
0/8
2.0
0.-
Gan 72
3.0
202
800
1000
200
600
800
1002
Figure 2. Accuracy results for three algorithms. The blue lines represent the algorithm in (Gao et al., 2022), the orange linc is Alg. 1, and the green line is Alg. 2. The x--aiis represents the number of samples m. and the y-axis represents the accuracy. Figure 2 shows the results of
the three Algorithms with an increasing number of nodes
degree d = 4, bmin = 0.5, Imax = 1, and an increasing number of samples m = 200, 400, 600, 800, 1000. For each combination of parameters n, d, mn the accuracy is tested by generating 45 random graphs, running all three algo- rithms them, and calculating the fraction of successful recoveries of the graph.
We also present a plot (Figure 3) of the average number of false positive (reported but not in the ground truth) edges for the graph with the 95% confidence intervals. We observe that the number of false positive edges in Algorithms 1, 2,
The experiments were run using Python with numpy and scikit -leann package on an 11th Gen Intel Core i7- 11800H 2.30 GHz CPU with 16GB of memory.

---

Learning Gaussian DAG Models without Condition Number Bounds
Impact Statement
Duchi, J. Lecture notes for statistics 31 //lleccrraaa en- gineering 377. URL: https://stanford. edulciass/s-
This paper presents work whose goal to advance the field of Machine Learning. Specifically, this paper presents new algorithms for learning a specific type of Bayesian Network, graphical model that occurs in many disciplines. Although these algorithms are novel, this paper contains theoretical results. We have discussed the related work and open direc- tions, which mainly focus on the theoretical aspects of the problem. The data generated for the empirical experiments are synthetic and do not have any privacy issues. We do not exclude the possibility that there are potential societal consequences of our work, none of which feel must be specifically highlighted here.
tats311/Lecturs//full notes. pdf. Last visited on, 2:23, 2016.
Duncan, O. D. Introduction to structural equation models. Elsevier, 2014.
Friedman, J., Hastie, T., and Tibshirani, R. Sparse inverse covariance estimation with the graphical lasso. Biostatis- tics, 9(3):432-441 2008.
Gao, M. and Aragam, B. Efficient bayesian network struc- ture learning via local markov boundary search. Ad- vances in Neural Information Processing Systems, 34: 4301-4313, 2021.
References
Gao, M., Ding, Y., and Aragam, B. A polynomial-imme algorithm for learning nonparametrie causal graphs. Ad- vances in Neural Information Processing Systems, 33: 11599-11611, 2020.
Anandkumar, A., Tan, V. Y., Huang, F., and Willsky, A. S. High-dimensional gaussian graphical model selection: Walk summability and local separation criterion. Journal of Machine Learning Research, 13(Auuust), 2012.
Gao, M., Tai, W. M., and Aragam, B. Optimal estimation gaussian dag models. In International Conference on Artificial Intelligence and Statistics, pp. 8738-8757. PMLR, 2022.
Ao, R. and Chang, L. Influencing mechanism of regional ageing china based on the structural equation model. Acta Geographica Sinica, 1570811157-1158, 2020.
Bresler, G. Efficiently learning ising models on arbitrary graphs. In Proceedings of the sorty-seventh annual ACM symposium on Theory of computing, pp. 771 2015.
Gao, M., Tai, W. M., and Aragam, B. Optimal neigh- bourhood selection in structural equation models. arXiv preprint arXiv:2306.0222
2023.
Byrnes, J. E., Reed, D. C., Cardinale, B. J., Cavanaugh, K. C.. Holbrook, S. J., and Schmitt, R. J. Climate-driven increases in storm frequency simplify kelp forest food webs. Global Change Biology, 1713)251332 2011.
Ghoshal, A. and Honorio, Information -hoooation-thhoreticc limits of bayesian network structure learning. In Artificial Intelli- gence and Statistics, pp. 767-775. PMLR, 2017a.
Ghoshal, A. and Honorio, Learning identifiable gaussian bayesian networks in polynomial time and sample com- plexity. Advances in Neural Information Processing Systems, 30, 2017b.
Cai, T., Liu, W., and Luo, A constrained minimization approach sparse precision matrix estimation. Journal of the American Statistical Association, 106(494):594- 607, 2011.
Ghoshal, A. and Honorio, J. Learning linear structural equa- tion models in polynomial time and sample complexity, In International Conference Artificial Intelligence and Statistics, pp. 1466-1477 PMLR, 2016 88
Cai, T. T., Liu, W. and Zhou, H. H. Estimating sparse pre- cision matrix: Optimal rates of convergence and adaptive estimation. 2016.
Chen, W., Drton, M., and Wang, S. On causal discovery with an equal-variance assumption. Biometrika, 106(4); 973-980, 2019.
Grace, J. B., Michael Anderson, T., Smith, M. D., Seabloom, E., Andelman, S. J., Meche, G., Weiher, E., Allain, L. K., Jutila, H., Sankaran, M., et al. Does species diversity limit productivity in natural grassland communities? Ecology letters, 10(8):680-689) 2007.
Chickering, D. M. Learning bayesian networks np- complete. Learning from data: Artificial intelligence and statistics V. pp. 121-130, 1996.
Heckerman, D... Geiger, D., and Chickering, D. M. Learning bayesian networks: The combination of knowledge and statistical data. Machine learning. 20:197-243. 1995.
Chickering, D. M. Learning equivalence classes of bayesian- network structures. The Journal of Machine Learning Research, 2:445-499 4488, 2002a.
Johnson, C., Jalali, A., and Ravikumar, P. High- ligh-dimmnssonal sparse inverse covariance estimation using greedy meth- ods. In Artificial Intelligence and Statistics, pp. 574-582. PMLR, 2012.
Chickering, D. M. Optimal structure identification with greedy search. Journal of machine learning research, 3 (N0v):507-554, 2002b.
10

---

Learning Gaussian DAG Models without Condition Number Bounds
Kalisch, M. and Bühlman, Estimating high-dimensional directed acyclic graphs with the pe-algorithm. Journal of Machine Learning Research, 8(3), 2007.
Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A., and Jordan, M. A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7 (10), 2006.
Kelner, J., Koehler, F., Meka, R., and Moitra, A. Learning some popular gaussian graphical models without condi- tion number bounds. Advances in Neural Information Processing Systems, 33:10986-10999 10998, 2020.
Spiegelhalter, D.J., Dawid, P., Lauritzen, L., and Cow- ell. R. G. Bayesian analysis in expert systems. Statistical science, pp- 219-247, 1993.
Spirtes, P., Glymour, C.,, and Scheines, R. Causation, pre- diction, and search. MIT press, 2001.
Kelner, J., Kochler, F., Meka, R., and Rohatgi, D. Lasso with latents: Efficient estimation, covariate rescaling, and computatioonllssatisiiiaaa gaps. arXiv preprint arXiv:2402.154002 2024.
Tan, V. Y., Anandkumar, A., and Willsky, A. S. Learn- ing gaussian tree models: Analysis of error exponents and extremal structures. IEEE Transactions on Signal Processing, 58(5))7001-2711 2010.
Klivans, A. and Meka, R. Learning graphical models using multiplicative weights. In 2017 IEEE 58th Annual Sym- posium on Foundations of Computer Science (FOCS), pp. 343-354. IEEE, 2017.
Tan, V. Y., Anandkumar, A., and Willsky, A. S. Learning high-dimensional markov forest distributions: Analysis of error rates. The Journal of Machine Learning Re- search, 12:1617-1653, 2011.
Lauritzen, S. L. Graphical models, volume 17, Clarendon Press, 1996.
Uhler, C., Raskutti, G., Bühlmann, P., and Yu, B. Geometry of the faithfulness assumption in causal inference. The Annals of Statistics, pp. 436-463. 2013.
Lokhov, A. Y., Vuffray, M., Misra, S., and Chertkov, M. Optimal structure and parameter learning of ising models. Science advances, 700791, 2018.
van de Geer, S. and Lederer, J. The bernstein-orlicz norm and deviation incqualities. Probability theory and related fields, 157(1):225 250, 2013.
Malioutov, D. M., Johnson, J. K., and Willsky, A. S. Walk- sums and belief propagation in gaussian graphical models. The Journal of Machine Learning Research, 7:2031- 2064, 2006.
Vuffray, M., Misra, S., Lokhov, A., and Chertkov, M. Inter- action screening: Efficient and sample-optimal learning ising models. Advances in neural information pro- cessing systems, 29, 2016.
Meinshausen, N. and Bühlmann, P. High-dimensional graphs and variable selection with the lasso. 2006.
Vuffray, M., Misra, S., and Lokhov, A. Efficient learning of discrete graphical models. Advances in Neural Informa-
Misra, S., Vuffray, M., and Lokhov, A. Y. Information theo- retic optimal learning of gaussian graphical models. Conference on Learning Theory, pp. 2888-2909. PMLR, 2020.
tion Processing Systems, 33:13575-13358 2020.
Wainwright, M. J. High-dimensional statisties: A non- asymptotic viewpoint, volume 48. Cambridge university press, 2019.
Pearl, et al. Models, reasoning and inference. Cambridge, UK: CambridgeUnivvertttyeess 19(2):3, 2000.
Wainwright, M. J., Jordan, M. I., et al. Graphical models, exponential families, and variational inference. Founda- tions and Trends@ Machine Learning. 1(1-2):1-305. 2008.
Peters, J. and Bühlmann, Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1):219 228, 2014.
Peters, J., Mooij. J. M., Janzing, D., and Schölkopf, B. Causal discovery with continuous additive noise models. 2014.
Wang. W., Wainwright, M. J., and Ramchandran, K. Information ceemmation-theortiic bounds on model selection for gaus- sian markov random fields. In 2010 IEEE International Symposium on Information Theory, pp. 1373-1377. IEEE, 2010.
Robins, J. M., Scheines, R., Spirtes, P., and Wasserman, L. Uniform consistency in causal inference. Biometrika, (3):491-515,, 2003.
Wu, S., Sanghavi, S., and Dimakis, A. G. Sparse logistic regression learns discrete pairwise graphical models. Advances in Neural Information Processing Systems, 32, 2019.
Rudelson, M. and Vershynin, R. Sanson-wright inequality and sub-gaussian concentration. 2013.

---

Learning Gaussian DAG Models without Condition Number Bounds
Yu, B. Assouad, fano, and cam. In Festschrift for Lucien Le Cam: research papers in probability and statistics, pp. 423-435, Springer, 1997.
Yuan, M. and Lin, Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19-35, 2007.
Zhang, J. and Spirtes, P. L. Strong faithfulness and uni- form consistency in causal inference. arXiv preprint arXiv:1212.2506, 2012.
Zheng, Y., Huang, B., Chen, W. Ramsey, J., Gong, M., Cai, R., Shimizu, S., Spirtes, P., and Zhang, K. Causal- learn: Causal discovery in python. Journal of Machine Learning Research, 25(60):1-8, 2024.
12

---

Learning Gaussian DAG Models without Condition Number Bounds
A. Proof of Theorem 2.1
The Algorithm consists of two parts: finding the topological order and finding the parents. We use the notation (i ) denote that comes (strictly) before j in the topologiaal ordering. ž and > are defined similarly.
We will sometimes the notation b when want apply some Lemma that has variable b with the value a. e.g. in the proof of Theorem 2.1, we use Lemma A.3, we are going to take 8 in the Lemma A.3 as (2N), where this d is the ò in Theorem 2.1
As we have explained in Section 3, the general approach in both phases of the algorithm is to use the Mean Squared Error (MSE) of linear regression to estimate the conditional variance of a node given some subset of nodes that appear before it in the ordering. The variance gap between correct and incorrect subsets of nodes will be different each phase, which gives rise different sample complexities for finding the ordering and finding the parents. We give the details for each phase in the above discussion.
below. Without loss of generality, 2 2, bmin I /3
A.1. Finding the ordering
Algorithm I builds the topological ordering iteratively by finding at each step the next node in the topological order. Suppose inductively that the first t nodes form a valid topological ordering, for some t < n and let T be this subset of nodes. We will show that Algorithm correctly identifies the + 1-th node in the ordering. To achieve that, it considers all other nodes i [n]] \ T' and calculates the conditional variance of given subsets J _ T with [J] S d. This is the step where previous analyses ((Chen et al., 2019; Gao al., 2022) introduce condition number assumptions, in order estimate this conditional variance. Instead, will use the MSE of linear regression to estimate this conditional variance. Indeed, can show that the empirical and true MSE are close to each other using number of samples that is independent of the condition number.
The Empirical mean squared error consists of two parts. We can write it as
MSE- - Σ(x
-x B,)²
3=1
1 BJ m X
m -Xx() Br xx(18))
variance error
beta error
In the above, we have identified the two key quantities the error, which we call the Emprirical Variance Error (call it V E) and Empirical Beta Error (call it BE) accordingly. By Schwartz inequality, we have
VE + BE 2VVE BE MSE VE + BE + 2VVE BE
In order to decide whether i can be the next node in the ordering, we have to determine whether all parents of are contained in T' or not. The following Lemma shows that, in the population case (infinite samples) there is way to distinguish between these two cases by looking at the conditional variance.
Proposition A.1. (1) T contains all the parents, then there exists some J C T and [J] S d such that the conditional variance of X1|æs is 1.
IfT does not contain all the parents, then there exists a parent that is not in T, so for any coefficient vector By each X, Br is a Gaussian random variable that has variance at least 1 + 12 bii Specifically, for any J C T, for any coefficient Br, X; X333 also has variance > + -2
Proof. This property is also proved in (Gao et al., 2022), but for completeness we provide a proof. We first prove Property A.1. The first point follows immediately from the definition. For the second point, let k be the last parent of in the topological order. By the preceding discussion, it follows that & T. Let Ti be the set of nodes that is prior to k. So
13

---

Learning Gaussian DAG Models without Condition Number Bounds
T C T3 Then, have
Var( X; - Fr = -ExT, [Var(X; -
Var(B(X, - Xr
(3)
Bxxx Var(XX XTxT//XXTT = Exr, [Var(X;
(4)
Exr, Ex,
Ex,, [Var[E, (X,
(5)
= - + Exr, (Var(be X = 1 1 b ++ b
(6)
We explain each line: (4) is by the law of total variance, (5) because the variance is greater then zero and XT is constant under the condition of Xr, because T S Tk. (6) is again by law of total variance, (7) is because while conditioning on Xr, and Xg. all parents of are conditioned. Thus the variance of X; is because of the model. Also, the expected value of X; condition on XT, and X is X; = Ej bj = Ei d bi-iX; ; + b Since conditioning on XT. that all prior nodes than k the value Ejfu X is a constant, so the conditional variance of (X))) is the same as that of bR. which by the definition of the model is b/_ times 1, which is 6i
The preceding property holds in the population setting. We would like to establish similar separation using finite samples. This is done using the following Lemmas. The first one quantifies how large BE is when the set might might not contain all parents of node i. To state it, it will be helpful to we extend the meaning of the ground truth coefficient: since
Ens
Eji
the joint distribution Xi) is a multivariate Gaussian distribution with covariance matrix write
we can always
EN
Ens
E =
where 3 = En. We call this vector 3 the population coefficients of X; with respect to X,, The reason is that if we performed linear regression of X; on X, with infinite samples, the result would be 8ÿ and the conditional variance is En Σe>DoドE J contains all of i's parents and none its descendants, then we have the population coefficiens of X; wrt X, are exaclty the coefficients bj-ts-
Lemma A.2. (For beta error.) In Algorithm 1, let S [m]] Vi) with |J|= k S d and Var(X,|Xy) = V2, B* be the populations coefficients of X; w.r.t X, and let 33 be the result of linear regression of X; on X, using m independent samples (1), ,..,, ((m) drawn from this distribution. Then, there is a universal constant CB = 4, such that for any E< </1/ 1/2, if we are given m C3(Olg(11/))) (16))// samples, then
m
-x²) m
<VV =
with at least 1 probability.
Thus, this shows that BE will be relatively small if J contains all parents. The next Lemma will be used to argue that VE will also be small in that case. It's proof is standard and uses the subexponential property for the x² distribution.
Lemma A.3. (For variance error.) Let Y1, Y2..... ... Y. be i.i.d. Gaussian random variables with variance o². Then, there is universal constant Cv 6 such that for any 1/2,6 8 < 1, if we are given 2 Cv (log(2//6)/ee) sumples, then with at least - å probability, we have
- 8) <
y2 < (1+00) + m
m
The proof of the Lemma A.2 is presented in Section D.3 and the proof for Lemma A.3 is in Section D.4. We are now ready to prove the first assertion of Theorem 2.1, which says that Phase of Algorithm correctly identifies the topology with high probability. Aided with these two lemma, we return the proof of Theorem 2.1
Proof. Let's define N = (() + (1) +... (2)), which is an upper bound of the number of linear regressions performed in Phase of Algorithm 1. We would like to argue that all these regressions succeed with high probability, so we employ
14

---

Learning Gaussian DAG Models without Condition Number Bounds
union bound. This involves bounding log N, which we now do. Since have >> we have N 2n(2) and thus by Stirling Formula,
log < 1og (d 1) log dlog(d/e) = + 1) + + logd < 5dlog(n/d).
Now let's assume inductively that Algorithm has identified the first t nodes in the ordering correclty and let T be this subset of nodes. The Algorithm then examines every i [n] T determine if it is the next in the ordering. There are two possible cases.
T contains all the parents of i. 7 contains all the parents of i, then can be used as the / + 1-th node in the ordering. We prove that Algorithm will indeed select it. By definition, there exists a J S ' with |J| S. d that contains all the parents. We want upper bound the MSE of regressing X; on x,, To apply A.3, we notice that x x B3 are i.i.d. (0 1) random variables.
Suppose it holds that
VE - 1 1 / BE < 64 min//44 /64.
Then, wherever bmin < 2 (which can wlog assume the empirical MSE can be bounded
MSE S VE + BE + 2VVE. BE S + b2 //iin/ /8 + miia/ /64 + 2/ 12 /iiin// /8- /iin/ < 1 1 62 min!
(8)
Thus, the condition in the If statement of line 9 of Algorithm I would be satisfied, which means i would be selected the next node. We now establish that (7) holds with high probability. We first notice that Var Var(X;|XG)) = Var(e, = 1. Thus, we can apply Lemma A.2 by plugging 5 t- 6/(2N), random variables by definition of the model. Thus, By union bound, these two lemmas give that with probability at least 5/N (7) holds for all regressions in Phase 1, as long the number of samples satisfies
mmin/ We also notice that x x 33 are i... N(0, 1) apply Lemma A.3 by plugging in ō t 6/ E b2 /
m max (CD(k (2 /64, Cv l4
Therefore, if we choose C 5120 2 max(1280Cp, 6400) bound on m.
C(d Codlog((/d)) (ogllogg///)++o((()))))) will imply the previous
T does not contain all parents of We want to establish that will not be selected as the next node. Let J C T be any subset of T' that is considered by the algorithm for i. By Property A.I we have
Var(X:|XT) 2 1 1 02 1
Thus, the law of total variance yields
Var(X,|Xs) y [ar(XX +
2 Ex [Vr(( 2 + 62 min-
Letting 12 Var(X,]X,), our goal will be to lower bound the empirical MSE. Suppose the following holds
VE 2 (1 - b /8)v2
BE S (bminn/ /64)//
(9)
then this would imply (the third inequality is due to the monotonicity of VE.)
MSE 2 VE + BE - 2 BE > VE - 2VVE BE 2 1 b / - 2/ 12.in// / /8
8
Thus, since V2 > 1 62 - this would imply for bmin
v3 (again, we can w.l.o.g. assume this),
3 62 V2 V2
MSE 2 (1
1 / /2.
That would mean that the condition in the If statement of line 9 of Algorithm 1 is not satisfied, which means i will not be selected as the next node. Establishing that (9) holds with probability at least 1 - 8 then involves a similar application of LemmasA.2 and A.3 to the one we saw in the previous case, which we omit.
15

---

Learning Gaussian DAG Models without Condition Number Bounds
A.2. Finding the parents
The next phase focuses on finding the parents of each node i, assuming we know the subset I of nodes that come before it in the topological ordering. Here, we can make connection between the linear regression problem for the directed graphs and the one for undirected ones. In particular, we focus the marginal distribution of (X,, Xr). which is itself coming from a Bayesnet, which is obtained from the original Bayesnet by keeping all nodes up to i in the ordering.
We start by proving some preliminary properties of the inverse covariance of a Bayesnet model in the following lemma:
Lemma A.4. Let (X1; (X1,X2..... Xx) be a multivariate Gaussian, with covariance matrix E. and inverse covariance matrix Σ-1, distributed according to model 17. Then we have: for i + j.
Oij = bi-s + مز + Σ bi bi
for any i,
i =1 Σb D3-o3-
Proof: We can show it by calculating the likelihood of X = (X1, X2,...... Xn): it is exp(- because of the joint Gaussianity, Let
(n) be the topological order. By the Bayes rule, the likelihood of X = (X1, X X,)
is also
p(X1, X2...., X,) = II 1=1
pl X (x(())
X Xo(I)...x((a))))
exp
(Xo())
(a(i-)) (o))
exp
X (1)
ba ba(-1)) Xo (o)
2 i=1
/=1
exp
(X₃ -
Σ
2
im1
j-ti
So, we have
XToX = Σ(x, - Eomme
1=1
J-1
Therefore, Oạ is the coefficient of X2, in the above sum. We get a contribution of from the term (X, - (X,-Ex, by bxx)
and contribution of the sum of رنا from the sum of (X, Ex-, bn->j X,) 2 i.e. O.; =1 E, b²z
Similarly, O. is the coefficient of XX,X,. We have a contribution of by if i is a parent of j. coming from the sum of (X, - Ex_ bx-rxx,), or by- if j is a parent of i, coming from the sum of (X; Ex bx, Also, for each common child of j, we get a contribution bi-r1bjstt in the sum of (Xi - Σm 11) Overall, Oi b, + b; + E, -ib
Now we combine these observations with a lemma in (Misra et al., 2020), which is stated for arbitrary undirected Gaussian Graphical Models. It quantifies the accuracy of estimating the coefficients in OLS for a node, if we regress it with a small subset of nodes that contains its Markov Blanket. We will apply to the marginal distribution of (X,, XT).
Lemma A.S. (For general Lo constrained sparse linear regression) (Proposition 2 in (Misra et al., 2020)) Suppose we have a multivariate n-dimensional zero-mean Gaussian model with inverse covariance matrix O and m i.i.d. samples from this model, where
8 _2 log(()) + log( 2d ) 8
m 2 2d +
g2
16

---

Learning Gaussian DAG Models without Condition Number Bounds
For node i, let (i) be the neighbors of in the undirected graphical model induced by Q Suppose A S [n]] \ {i) with |A] I 2d and T()) S A. Then the OLS coefficients B when we regress X, on XA using m samples satisfy with probability at least 5 the following bound: Vj A,
O OH
J3 Os
Big
I
1
We now make the connection with our DAG model more precise. As discussed earlier, let i be some node and the set of the node before in the ordering T. We consider the marginal distribution of (XT, X:). which still Bayesnet. Let O be the inverse covariance matrix of (X),-X... Since is the last node in the ordering for this smaller Bayesnet, it is easy to that its Markov Blanket T (t), if the joint distribution is viewed as an undirected Graphical Model induced by θ is exactly equal to the set of parents pa(i) in the DAG. This can be inferred for example from Lemma A.4 by noticing that Oij is non-zero if and only by: * 0 (since has no children in this marginal distribution).
This allows us to apply Lemma A.5 to find the parents of and gives rise to the following Corollary.
Corollary A.6. Let's i be a node and T be the subset of nodes that come before in the ordering. Suppose we use m i.id samples to regress X; on X,, where J C T, |J| S d and pa(i) C J. Let B be the coefficients output by the regression. Suppose
m 2 64dlog(n) T bm2
Then, with probability at least 8 we have that for all j E pa(i).
13 > bmin/2
and for all j & pa(i))
15 bmin/4
Proof. By the preceding observations, we can apply Lemma A.5 to analyze the regression. If θ is the inverse covariance matrix of the joint distribution of XT), then by Lemma A.4, since X; is the last node, we have Oij is bj O =
and ,₁ = + Σ, j b7
Thus 10 0 = 1bul 2 min j E pa(i). Also, jf $ pa(t), then /0.//aal = 0. Thus, by applying Lemma A.5, we need to ensure E v 1 + Ej 1/ S /minn/2 2 for all í,j. So, it suffices to have E < Smin/ (2/7).222..)))) By union bound, we need, 8 +- ó/n. To determine the sample complexity, we plug these values Lemma A.5
8
4
m 2d + (bmin/ (bmin/2v7F) d dlog + (bmin/2 72/2
log( log(2dh) (2dn).
Where we take m (odlog(n) T min suffices, and this finishes the proof.
Below we give the details of the proof of the second half of Theorem
Proof. We first prove that for all and prior nodes T in the topological ordering, Algorithm correctly identifies a subset that contains pa(i). The strategy is similar to (Misra et al., 2020). For every candidate J C T with [T] d, iterate over all other K C T with K| = d and regress X; on X XJuk.
Suppose pa(i) C J. Then, for all other K, for each j K we have 16 < bmin/2 2 with prob > 1 by Lemma A.6, so the condition in line 17 of Algorithm 1 doesn't 't hold and J is accepted as true superset of the parents.
Suppose there is j E pa(i) with $ J. Then, there exists K with j E K and pa(i) C J U K. Then, regressing with J K will yield |3 bmin/ by Lemma A.6, leading to the condition in line 17 to be satisfied. Thus, set J is rejected.
After we find a superset J of the parents, we can simply remove the nodes with small 3 coefficients in the regression and appeal to Lemma A.6 again for correctness in an identical fashion.
17

---

Learning Gaussian DAG Models without Condition Number Bounds
B. Proof of Lower Bound (Theorem 2.3)
The standard testing setting involves a family distributions indexed by a set V. A random variable drawn uniformly at random from V and X is generated according to the distribution corresponding V. Note that X could be one multiple samples from this distribution. Then, an algorithm that takes as input X and guesses V is denoted by V. We use Corollary 9.4.2 from (Duchi, 2016), which is version of Fano's s Inequality for testing.
Lemma B.1. Assume that is uniform on V. For any Markov chain V X
P(V = V) 2 X(V:X) + log
log IV
In our case, the V corresponds to a of graphs G { ......, . Gk} that we would like to distinguish from each other. V denotes graph chosen uniformly at random among V. Let Qv denote the distribution of one sample drawn from graph V and Q denotes the product distribution m i..d. samples drawn from V. The mapping X + V is the estimator of the graph from the samples for model. From Equation (9.4.5) in (Duchi, 2016), we have
X(V
Σ V.VVEV
KL(QQ"IIQN)
Σ KL(Qv]IQv.)
DP
2
G.G'ES
(This is due to the tensorization of KI. for product distributions.) Therefore, we can construct a family of [g] graphs such that for any G, G' E G, have KL(Qa] Qa+) a, then need at least 2 a prediction probability 2 1/2. We consider the following ensembles:
lug(18))22-oog samples to obtain
Ensemble for log(n)/<m min- Assume 2 divides n, consider g = {Go, GI, G1, "a....,Grpph As Figure 4 shows, all graphs are defined on random variables Y1, Y2.... Yn/2 and Zi,..., Zn/2. Go is a direct matching of Y; - Zi. with each edge having weight byinn G3 is kth model and is equal to Go. except that the edge between Y& and Zg is reversed.
2
Y
Zi
Y2
Zz
Y2
Z
Y/3
Zz Y3"
Model Go
Model G
Figure Schematic graph for the ensemble. Every G, for j 2 is obtained from Go by reversing one edge the matching of model Go. e.g. G1 is just reversing Yi Zi in Go.
We notice that this model is a product measure of n//2 disjoint pairs of variables. Therefore, any two models differ in at most two edges. Let (Y, the joint distribution of a two variable DAG with edge Y + with weight bmin- Then, by the preceding observations, for any two models G3 Gy, their KL divergence is no more then Z)(( Next, calculate this distance.
First, the distribution of (Y1, Z1) has the multivariate Gaussian distribution of covariance matrix
1
Umii 62 62
Σ
=
biin
+
Similarly, (Z,Y) follows multivariate Gaussian distribution of covariance matrix
Σg
(
112 bmin
min
bmin
18

---

Learning Gaussian DAG Models without Condition Number Bounds
Therefore, the KL distance is
Σ Σ
1 (-2 TF 2+ Tr(
(1 11 82
b4 min
b3 min 6
2))) =
log
-2
ming -ba
b min
1
min
So, in this case, any two models has KI. distance at most iiin and this implies that with probability 2 1/2 we need at least (log(n/2)/2 log 2)/bmnn log(n/8)/ (28m)/)))))mnn))) samples to find the right model.
Ensemble for log(n)T/Bb2) Assume 3 divides n, consider G = [G,, G1, G1, G2....,, Gn/3)- Our graphs now consist of independent blocks on variables X,, Y, Z. As Figure shows, Go the topology of each triplet is X; - X; + Zi Y + Z, The weights are: X; + Zj with bmin weight, X; + Y; with B weight, and Y; + Z; with by weight for some B> bmin 0. For G, with > 1, the topology the same as Go but we change exactly one of the triangles Zz Bbain weight. it becomes X; Z; with no edge, X; Y; with B weight, and Y; Z₃ with bj + 1+B
X
Xi
bx
B
xx = b
bxt .
Bb
lng-2 2 bj +
mm
1 ++
B
Y1
Z
Y1
by 2 b
X2
X2
Y2
2
Y2"
Z2
X
Xg
Zz
Y3
Y3"
Model Go
Model G
Figure 5. Schematic graph for the ensemble. Every model G; with i 2 is obtained from Go by switching one triangle to a chain, ie. G1 just changing the triangle X, Zi to chain X; Yi Zi.
Let (X, Y, Z) be the distribution of the DAG with topology X + Z with bmin weight, X - Y with B weight, and Y Z with b1 weight. Let (X', Y', (2) be the distribution of the DAG with topology X' Y' with B weight, and Y' Z' I-B² weight. Similar to the observations for the previous ensemble, the KI. distance between any two models is
with bj + banin
bounded by 2 KL((X, Y, (XX', XX, Y', ZZ))..
The distribution of (X, Y, Z) has the multivariate Gaussian distribution of covariance matrix
B 1 1+ 1+B* B by + B2b + Bbmin
Bb + bmin
Σ
B
61 + B²bb + Bbmin 1+ (Bb, + bmin) ) + 83
Bb, + bmin
The distribution of (X', Y", Z²) is multivariate Gaussian distribution of covariance matrix
Bby + B2bam 1+B" 61 + B2b + 13bin
1
Vi =
B
1+ B²
Bbi
n B'bain 1+B*
Bºb Blmin
1 + (b: Bbmiy (B2 + +1)
Therefore, we can calculate
Er E52|
(because it is the inverse of product of two upper triangle matrix with
diagonal 1: Notice that
|
[([1 B) 1)
B3) -1 1((1 B) 1) 1(1 B) -1 1 1) and
+
banin
0
bmin buin B H? 1
6min
E Σn'Σ Σラ
=
bu
B+
1
:
Σng
1 Σ1 =
Banin Omn B+ + +
( (1+B4 B2)
buii. ( B+ 1322
(1+ Ba
.3 + N
by + By &
13
B+
-bain
1 132
1 +
14 13
62
19

---

Learning Gaussian DAG Models without Condition Number Bounds
6 Therefore,
Both of the traces are 3 + 1
Σ Σ2
(
62
XU Y', 2
log
+ Tr(::
= -3
= /0min(11 + B²)
1 B
So, in this case, any two models has KL distance at most m + B²) /inn TT. and this implies that with probability 1/2, we need at least (log(n/3)/2 log 2)7/b2 miin log(n/ 12)7 (2) samples to find the correct topology.
Thus, from both ensembles, Theorem 2.3 follows.
C. Analysis of Efficient Algorithm (Theorem 2.6)
give the full version for Algorithm below.
Algorithm 3 Efficient algorithm for recovering the topology
input Samples (X(1,, (2) X(m)
output Topology
1: T +- i: Var(X;) & (1 + P
2: while [ < n do
3:
MM 0
4:
forie E [n] do
5:
Choose x = O( 10Vooool log( n R/ /m)
6:
3 LASSO regression of X; on XT
7:
J [Al1B&kl 2 buin / if [] > d then go line 4.
8:
Perform OLS for X, on X,, get coefficient 8.
MSE < m E, ( x
9:
10:
if MSE 1 11 min /2 then
11:
M + U {i}
12:
pa(i) + [j E JJ 2 buin//)
13:
end if
14:
end for
T Taappend(M)
end while
In line 6, we perform the LASSO algorithm to find the next vertex for the topological order as well as to find the parents. For any node and subset that appears before in the ordering, we are optimizing the Lagrangian LASSO:
8 = argmin AERITI
Σ -3 + A||3||
(10)
2N
For vector 3 and subset S' of coordinates, we use B3 to denote the subvector corresponding to the coordinates in S. As usual, let 3r be the population coefficients of X; with respect to Xy. The following Theorem from (Wainwright, 2019) provides guarantees about the estimation error of the LASSO.
Lemma C.I. (Combination of Theorems 7.18, 7.19 in (Wainwright, 2019)) Consider the Lagrangian Lasso where X is generated i..d. from.N(O. 2). The independent noise for each sample is wi() = X - -3- x(I) Choose the regularization
parameter 2 2 maxkeT Σn wl() x(1) For any 6* E gd, with probability at least T m/5z any optimal solution
-m/33
θ satisfies the bound
3 B 3/ < 9216
1135-11 + 12800 22(E) log n
S
+ 128 R
11
(1))
R2
20

---

Learning Gaussian DAG Models without Condition Number Bounds
valid for any subset S with cardinality |S| the smallest eigenvalue 2. >
S
Here, (Z) is the maximum diagonal entry of E and
25600 p=(2)) log n
There is a twist in the algorithm after computing the LASSO. First of all, if the output vector has more than d "large" coefficients, this clearly means that ''' does not contain the full neighborhood and thus is not the next node in the ordering. If it has at most d "large" coefficients, instead of just calculating the MSE using the output of the LASSO, we run an additional OLS to verify whether i is the next node in the ordering. Then, will use the property A.1 to judge whether have the correct parents. Specifically, if is the next node after T, then LASSO will find the correct parent set J (not superset) and the expected MSE after the OLS is 1. If i is not the next node after T, then for any J C T, the expected MSE after the OLS is 1 + b Since the sample complexity for Theorem 2.1 is upper bounded by that of 2.6, this verifying process is going to be successful with high probability; we can see this by repeating the proof of Theorem 2.1. The reason we use OLS instead of Lasso is that it is harder to argue about the MSE of Lasso when not all parents are involved.
Therefore, the only task is to guarantee that LASSO will determine the true parents for any i that is next in the topological order of T. In order to have that, we need the C2 norm of the distance of 3 and 3 to be smaller than 2 /2. Applying Lemma C.1 above for LASSO in our algorithm, we need to select À and m appropriately for this to hold. Therefore, it suffices to choose X, m so that the following inequalities hold simultaneously:
R
m
2|51 2 d
(12)
25600 p² ( log 72
A2
12
9216
S
min
(13)
R2
16
X
12
128
W1 S
min
(14)
R
16
R2 log
2
12800
W <
min
(15)
RN
16
Equation (12) requires |S|] to be large enough, and Equations (13), (14) and (15) require || 3- smaller than //min/4 /4 which achieves the guarantee. Also, here introduce W1 be the maximum L1 norm of the coefficient vector correspond- ing to the parents of each node: that is,
Wi = max Σ 1bj-ss|-
Now, we bound the variables one by one. By our assumption, we have 2²( < R2, W, S dVF (because bi- S VT) and by Cauchy-Schwarte inequality, R 2 Ta+1)-= 2 1
2d-7. We apply the following lemma for all XT since the sub-Bat uub-Bayesnet has
the same model as the whole Bayesnet:
Lemma C.2. Assume we have the Bayesnet as Equation (17). Let R is the smallest eigenvalue of the covariance matrix of X1,X2, X,,..., X,, Then, we have R 2 (d+1)r'
Proof. We argue by Cauchy-Schwand inequality. Let Σ be the covariance matrix of X2, Xn, and θ = Σ'1 So the least eigenvalue of Σ is the largest eigenvalue of è. Let z = (I1, .nn) and ||>||| = 1, so that we have TT EL EL I; - bj ;;,; (by the proof of Lemma A.4.) By Cauchy-Schwartz inequality, we have El= (E) Fi by E
Ej= 21 E bj b), I (d 1)7. Thus, we have |[0|l.p (d 1)T. Since E = 0 we
have the least cigenvalue of Σ is no less than (a+1)=
Also, to bound X, we need the following Lemma, whose proof is given in Section D.5.
Lemma C.3. Let Z; be random variable of X Y; where X,, Y, are independent (O, 1). Then, there is universal
constant Cx = 32, such that for any < 1/2, 5 < 1/2 given m 2 Cx (x(log(//6))/2² we have that the following holds:
1
P
Σ
2
>E
<
21

---

Learning Gaussian DAG Models without Condition Number Bounds
First, calculate make all inequalities hold for |S| d that satisfy equation (12). This requires that m 51200 d² R2 log(n).
Also, since we have assumed that for each i, the variance of X, is E R, thus, for each k. ΣØ= w( x()xx has each w sampled independently from (0, 1) and x(") is sampled independently from (V Var( (ar((x)), where Var(X*) S Rº. Therefore, because we need to exccute n² LASSOs, so we need to substitute +- 6/n² in Lemma C.3, in order to apply a union bound. Therefore, by Lemma C.3 and union bound, we substitute ε < 2/(2R) and 6 t 8/4n² (the reason for the extra factor of is to accommodate some failure probability for the accuracy of r, VE, BE and for the failure of Lemma
/ /18 Ing(8n? 16)R
C.I). Thus, it suffices to have m 2 128 log((on /8)/>². Equivalently, this can be written as À <
mm
Solving (13), (14), (15) and substituting the bound for R, |S], W1, we have
128 log(8n" /6) R (Adr) -
/
02
9216
min
m
16
128 log (8m2 /
82
128
d² -1.5 _
min
16
R2 log n
82
12800
d3,"
min
16
Thus, we can choose m 2 22.2. log(n/8)) d ...73 RR. b-d to satisfy the bounds above. This finishes the proof for Theorem 2.6.
D. Additional Proofs
D.I. Proof for Lemma 2.2
By Lemma A.4, the diagonal entry of θ is
i 1+ Σ bE
So, by choosing x to be standard basis vector, get that there exist æ1, 22 such that
min O
wig Owa max O
is clear that
max O 2
by definition. Also, by choosing the diagonal element of the last node in the topological ordering, it is clear that
min Oal
Therefore, the condition number of θ is at least and the same holds for Σ.
D.2. Proof for Lemma 2.5
For (1) Let Σ be the covariance matrix, then consider all unit vector x, x Ea ranges from 1 (choose the first vertex in the topological order) to max; Var(X,) = R. Also, we know that the smallest and largest eigenvalues of Σ (Amin and Xmax) satisfy that Amin =in x Ex and Amax max |-||=1 x² Exr. So we have Amin < 1 and Amax 2 R and thus the condition number of Σ is at least R.
For (2) We consider a directed rooted tree with branching factor with variables Xi,...,,,),, Consider root X1 ~ (O, 1). And for any i, X, has two children X and X2++1 with recurrence X2: = Az; + E2F and X = Ax; + E2i+1 where E2s, E2s+1 are i.id. 1). Assume the tree has height h, and thus n 2 Choose some X E ( (122 1).
For any leaf (X, for 2*), we can write X, Es + ( 1 + X² + ... + 2h --
X2 -... X Therefore, its variance is
22

---

Learning Gaussian DAG Models without Condition Number Bounds
On the other hand, lower bound the condition number, need find different values of x: Ex for some | we choose r to be the unit vector with X1 coordinate to be 1, the variance is 1. We can also choose a to be 1/v24 for coordinates corresponding to leaves and otherwise. Then, x - Ex corresponds to the variance of the linear combination V2h (X 12)- We are going prove that the variance is
Var +.. X2h-1-1) = 1 (2X²) + (2A7) + +... (2
(16)
If this holds, the condition number is at least desired bound.
(2A7) (21--1((1-2)
By choosing À 2(a-1)/2 get the (21--1)(1-2)-
Now we prove (16). Suppose Ve is the conditional variance of 1 (Xze +... Xg(+-11 conditioned on X1, for any value /22
of l. We have V1 to be the conditional variance of (X2 + X3) on X1, which is 2/2 1. We consider recurrence. By the law of total variance
Var( 1 V2"
(X2n 4... X2ht1-))XX1)
=Var( 1
1
(X2h +... X2h X2 X3) Var(
[E[ + ... 1
2
/
The first part we can write as
Var (X2a + + X2 + X2h42h-1 + ... -1),,, V2" 2
Notice that X..,..., X2 12h-1-1 are the leaves of the subtree rooted X2, and 2hh-1,, Xgi+a-1 are the leaves of the subtree rooted at X3. Both subtrees with root X2 and with root X3 have the same topology tree of height h 1. So this term can be written as
Var( h- (X2a + ... X21420-1_1 + X2h ..... X2at1-1)|
2
12
1
1
(Var
22"h h- (X2a +... + X2 xx1)x)) + Var(
h-1 (X24424 1 1t ....., 1 (x))
V2"
2 (Vh-1+ Vn-1) -Vh-1
For the second term, have for every 2 <i< 2" 4 [(x X (the other noise this branch is mean
zero, and the coefficient of X2 is 1--1.) For identical reasons, for every 2 + -n < i 2 Xx] x X3. Therefore, we have
1
Var(
EE
X2n +.. X2tt1_1XX2 X3]
V
1
Var(
(E[X2) + .... X X2) + E[X 424-1 + Xe (Xa)))
V2"
Var(
1
x2 X2 + 2 -
Xz) (213)--1 - Var(V2 + V1) (2
2
So, by induction, we have V/, = (2X2)2 +... (222)--1 Finally, for all 2 si 2h+1 we have E(X,|X,) X1.
23

---

Learning Gaussian DAG Models without Condition Number Bounds
This yields
Var( V2"
(X2r +... + X2h+1-1))
1
-V/ + Var(
[B[(X2n +... -2ht1-1)](X1)
V2"
V/ + Var(
x (21)a X1) = 1/+ + (213)4 = -1+ + + ..... (2x²)).
V
This finishes the proof.
D.3. Proof of Lemma A.2
We will use the following standard result about concentration of the chi-square distribution. Lemma D.1. For Chi-square distribution, we have the following concentration inequality: if xi
+ & 8g//2
P(t > (1+ S e
and
P(t < (1 e)a: - ((1 )//2
This is folklore bound and can be proved by using the MGF of the chi-squared distribution. Now present the proof of A.2.
Proof of Lemma A.2. We know that 3 - 8r is multivariate zero-mean Gaussian with covariance matrix
12 x() )-1 This follows from the closed form solution of OLS. So know that the sum
Σ
((x))) (53 (83- B)) = - Σ j=1
x x56 x
(83 Bu)
1=1
follows X2 distribution. By Lemma D.I, take
4 1oR(1/6) 1, we have
P( (1+ S( k0122= (2 4log (1/6)/k
k/2
-. + 4 log (1/8)/K
2k/2
(1/81//)e
e/82/h
Let x = (1/8)2/k, we just need that 2 + 2 ln((x) S ea. This inequality holds for any x and can be proved by taking derivatives. So accordingly, m can be loggg1))) 2k)/e.
D.4. Proof for Lemma A.3
will use Lemma D.I. We know that E y2 follows the distribution 0² 2m so we just choose m such that the
probability of x 2 (1+ + E)m and x 2 (1+ are both most 8/2. Since (1 (1-s)e- Its we can focus on the upper
bound. We need to be (Lte) +ttym/22 < 8/2, which gives m 2 2log(2/8), Since 1/2, we have In(1 s) > _ So [n[
we have m 2 (log(2/6)/5)
D.5. Proof for Lemma C.3
We cite lemma of (van de Geer & Lederer, 2013):
Lemma D.2. (Theorem 1) Let Xi...., Xn be independent random variables with values R. and with mean zero. Suppose that for some constants a and K, one has
m! K
E X m <
m = 2.3.... 02
2
1=1
24

---

Learning Gaussian DAG Models without Condition Number Bounds
Then for all 1 0,
Kt Vni
P
X 2
20V2t+
2 exx((-t).
First we calculate the moments of Z. We can easily calculate B|Z,| = 3 21( 12
Let F(k) 2 1
n-
We have F(2) 5. F(3) उे
1/2, and also have that F(k 2) (A)
& < 1, have F(k) < 1 for all
k+2
integers & 2 2. This establish 1, K 1 in the Lemma. Therefore, by the lemma we have
2f
t
IP
Z
2exp(-1).
V11
m
m
We take t log(2/6) and m such that V +
Since
1/2, it suffice that < e/4, which implies
2 32 log(2/8)/e² suffices.
D.6. Comparison with (Gao et al., 2022)
As mentioned in Section 2.4, our lower bound is at least as good as the one provided in (Gao et al., 2022). To argue about that, we need to prove that M 1 2 dbzinn where S VM.
Without loss of generality, we assume the variance to be 1, or otherwise, we just scale the covariance matrix Recall that, the adjacency matrix B = {biz] is lower triangular after sorted by topological order, and the inverse covariance matrix = (I - (I - B) T Therefore, both 1g and (I T) are 1 since they are lower/upper triangular), and thus the smallest eigenvalue is at most 1.
inverse covariance matrix having as its determinant (because the determinant of
On the other hand, suppose X, has parents El X3, 1. 1. We can lower bound that
Xja.
We can write X; = Ei Xje bja + E. We take X? +
xx E bb
2X, Σ bo
k=1
By Cauchy Schwartz, we have the upper bound
(X; Xinbnnni)
2 5 (1+
Σ
62 (X? +
x2 Σx3) = 1+ Σ
(x,)=1+
1=1
1=1
k=1
k=1
This inequality can be achieved if
bij -1
X
Xja
v
1 Σd
k= 54-
1
kk=1
12
1
So, we proved that the maximum eigenvalue of θ is at least 1 + Ei 1ن 2 1 d12 eigenvalue of θ is at most 1, the ratio between the maximum and minimum eigenvalue of O is at least + all min- or, the condition number is at least db? db So, have concluded M 2 /1+ db² m
min Therefore, since the least
E. Discussion and Extension.
F. Identifiability for non-equal variance
By calculating the conditional variance, one can identify the topological order of the DAG and thus identify the topology. However, we here give an example where the equal variance assumption violated and one cannot find the topology the DAG, even up to the same Markov Equivalence Class.
25

---

Learning Gaussian DAG Models without Condition Number Bounds
Let b be small positive number. Consider two different DAGs G1, G2. variables X, Y. Z and X', Y', respectively. X,Y, The topology of G1 is x Y,Y - Z and the topology of G2 is X' Z' Z' Y", X' Y'. These two DAGs belong in different Markov Equivalence Classes. The variables in G satisfy the relations with equal variance:
== E1, Y = bX + 62, 2 = bY + E3
where 61,62, Ea are drawn i.idd. from ( distribution. Therefore, we can calculate the covariance matrix of the joint distribution (X, as
1 b
b
b
Cov(X, (X
b² +1
69 - b 63 +
b²
b² - b -
On the other hand, the variables in G2 satisfy
X -
X' +
Z' +
E3
12 +
82
where e S2, Ea are i.i.d. centered Gaussians, but the variances are 1, 1 + 1, F4I respectively. We can calso check that the covariance matrix of the joint distribution (X',Ya Y', Z') is
1
b
6
Cov Z") b
ba +
3 b
b2
b 1 b
b
+ b
follows that (x, Y, 2) and (X', Y', Z') follow the same distribution, hence the two DAGs are indistinguishable Notice that the minimum strength of the edges among the two networks is b², whereas the ratio of the variances in the second network is 1 + O(P). Thus, we have concluded our example to show pair of graph that is not distinguishable if the variance ratio is at least + (bmin).
E.2. Comparison to the PC Algorithm
A standard approach for structure recovery is the PC algorithm (Spirtes et al., 2001), which is based on conditional independence testing. In order to provide rigorous guarantees for correct recovery. most works rely on some form of the strong faithfulness assumption. It has been noted in prior work (Uhler et al., 2013) that this assumption might be too restrictive. Indeed, we next demonstrate that this assumption need not hold in our setting.
More formally, the strong faithfulness assumption requires that for any nodes a, and subset S S V, if, by the graph's structure, a, are not independent conditioning on S, then the conditional correlation of a, b conditioning on S is lower bounded. However, in our model this does not necessarily hold. We give the following example: consider the graph with three nodes X, Y, Z, and we have the model as
X = EI, 6Xi E2, = bY -62 +
Here, 83, are i.i.d. (O 1). We can derive that Z and X are independent, as we could write Z = bE S3. This happens due to path cancellation.
We complement the theoretical observations with experiments that show how our method compares to the PC algorithm. We use the PC algorithm package causal-lea (Zheng et al., 2024). We use the same setting as in the simulation experiments of Section 4, with the hyperparameters n 5, 10. d = and with bmax 1. We set 100 trials and the figure shows how much portion of graphs are correctly identified. All the graphs are random graphs We note that in general, the PC algorithm is only guaranteed to find the correct Markov Equivalence Class, since some orientations depend the order of increasing variance, which is extra information about the model that PC doens't "iees'ttnn use. Therefore, we consider that the PC algorithm correctly identifies the graph if it outputs any graph in the same Markov Equivalence Class. This means that it detects all colliders (i.e. a child with at least two parents) and the skeleton (edges without directions). In Figure 6, we present the percentage of times that the PC correctly identified the Markov Equivalence Class, as well as the percentage of times that it identifies the correct skeleton, which a relaxed goal. The result is as follows.
26

---

Learning Gaussian DAG Models without Condition Number Bounds
dm2, n
dm2, nm
0.
0.8
Algonthi
Algorithm
Algoritim2
Algoritim
Algoritim, edges correct
PC Algorthm, odges correct
Algorithm
Algorthm
0.3
0.4
0.-
0.3
200
300
600
500
600
700
800
900
1000
200
300
400
500
600
700
800
900
1000
d=3, n
d=3, n-10
1.0
1.0
0.8
Gso et al 22
Cap 22
Algonthm
Algenttim
Algonthm2
Algenthm
Algorittm, edges
Algorithm, edges correct
0.4
PC Algorith
PC Algorthm
0.
0.2
200 400
900
300 300
800 1000
Figure Comparison with the PC algorithm
In Figure 6, we notice that the PC algorithm behaves poorly even for the task of finding the skeleton of the graph, if the number of nodes becomes larger. This is explained by the observaton that the faithfulness assumption is less likely to hold if the number of nodes is larger. Note that we generate the edge weights with random signs, so as the number of nodes increases, the possibility of cancellations, such as the ones we saw in our previous example, also increases.
E.3. Adaptivity of Algorithm 1: the case of unknown 2² ,min or d
One reasonable question concerns the situation when one more parameters of the problem is unknown when we run Algorithm 1. Below, we explain some simple modifications that ensure the algorithm adapts to the unknown parameters. Essentially, we would like to argue that even we do not the parameters priori, if the algorithm is given enough samples it will succeed in identifying the graph.
Unknown 2 We notice that with the equal-variance assumption, min; Var(X,) o² Therefore, we can find the variable with empirical smallest variance, and take the empirical variance as 8², As we're given O(bmt))) samples, we can have with high probability, d² /02 1 to be + /10 (proven by Lemma A.3.) Thus, the error is small enough to make the Proposition A. hold, by slightly shrinking the threshold of judging whether node is the immediate parent. In particular, the following modification of Proposition A. holds, with a similar proof.
Proposition E.1 (with unknown )²). (1) T contains all the parents, then there exists some C and |J] S such that
27

---

end for end for Learning Gaussian DAG Models without Condition Number Bounds
the conditional variance of X is o². With high probability, it is no more than G2 /imii/ /10)
b2
(2) If'T does not contain all the parents, then there exists parent that is not in T, so for any coefficient vector ßT each X; - ST is Gaussian random variable that has variance at least (1 + biinn Specifically, for any J S T, for any coefficient BJ, X;- XJ also has variance 2 1 + bain- With high probability it is no less than o² + 96 /10).
Therefore, we can assume that x² is known and furthermore assume o² 1.
Algorithm 4 Estimating the DAG Information Theoretically input Samples (X(1), ,x))), ((()))....(())
output Topology
1: Phase 1: Finding the Topological Ordering, not knowing d and bmin-
2: T=
3: Find X1 argminy Var(X) and put i T
4: d=1
5: while [T' do
6: .M + 0, FOUND <- False
7:
for prriE e [n]][ do
8:
for / T, | < d do
9:
3J - Regress X; on X,
MSE m2 E, ( -X
10:
11:
if MSE < 1+ Olv
dlog then
12:
+ -appendd(
13:
FOUND + True end
14:
15:
if FOUND = False then
18: de 1-d+1
19: end if
20: end while
21: Phase 2: Finding the parents, knowing but not Dmnin
22: bmin +-
for E T do
S f- nodes in T before
25: for J C S, J| min( (()), do
26:
PASSED +- True
27:
for K E SKJ.KK < min([S] d) do
28:
BJUK + Regress X; on XJUK
29:
if E K : 1Bx 2 bain or 13x1 S b /40 then
30:
PASSED &- False, break.
end if end for
32:
33:
if PASSED = True then
34:
$ f- regress X; X,
35:
pa(i) + (j E J 13, bmin break.
36:
end
end for
if PASSED False then
39:
bmin bmin /2.go line 23.
40:
end if
41:
end for
28

---

Learning Gaussian DAG Models without Condition Number Bounds
For unknown Umin and d. Without knowing bmin- we need consider the following issue: because of the noise, the strength of the edge could be arbitrarily small. Therefore, we could never distinguish the case with faint signal versus no signal. Therefore, when sweeping through the possible values of 6min. we may enter region where all the detected edges have strength more than bmin. and the rest have much less. We may not infinitely search for mmiu- Our solution is to change the algorithm as Algorithm 4.
We first explain what is changed from Algorithm In Phase I, we add a variable FOUND to determine whether the current estimate d is enough to explain the variance of the next node in the ordering. We changed the threshold for MSE in line 11 so that it does not include d, min-- If we are provided sufficient samples, by Proposition A.I, we take the next node with small MSE. If no node has small enough MSE, then FOUND is False, and there is next node. We can conclude that the number of neighbors is not sufficient to make MSE small enough (there are missing parents), we need to increment d to d ++] 1 to fit.
In Phase 2, we changed the threshold in line 29. In line 29, we need to find all the edges to have strength either ômin or S 0minn400 for some K. we are given a sufficient number of samples, then with high probability, the ground truth for all those weights is either /min/ 22 or S /20 If for some min we cannot find such a K (so PASSED is False), it means that even if J U K covers all the parents, the coefficients are not within that range. Therefore, we need to make 8min smaller and try the algorithm again.
We can perform Phase of Algorithm 4 without knowing d and bmin- However, we need to perform Phase 2 of Algorithm 4, but not mmnin And we may change our guarantee to the following, which is comparable to Theorem 2.1. Notice that the number of samples doesn't change in Informal Theorem E.2 and only slightly increases in Informal Theorem E.3.
Informal Theorem E.2. Suppose we run Algorithm 4 using m independent samples generated from a DAG G according to (17). If Assumption 1.1 holds, then there exists an absolute constant C > 0, such that the following guarantees hold. Phase 1 of Algorithm succeeds in finding a correct topological ordering of the nodes with probability at least 1 5, provided that
1 (dlog()/d) + log(1/6))
m 2 C
1 min
Informal Theorem E.3. Suppose we run Algorithm 4 using m independent samples generated from DAG G according to (17), For a given node, we define a good set of parents as follows: there is a > O, such that for all reported i - j. |btg] 2 b and for other not reported biy. |b,yl S b/10. If Assumption 1.1 holds, then there exists an absolute constant > 0, such that the following guarantees hold. Provided that Phase of Algorithm succeeds in finding correct ordering, Phase finds the good set of parents for each node with probability least 1 6, provided that
2 C. (G)
(dlog(n/d) log(
8 1log(1/bmin) (bmin))))) )
1
122
min
We give sketch for proving the algorithm's correctness.
1. Sketch proof of Theorem E.2. From Lemma A.3, we know that the ratio between MSE and the variance is 1 X2 Therefore, the largest deviation should be
which can be a rubric for the conditional variance gap (i.e., if it is not the next node, then the conditional variance is O(hinn) which is O( O(Vdlog n/m). by Proposition A.1). Therefore, with high probability, this algorithm will terminate when d = d and finally find the true parents. The running time is at most d times the Phase of the original Algorithm 1.
2. Sketch proof of Theorem E.3, Since we are given a sufficient number of samples, by the proof of Lemma A.5 and Corollary A.6, we can derive that if J U K contains all the parents, then with high probability, the error of By - 3, is always & Obminn). ). Then, if we have all weights either bij 2 or bij S b/40, then with high probability |by 2 - / 80 or |bal 2 b/40 + bminn/8 which is [byl] 2 b/2 or |biy] < b/20, so this achieves our guarantee. Notice that this algorithm terminates when brmin < bmin/ 2, so the number of trials of 6min will be OOlo((1/ (()min)) and by union bound, we can have the algorithm succeed.
E.4. Extension to Sub-Gaussian Noise.
One further extension of the model is that we relax the structure of the noise Ey: the noise is Gaussian. We assume that for all the samples (x(1) x(() :..., x!k)), the model satisfies:
29

---

Learning Gaussian DAG Models without Condition Number Bounds
x
Σ )E
byixt((()) + E
(17)
(R) independent. Furthermore, for all the noise = ,.(m)) are independent sub- cau-Gaassinn matrices, Where all the noise Eq
where Varle, [ = 1 and the uub-Gaussian norm no more than 1r. Now, we show that we can perform Phase 1.
We start with the proof of Theorem 2.1 Section
finish the proof of finding the ordering, we take the same proof by measuring V E and BE. It boils down to proving that with this number of samples, V E and BE are upper bounded by VE < (1 (bbm))) Var[[x [x((r)) [ar]] and also BE S (1 Var[X][a.5]] The correctness can be guaranteed by Lemmas A.2 and A.3, so it just suffices to prove these two lemmas when Gaussian noise is substituted by sub-Gaussian noise, and the inequality holds for some different Cv and CB.
Before that, we use the following Haason--right Inequality (Rudelson & Vershynin, 2013).
Lemma E.4 (Hanson- (Hanson-Wright), Let X = (X1,..., "51,...,,,,,, Xn) E RR" be a random vector with independent components X₃ which xatisfy E[X:] = and (xll Let be an a x n matrix. Then, there a universal constant such that for every t t2 0.
2
P AX - E[X AX] > 1 S 2xxp -emin cmmn
KI
K2
Proof of Lemma A.2. We know that / 3 is multivariate zero-mean Gaussian with covariance matrix Xf )-1 This follows from the closed-form solution of OLS. So we know that
m [((x))) - - x (003 By). -B).
j=1
Let X be the matrix obtained by vertically stacking x5" x(² ..... x(m) and let Y be the column vector y
(x.") xx(3) ((m)))
m)) The error vector is also column vector, E = (5, ..... (mm) By OLS, we have
(x xx=
(XB* +0),
which implies
m
Σ
-,))² = =x((XX -
Therefore, let = XXX-XXXXXXXX which is projection matrix. We know that ||A|1 k. ||AAllop 1, and K th. Also, since 8 are independent noise with variance 1, we have Ele As] = k. Accordingly, we can set
((log((1/0)))) k) ve/e),
and that suffices.
Proof of Lemma A.3. The sum of squares of all random variables is a special case of the Hanson- Hanson-Wrigh Wright Inequality (Rudelson & Vershynin, 2013) with A = Im- Therefore, for the expected sum of squares of all y-(k'ss we have
P
(A)2 mV > +V2 < 2eexp - min
V m w3 $1 )
Hence, we know that
(1og 23)
suffices.
30

---

Learning Gaussian DAG Models without Condition Number Bounds
F. Further Simulation Results
The data setting is the one described Section 4. We provide here some more simulation results. We have two confi gurations of data. We vary different d, n, maxx N parameters and test the results.
d-2, n-1
d-2, n=15
d-2, n-20
d-2, nm25
d-2, n=30
1.00
1.0
1.0
1.0
1.0
0.9
0.8
0.8
0.8
0.90
0.8
0.7
0.6
0.6
0.80
600 et 22
530 ct 22
0.
600 et '22
630 al
G20 et
Algonttmm
Algorihmm
Aggnittm
Alporithm
Algorithm
Algorthi 2
0.6
Alçorthm 2
Aigorthm 2
Algorthm
Algorthim
0.5
200
400
600
500
1000
200
400
600
500
1000
200
400
500
800
1000
200
400
600
800
1000
200
400
500
800
1000
1.00
1.0
1.0
1.0
1.0
0.93
0.8
0.0
0.8
0.8
0.6
0.00
0.4
0.75
Geo et al '22
Aggonthm 1
Gao et '22
Alçonthm
0.2
Geo at '22
Aignnthm
0.7
Gao et al
Alporthm
0.2
Gao et al
Alporitmm
Algorthn
Algorthi
Aigonthm
Algorthi
Algorthi
0.0
0.0
200
400
600
800
1000
200
400
600
500
1000
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
1.0
1.0
1.0
3.0
1.0
0.9
0.8
0.6
0.8
0.0
0.6
0.4
0.4
0.4
Gao '22
Algorthm
Algorithm
0.2
Gao et '22
Alporthm
Alcorithm
0.2
Gao et '22
Aigorthm
Algorittim
Gao et al '22
Alporithm
Oso al
Algorthm:
0.0
600
800
1000
200
600
100
1000
200
400
500
880
1000
200
400
830
800
1000
200
400
500
800
1000
Figure 7. Figure for dmax
Figure shows a of simulation results with dmax = The test metric is test 100 different random graphs, and for each graph test for those graphs whether (Gao al., 2022), Algorithm and Algorithm discover the same topology.
We also performed another set of experiments for bmax 2 and the results are shown in Figure 8.
The test metric is to test 100 different random graphs, and for each graph test for those graphs whether (Gao et al., 2022), Algorithm 1 and Algorithm 2 discover the same topology. Also, we add the number of false positive edges with 95% confidence interval (centered with mean and 4:20 error region) as in Figure 9. Because of the long runtime, for the case when d = 4 and n = 30, the experiment Algorithm is not performed.
For all the cases we can find that eventually, when the number of samples m grows, Algorithm and Algorithm converge at faster rate than (Gao et al., 2022). We
that for d 1 case all the algorithms perform better then for bmax = 2, which is expected from the theoretical results, since T(G) and R grows much faster (as an exponential function) when b is larger than 1. For both (Gao et al., 2022) and Algorithm 2, we see degradation in performance as n increases, while comparatively, Algorithm is not affected as much. This can be explained by the growth of т, which seems to be sublinear in n, in contrast with R which grows linearly and the condition number which grows superlinearly (see also Figure 1).
31

---

Learning Gaussian DAG Models without Condition Number Bounds
d=2, n010
d=2, n-15
d=2, n-
d=2, n-25
1.3
10-
10
10-
0.9
0.8
0.8
0.3
0.6
0.6
0.6
0.7
0.4
0.4
0.5
Gao
Algoritim
'22
Gao al
Algorithm 1
Geo al '22
Algorithm
0.-
Gao al '2
Algoritim
0.5
Algorthm
0.2
Algorthm
Algorthm 2
0.0
Algorthm 2
600
800
1000
200
600
800
1000
200
800
1000
200
400
800
1000
d=3, n=10
d=3, n=
d=3, n=
d=3, n=
1.3
1.0
0.8
Gao et '22
0.3
Algoritmm
Algorthm 2
0.8
0.0
0.6
0.7
0.4
0.5
0.4
0.5
0.2
Gao al Algonitmm Algorthi
Algorihim
Algorithm 1
Algorthm
Algorthm
0.0
600
800
1000
200
800
1000
200
800
1000
200
400
1000
d=4, n=10
d=4, n=15
d=4, n=20
d=4, n=25
1.3
0.8
Gao et al '22
0.6
Gao et '22
0.6
Gac et '22
Algorithm
0.5
Algorithm
Algorithm
0.3
Algorthm ?
Algorthm ?
Algorthm 2
0.6
0.4
0.5
0.3
0.4
0.2
0.-
Algorithm 1
0.1
0.2
Algorthm
0.0
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Figure 8. Figure for Umax 2
dm2, nm10
d=2, n515
d=2, no
d=2, n-25
dm2, n=30
1.0
Gao et al 122
2.5
cao et 22
Gao al '22
Cas et '22
Gio et al 22
123
Algonthm 2
Algenithm
2.4
algorhhm
Algonthm
Algorthm 1 Algorthm
Algonitim 1
Algorthm
Algonthm
Algonthim
1.5
0.75
3.0
0.0
025
0.3
1-
200
600
500
1000
200
400
600
800
1000
200
400
600
800
2000
200
400
600
1000
200
500
800
2000
d=3,
6-3, n0
d=3,
al Algorithm
Gao 22
Gso '2
Cao et '22
15.0
al
Alpor thm
Algorittm
Algorithm:
Alporitimn
2.9
Algorthm
Algorthm 2
Algorthmm 2
Algorthm 2
Alporthm
10.0
6-
5.
1-
2-
2-
0.0
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
2000
200
400
600
500
1000
200
400
600
800
1000
d-4,
n=10
d-4, n-
15
0-4, 2-
d-4, n-25
d-4, n-30
Gan al 22
Gao '22
'22
30-
al
Algorithm
Algorithm
Algoritim
Algonttim
Algorthi
2.5
Algorthen
Algorthmm
0-
15
3-
0-
5-
0.9
0-
10000
2000
Figure 9. Figure for bmax = 1
32

---

Learning Gaussian DAG Models without Condition Number Bounds
G. Code
import numpy as np
from Selleann.iinee model import Lasso
from sklearn. linear_modell import LinearRegression
import random
def [mmple_d(stt d)
Sample <=d things in the list n len (1st)
if n==0: return []
10
parents = 1I
propability (d/ (n+2) 1/22) for node 1st:
12
13
if random. random () probability:
14
parents. append (node)
return parents[:d]
17 def
d, b o_min, b max)
18
Sample a random Bayesnet condiguration of n vertices with <=d indegree order = 1ist (range (n))
19
random. shufffl
coefficient (i:() for i range (n)
22
parents ():set() for in range (n)
23
B [[0] *n for in range [n)]
24
for i in range (n):
node orderlil
26
prt sample d(orderl:i],dt _ample_d(orderr]
27
parents nodel set (prt.copy
28
for p in prt:
abs_value (b m_max-b_minn rrandon. ()
30
(0,1)-1)) sign ra*damo bij abs_value sign
31
32
[[p] coefficient [node] [p] bij
return order, coefficient, parents, np.array (B)
35
def generate_samples (B, N)
36
37
From the sample configuration sarple set of values. N is the number of samples, B is the coefficient matrix. //'
40
n B.shape
raw = np. random. normal (0, 1, [N, 41
42
return rawênp. linalg. inv (np.eye (n)-B)
44
def combination (1st, d):
45
List all subsets (as list) with <=d clements in 1st if == 0 or len (1st) == 0:
46
return [[]]
else:
19
return combination (1st((1;),,,
[[ist [0]]+] for in combination (1st [1:), d-1)
50
def choose (1st, d):
List all subsets (as list) with =d elements in 1st if 0: return [[]] elif d > len (1st):
54
55
return []
else:
58
return choose (1st [1:], d) [[]st[0]]+] for in choose (1st [1:1, d-1)
59
60
def 1r (data, cand, index)
///
Do OLS linear regression, given the data, candidate nodes (cand)
63
regress on and the index of the node whom to regreas on. Return the
33

---

Learning Gaussian DAG Models without Condition Number Bounds
coefficient (if cand nonempty) and the MSE
65
66
assert index not in cand
N = data. data. shape [0]
datal;, [index]]
69
data(;,
70
cocf np. mp.linalg. inv (X.T3X) 8 (X.IeY)
71
mse ( ((Xeceef-Y) **2) sum () /N
return ccef.I,
Listing 1. Data-sampling algorithm
def
topo (data, d)
Kind the topelogical order given the data and d. Return the topo og cal order list of nodes. ...
N, data. shape
rem nodes = list (range (n))
top [np. awg.in (sum (d)
ren_nodes.remmoeecooo (0])
10
for range (n-2):
11
best_mse = 1_000_000
12
next_ncde -1
for in rem rem_nodes:
1_000_000
15
for cand combination d)
16
coef, msel = ir (data, cand, F] nse = min mse))
17
if mse best_mse:
19
best_mse mse
20
next_node r
21
rem sem_nodes: remove top append (next_ _os_____dee_ top.apppe ano.appendddcc reee_modes____eeeeeeeeeeeeessees [0]) return top
remove (next_node)
24
Listing 2. Phase 1 of Algorithm topological order
(Gao et al., 2022)
def learning parents_g (data, order, d)
Using Gac' s algorithm to find the patents given the topological order. Return the partents for each node as dictionary. //,
n = shta. [1]
parents = (i:set for i in range (n))
for i in range (n)
if i==0:
10
continue
11
prev = order [:i]
12
node = order [il
mse 1_000_000
carent_superset
15
Find c_j to be argmin (v_jC) for cand in combination (prev, d):
16
ccef, msel ir (data, cand, node)
if msel mse:
19
mse msel parent_superset = cand
20
parent sct()
for pa in parent_superset: parentml
copy ()
24
parentm1 remove (pa)
34

---

Learning Gaussian DAG Models without Condition Number Bounds
25
ccef, msel (data, parentml, node) if abs (mse-msel) gamma:
26
27
parent.ade (pa) parents [node] parent
parents
Listing 3. Finding the parents using (Gao al., 2022) algorithm
def leaming_parens__vddte order, d):
Using inefficient algorith thm find the patenta given the topo ogi order. Return the partents for each node as dictionary.
= shapeel
parents (1:set for i in range (n)
for range
cache - ()
10
if 1==0:
continue
12
prev order[:i]
13
node order [i]
14
parent_superset = None if i <= d:
15
16
parent_superset - prev
else:
18
for J thoose(prey) d):
19
passed True
others = list (set (prev) -set (J)
for K in choose (others, min (d,i-d))
22
JKsorted
(sortedd
23
if JKsorted not in cache. keys ()
24
I print ("new", JKsorted)
coefJK, msel Ir (data, sorted (J+K), node) cache [JKsorted] coefJK
27
# print (coef)
28
else:
# print ("cache")
coeЛJK cache [[oorted]
ck-[abs (Cc]] [m]) for range (1en (J+K)) if [Ksorted[m] in K]
12
if max (cK) >- b_min/2:
33
passed False break
if passed: break
36
37
parent_supersee =
38
cocf, msel = 1r (data, parent_suporset node) parent set (n for (c,n) in zip (coef [0], parent_suporset) if abs (c) >= b_min/2)
40
parents [node] parent return parents
41
Listing 4. Finding the parents using the Phase 2 of Algorithm 1
def lasso (data, cand, index):
data. [0]
Y data [;, [index]]
X = datal:, cand]
YI datal:, [index]]
X1 datal: cand]
lambda_n 0.01
if cand = []:
return [], (Y++2) sum ()/Y.size
reg Lassc (fit_intercept False, alpha lambda_n, copy_X Faise) fitt (X, Y)
coef reg.
12
return ccef, ((Xenp ([ccef]) .--)) x*2) sum () /N
35

---

Learning Gaussian DAG Models without Condition Number Bounds
13
14
def
(data,, d)
15
Learn the parents by LASSO, the efficient algorithm
N, = shape
rem_nodes list (range
18
top anp.argmin (sum (data**2)
19
rew.nodes.rcmove..oo (top (0])
20
parents = (1:set for i in range (n)
for in range(m-1)
22
best_nede -1
23
parent set
24
best_mse - 1_000_000
for node in rem_nodes:
coef, if mse best mse:
lasso (data, top, node)
27
28
best_mst mse = mse
29
best_node node
parent
(n for in zip (coef, top) if abs b_min/2) parent set (list (parent) [:d]]
32
(best_nodee_ parent top. app.append (best_node) rem_ncdes. remove (best_node) return top, parents
33
Listing 5. Algorithm 2: the topological order and parent
5 $ Choose any want Choose any you want
b_min 0.5
max 1 Choose any max you want
sigma = 1
gamma b_min++2/2
Gao' et al 2022 repeat 45 Choose the number of trials you want
def doing_regular (tup):
Genetare process to perform all three of those algorithm tup is of integers (i, N) where is a dunmy parameter and the number of the samples.
12
13
14
Return the result of whether they have find the correvt parents for three algorithm.
17
i, N tup
18
if j&20 = 0: print (j)
19
order, coefficient, parents, generate_contig (a, d, b_min, b. max)
data ap.array (generate sample (2, N) order (dattd)
21
22
g Tearningg_areent__ (data, order, d) parents
23
learning_parentt___vde (data, order, d) -- parents
learming_parents_e (data, d) [1] parents
return g, v, Gao al's, inefficient, efficient
Listing 6. Paremeter settings and function to do the experiments
36


Reviews:
[{"review_id": "Y5iQDYGUEm", "reviewer": "Reviewer_QEWp", "summary": "This paper studies the sample complexity of learning linear Gaussian DAGs under equal variance assumption. It proposes a new algorithm and provides the graph recovery guarantee independent of the condition number. Both the upper and lower bounds of sample complexity are proved. The authors also provide simulation to verify the proposed method and demonstrate improvement compared to prior results.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their encouraging words and detailed feedback. Below we answer the main points that were raised.\n\n# For comments:\n\n>For figure 2, it would also be good to include the error bar.\n\nThank you very much for this suggestion. We will make sure to include the error bars in the final version.\n\n>In Algorithm 1, it would be clearer to state how T is sorted and how to handle the topological order when two nodes have the same MSE.\n\nThank you very much for this question, it is very important to clarify how $T$ is updated since that might not be clear from the pseudocode. Essentially, $T$ is a list of nodes that is initialized as empty and every time we find the node with the smallest MSE we append it at the end of the list. If two nodes happen to have the same (smallest) MSE, then we choose either of them to append to $T$. This is because having the same MSE indicates they are both valid to be the next node in the topological ordering since they have the same conditional variance given the previous nodes, so we can pick any of them as the next node.\n\n>For the experiments, it would be useful to report the TPR and FPR of the edge discovery.\n\nAnswer: We agree this would be an excellent benchmark for performance, we will make sure to report these in the final version of the manuscript.\n\n# For questions:\n\n>To run Algorithm 1,2, knowing $b_{\\min}$ is required. How is this value estimated in practice?\n\nWe thank the reviewer for this insightful question. If we interpret the question correctly, this question is about how the algorithm would operate without knowing $b_{\\min}$ beforehand. We have included a detailed discussion in appendix E.3 about adaptivity when various parameters (including $b_{\\min}$) are unknown. We will make sure to make this discussion more visible in the main text. \n\nEssentially, without knowing $b_{\\min}$, the algorithm would proceed with an initial estimate $b$ of $b_{\\min}$. Then, we run the Algorithm assuming $b_{\\min} = b$ and in the end, we check whether the estimated strengths of all the edges are either less than $b/10$ or larger than $b$. If we detect some edge that has strength between $b/10$ and $b$, we half $b$ and run the process again. We can show that if our number of samples is the one specified by our main Theorems for Algorithms 1 ,2, then such a process will eventually terminate and return a $b$ such that all edges are either smaller than $b/10$ or larger than $b$. \n\n>In Lemma 2.2, the authors showed that $\\tau(G)\\le\\kappa$. Is there an example that $\\tau (G)$ has an order of magnitude smaller than $\\kappa$? Since Theorem 2.1 depends on  $\\tau (G)$. The theorem would be more convincing if there is an example demonstrating that  $\\tau (G)$ is of magnitude smaller than $\\kappa$.\n\nThank you very much for raising this point, we absolutely agree that such an example would help demonstrate the improvement over previous bounds. Such an example is essentially provided in Lemma 2.5 (ii), where the topology is a binary tree, where each edge has $2^{-1/4}$ weight (which makes $\\tau$ a constant), but the condition number grows as $\\kappa = \\Omega(\\sqrt{n})$. We will make sure to highlight this example immediately after presenting Theorem 2.1, so that it is clear to the reader why the dependence on $\\tau$ is preferred over $\\kappa$. \n\n>The simulation results are in the low-dimensional regime. It would be nice to see the result in higher dimensional setting (when the ratio of n/m is smaller, and n is larger).\n\nThank you for pointing out. We will make sure to include simulations with higher $n$ that better reflect the high dimensional nature of the problem.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes, it is evaluated on simulated dataset. The data generation process of the simulated data and the evaluation criteria makes sense to me.", "theoretical_claims": "I checked the proof of Lemma 2.2. It seems reasonable.", "experimental_designs_or_analyses": "Yes, I checked the simulation setting. The setup seems reasonable.", "relation_to_broader_scientific_literature": "Learning graphical models can help drive deeper understandings for the physical process of biology, neuroscience and so on. Applications in this domain often have limited sample size. Hence, understanding the sample complexity of the algorithm helps better assess the quality of the estimation and understand its limitation.  Furthermore, this paper proposes a new algorithm that is shown to recover the graph without dependence in condition number, which often grows with respect to the number of parameters. This makes estimation in high-dimensional setting with small sample size possible.", "essential_references_not_discussed": "Not that I am aware of.", "other_strengths_and_weaknesses": "1. This paper is clearly written and the problem it studies is interesting and contributes to the theoretical aspect of graphical models.\n2. Algorithm 1 suffers from computational complexity. Algorithm 2, while more efficient, has some loss in sample complexity.", "other_comments_or_suggestions": "1. For figure 2, it would also be good to include the error bar. \n2. In Algorithm 1, it would be clearer to state how T is sorted and how to handle the topological order when two nodes have same MSE.\n3. For the experiments, it would be useful to report the TPR and FPR of the edge discovery.", "questions_for_authors": "1. To run Algorithm 1,2, knowing $b_min$ is required. How is this value estimated in practice?\n2. In Lemma 2.2, the authors showed that $\\tau(G)\\leq \\kappa$. Is there an example that $\\tau(G)$ has an order of magnitude smaller than $\\kappa$? Since Theorem 2.1 depends on $\\tau(G)$. The theorem would be more convincing if there is an example demonstrating that $\\tau(G)$ is of magnitude smaller than $\\kappa$. \n3. The simulation results are in the low-dimensional regime. It would be nice to see the result in higher dimensional setting (when the ratio of n/m is smaller, and n is larger).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xuliDUHNJu", "reviewer": "Reviewer_vPQm", "summary": "This submission investigates the estimation of Gaussian DAG structure from i.i.d. samples and under equal-variance assumption. The main findings are two-fold: the authors give a polynomial time algorithm to recover the structure based on ideas of sparse regression; the authors shows that the sample complexity can be expressed in terms of the maximum sump of squares of the coefficients for the out-edges. Furthermore, they show that this complexity parameter improve upon standard analyses based on conditional number, as this latter may grow polynomially with the dimension. Hence, their method is particularly interesting in high-dimensional settings where the number of nodes of the graph is large.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their encouraging feedback and for pointing out these omissions, we will make sure to correct them in the final version.", "claims_and_evidence": "The submission developed the theoretical aspects under reasonable assumptions. The authors proved upper and lower bounds, and sample complexity to support their analysis. Numerical experiments on synthetic data confirms the benefits of their method in high-dimensional settings.", "methods_and_evaluation_criteria": "Evaluated on synthetic data", "theoretical_claims": "I have checked the sample complexity of Algorithm 2 (Theorem 2.6), various lemmas in the main document.", "experimental_designs_or_analyses": "I did not run their code but I believe that their reported experiments are sounds.", "relation_to_broader_scientific_literature": "This paper investigates DAG models, applications of DAGs can be interested in these results and algorithms.", "essential_references_not_discussed": "Not to my point of view.", "other_strengths_and_weaknesses": "The submission is well written and delivers a quite exhaustive analysis going from information theoretic result to polynomial time algorithm with guarantees.", "other_comments_or_suggestions": "Typos:\n024 space between : PC-algorithm(Spirtes et al., 2001)\n070 space between : model selection(Spiegelhalter\nother spaces between word+ref, please check.\n149 variance should be one: \\sigma^2=1, wlog \n368 \\in missing : jpa(i)", "questions_for_authors": ".", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DmlkxH5BKr", "reviewer": "Reviewer_TgLM", "summary": "In this paper, the authors revisit the problem of learning an $n$-variate Gaussian graphical model from samples. This problem is known to be solvable in $O(d\\log n)$ samples where $d$ is the degree. However, there is a hidden polynomial dependence on the condition number of the covariance matrix of observations which can be polynomial in $n$ in the worst case. \n\nInstead, the authors make an assumption that the sum of squares of the linear SEM coefficients ($\\tau$ parameter) of a node is bounded. Together with this, they also have the usual assumptions such as each noise is standard Gaussian and the coefficients are not too close to 0 ($b_{min}$ parameter). Under this assumption, they manage to give an algorithm that runs in $O(n^{2d+2})$ time and $O(d\\log n)$ samples. They also show using Fano's method that their dependence on the previous two parameters is also optimal. Although, there is a gap of $d$ between the upper and lower bounds. \n\nSubsequently, they give an alternate algorithm that runs in poly(n,d) time, but its sample complexity is worse than that of the first algorithm. They further show that their first assumption is weaker that the condition number assumption. Technically, their algorithm is split is two phases: 1) topological order recovery using a greedy method that picks the lowest variance first  2) parent set recovery using regression using an existing result for undirected graphs.\n\nThe authors perform experiments with synthetic data to confirm their findings. The results show that their accuracy for graph recovery is better than the existing methods once the sample size becomes large.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their detailed response and thorough examination of our work. We answer the main points that were raised below.\n\n>Although, it looks to me too few (only one) prior work has been compared against. Section E in the appendix had some more such comparisons. Maybe some of it could be brought in the main body.\n\nThe reviewer is absolutely correct, the manuscript would benefit from some further discussion of prior work in the main body. We will make sure to bring some of this discussion earlier instead of only appearing in the appendix.\n\n>It would have been more convincing if they also had experimental results for benchmark datasets along with synthetic datasets.\n\nThank you for raising this point. We agree that including benchmark datasets in the simulations would complement the theoretical picture even more. Our focus in this work was to theoretically establish the optimal sample complexity for this problem and to demonstrate how that translates to improvements over previous algorithms. Since the algorithms we were comparing against were tested on synthetic datasets, we made the same choice. \n\n>In section 2.4 the authors interpret their theoretical results in comparison with the existing literature. It would have been also nice if they can point out their main technical novelty that gave them improved experimental result. What is the new thing they did algorithmically that gave the improvement? If they discuss this in more detail, it will make the paper more understandable.\n\nThank you for raising this important point. We will make sure to clarify the experimental improvement in the final version.  Our algorithm for finding the parents of a given node $X_i$ is based on a connection that we observed between our problem and that of learning undirected gaussian graphical models. Thus, it proceeds by choosing a candidate neighborhood $S$ of size $d$ and then regressing $X_i$ with $X_{S\\cup T}$, where $T$ ranges over all other possible neighborhoods of size $d$. It only accepts a candidate neighborhood $S$ as the true one if for all these regressions we never find a node outside $S$ with large coefficient. In contrast, the previously known algorithm chooses the subset $S$ with the smallest Mean Squared Error (MSE) in the regression of $X_i$. This is a superset of the true neighborhood, so they find all nodes in $S$ such that removing them from $S$ increases the MSE by less than $\\gamma$ and delete them. Unfortunately, their choice of $\\gamma = \\Theta(b_{\\min}^2)$ is too big, resulting in many true parents being deleted. In our analysis, we quantify the correct threshold for deletion and show that it also depends on $\\tau$. This difference has a significant impact on performance.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes. \n\nAlthough, it looks to me too few (only one) prior work has been compared against. Section E in the appendix had some more such comparisons. Maybe some of it could be brought in the main body.", "theoretical_claims": "I briefly looked at section A and B. It looks okay at a quick glance. I didn't go through the proofs in detail.", "experimental_designs_or_analyses": "It would have been more convincing if they also had experimental results for benchmark datasets along with synthetic datasets.", "relation_to_broader_scientific_literature": "I think the claim that $\\tau$ is a better assumption has been theoretically justified. The proposed algorithms is also performing better empirically.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "In section 2.4 the authors interpret their theoretical results in comparison with the existing literature. It would have been also nice if they can point out their main technical novelty that gave them improved experimental result. What is the new thing they did algorithmically that gave the improvement? If they discuss this in more detail, it will make the paper more understandable.", "other_comments_or_suggestions": "NA", "questions_for_authors": "see above", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper studies learning of Gaussian linear structural equation models under equal variance. Algorithms are proposed. For n nodes and degree d, theoretical analysis shows that the sample complexity (lower and upper bound) is O(d log n) which has been observed before in Gao'22. The truly novel result is that the bound is now independent of the condition number.

The reviewers agree to accept this paper, although the time complexity for the proposed algorithm is O(n^(2d+2)). A more efficient LASSO estimator is provided, but for the new algorithm the sample complexity is O(d^4 log n).

Please take into the account the comments from all reviewers regarding for instance experiments with larger number of nodes n, or point out the main technical novelty that allowed for improved experimental results.
Decision: Accept (poster)