We thank the reviewer for their detailed response and thorough examination of our work. We answer the main points that were raised below.

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.

The 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.

It would have been more convincing if they also had experimental results for benchmark datasets along with synthetic datasets.

Thank 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.

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.

Thank 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.

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.

We thank the reviewer for their encouraging words and detailed feedback. Below we answer the main points that were raised.

For comments:

For figure 2, it would also be good to include the error bar.

Thank you very much for this suggestion. We will make sure to include the error bars in the final version.

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.

Thank 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.

For the experiments, it would be useful to report the TPR and FPR of the edge discovery.

Answer: We agree this would be an excellent benchmark for performance, we will make sure to report these in the final version of the manuscript.

For questions:

To run Algorithm 1,2, knowing $b_{\min}$ is required. How is this value estimated in practice?

We 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.

Essentially, 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$.

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$.

Thank 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$.

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).

Thank you for pointing out. We will make sure to include simulations with higher $n$ that better reflect the high dimensional nature of the problem.