Fast Scalable and Accurate Discovery of DAGs Using the Best Order Score Search and Grow Shrink Trees
Authors: Bryan Andrews
Joseph Ramsey
Ruben Sanchez Romero
Jazmin Camchong
Erich Kummerfeld
Venue: NeurIPS 2023 poster (2023)
Keywords: Causal Discovery
Directed Acyclic Graphs
DAGs
fMRI
Graphical Models
High Dimension
Densely Connected
TL;DR: We present new methods for causal discovery that are fast scalable and accurate. These methods are validated in many conditions including fMRI data.

Abstract: Learning graphical conditional independence structures is an important machine learning problem and a cornerstone of causal discovery. However, the accuracy and execution time of learning algorithms generally struggle to scale to problems with hundreds of highly connected variables---for instance, recovering brain networks from fMRI data. We introduce the best order score search (BOSS) and grow-shrink trees (GSTs) for learning directed acyclic graphs (DAGs) in this paradigm. BOSS greedily searches over permutations of variables, using GSTs to construct and score DAGs from permutations. GSTs efficiently cache scores to eliminate redundant calculations. BOSS achieves state-of-the-art performance in accuracy and execution time, comparing favorably to a variety of combinatorial and gradient-based learning algorithms under a broad range of conditions. To demonstrate its practicality, we apply BOSS to two sets of resting-state fMRI data: simulated data with pseudo-empirical noise distributions derived from randomized empirical fMRI cortical signals and clinical data from 3T fMRI scans processed into cortical parcels. BOSS is available for use within the TETRAD project which includes Python and R wrappers.

Fast Scalable and Accurate Discovery of DAGs Using the Best Order Score Search and Grow-Shrink Trees
Bryan Andrews
Joseph Ramsey
Department of Psychiatry & Behavioral Sciences University of Minnesota Minneapolis, MN 55454
Department of Philosophy
Carnegie Mellon University
Pittsburgh, PA 15213
andr1017@umn edu
jdramsey8andrew. cmu edu
Rubén Sánchez-Romero Romero
Center for Molecular and Behavioral Neuroscience Rutgers University
Newark, NJ 07102
ruben. www.orutgrrrrrrr.ooooooo edu
Jazmin Camchong
Erich Kummerfeld Institute for Health Informatics University of Minnesota Minneapolis, MN 55454 eri chk@umn edu
Department of Psychiatry & Behavioral Sciences
University of Minnesota Minneapolis, MN 55454 camch0020umn edu
Abstract
Learning graphical conditional independence structures is an important machine learning problem and cornerstone of causal discovery. However, the accuracy and execution time of learning algorithms generally struggle to scale to problems with hundreds of highly connected variables-fof instance, recovering brain networks from fMRI data. We introduce the best order score search (BOSS) and grow- shrink trees (GSTs) for learning directed acyclic graphs (DAGs) in this paradigm. BOSS greedily searches over permutations of variables, using GSTs to construct and score DAGs from permutations. GSTs efficiently cache scores eliminate redundant calculations. BOSS achieves state-of-the and execution time, comparing favorably a variety of combinatorial and gradient- based learning algorithms under broad range of conditions. To demonstrate its practicality, we apply BOSS to two sets of resting-state fMRI data: simulated data with pseudo-empirical noise distributions derived from randomized empirical fMRI cortical signals and clinical data from 3T fMRI seans processed into cortical parcels. BOSS is available for within the TETRAD project which includes Python and wrappers.
performance in accuracy
1 Introduction
We present a permutation-based algorithm, Best Order Score Search (BOSS), for learning a Markov cquivalence class of directed acyclic graphs (DAGs). novel method for caching intermediate calculations for cermutation-based algorithms, Srrw--hrrikk Trees (GSTs), is also developed and presented in this paper. Our implementation of BOSS using GSTs scales well in both accuracy and time to higher numbers of variables and graph densities. We demonstrate that in particular, this method scales least the complexity of dense cortical parcellations of fMRI data.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

In real world systems, it is not unusual to have hundreds or thousands of variables, each of which may be causally connected to dozens of other variables. Models of functional brain imaging (fMRI), functional genomics, electronic health records, and financial systems are just few examples. In order for structure learning methods have an impact on these real world problems, they need to be not only highly accurate but also highly scalable in both dimensions.
This work follows earlier research on permutation- structure learning 12
22
24
26 A
previous permutation-base method, Greedy Relaxations of the Sparsest Permutation (GRaSP) 12). demonstrated strong performance on highly connected graphs, but struggled scale to sufficiently large numbers of variables. By using BOSS and GSTs, we are able to overcome this challenge, while maintaining nearly identical performance GRaSP, As our simulations show, other popular methods fall short on accuracy, scalability, or both, for such large, highly connected models.
In the remainder of this paper, we provide background for our approach in Section followed by discussion of GSTs in Section We then introduce BOSS in Section and validate it by comparing several alternatives. We show that BOSS using GSTs compares favorably to a number best- performing combinatorial and gradient-based learning algorithms under a broad range of conditions, up 1000 variables with average degree of 20 in Section [S]
Finally, to demonstrate its practicality, we apply BOSS to two sets of resting-state fMRI data: simulated data with pseudo-empirical noise distributions derived from randomized empirical fMRI cortical signals, and clinical data from 3T fMRI scans processed into 379 cortical parcels in Section Section 7 provides a brief discussion.
1.1 Our Contributions
Our novel contributions to the field include the following:
1. We present new data structure, Grow-Shrink Trees (GSTs), for efficiently caching results of permutation-bant structure learning algorithms. Using GSTs dramatically speeds up many existing permutation-ban based methods such as GRaSP.
2. We present new structure learning algorithm, the Best Order Score Search (BOSS), which has similar performance GRaSP (and thus superior performance other DAG-learning methods), but is more convenient to use than GRaSP because it has fewer tuning parameters and faster and more scalable.
3. We prove that BOSS is asymptotically correct and provide implementation within the TETRAD project 17
4. We present validations of BOSS's finite sample performance via standard graphical model simulations and on both real and simulated fMRI data.
2 Background
The following conventions are used throughout this paper: V denotes a non-empty finite set of variables which double vertices in the graphical context, lowercase letters a, b, c.... V denote variables or singletons, uppercase letters A, B, C,... C V denote sets, and X denotes a collection of random variables indexed by V whose joint probability distribution denoted by P.
Probabilistic conditional independence between the members of X corresponding to disjoint sets A, B, C C V is denoted A IL B C P] and reads XA and X3 are independent given Xc.
2.1 Permutations
permutation is sequence of variables π = (a, c...... Let e V and i N < [V]). If is permutation of V, then is equipped with two methods: n. 7.index(v) which returns the index of in म, and ग moveeu, i) which returns the permutation obtained by removing v from я and reinserted at position i. Furthermore, pre, (v) {w V x.iidexxww) < (..index()) (1)) is the prefix of v.
2.2 DAG Models
directed acyclic graph (DAG) model is set probabilistic models whose conditional independence relations are described graphically by a DAG. In general, the vertices and edges of DAG represent

---

variables and conditional dependence relations, respectively. We use the notation 9 = (V, EE) where is a graph. is vertex set, and E is an edge set. In a DAG, the edge is comprised ordered pairs that represent directed edge and contains no directed cycles.
Let (V, F.) be DAG and V:
pag {w V - in 9} chg (v) {w E V
in 9}
the parents and children of v, respectively.
Graphical conditional independence between disjoint subsets A, B, C can be read off DAG using d-separation and is denoted A TT B [S]. This criterion admits the concept of a Markov equivalence class (MEC) which (in the context of this paper) is collection of DAGs that represent the same conditional independence relations.
2.3 Causal Discovery
DAGs are connected to causality by the causal Markov and causal faithfulness assumptions. A DAG is causal if it describes the true underlying causal process by placing directed edge between a pair of variables if and only if the former directly causes the latter.
Causal Markov: If 9 (V, E)) is causal for collection of random variables X indexed by V with probability distribution P. then for disjoint subsets A, B, C V:
A I1 B C
C|
More generally, this implication is called the Markov property and say DAGs satisfying this property contains the distribution. Moreover, a DAG 9 is subgraph minima no subgraph of 9 contains P.
Causal faithfulness: If = (V, is causal for a collection of random variables X indexed by V with probability distribution P, then for disjoint subsets A, B, C C.V.
A A
Al C[]].
Accordingly, under the causal Markov and causal faithfulness assumptions, finding the MEC of the causal DAG equivalent identifying the simplest model that contains the data generating distribution. From this point of view, causal discovery is standard model selection problem which
we address using the Bayesian information criterion (BIC) 7 21)
Let 9 be a DAG and X jid ~ P be a dataset drawn from member of a curved exponential family:
BIC X)= = Σ
BIC X., Xpag (w))
VEV
Σ VEV
lv (vg(()) (v) (fmiee
x)
1mi log(n)
where > 0 is penalty discount, é is the log- ligeliiel function, and 6 is the maximum likelihood estimate of the parameters; for details on curved exponential families, we refer reader to (9. Importantly, the BIC is consistent.
Proposition 1. Haughton 7 If P a member of curved exponential family and X ind P. then the BIC is maximized in the large sample limil by DAG models containing P that minimize the dimension of the parameters space.
Notably, Gaussian and multinomial DAG models form curved exponential families in which the causal DAG minimizes the dimension of the parameters space, and thereby maximize the BIC [IT] For these distributions, DAG models belonging to the same MEC correspond to the same DAG model,
'The edge v +- w corresponds to order pair (u, w))
This concept is also known SGS-min (GS-minimality 28
3

---

so it is common for algorithms to return a graphical representation of the MEC called a CPDAG rather than a DAG. Choosing between DAG models in a MEC requires additional information.
We greedily search over permutations of variables using grow-shrink (GS) 13 to project and score DAGs from permutations. Algorithm [ while Algorithm 3 project) details the process of projecting a permutation to a DAG.
(grow) and Algorithm (shrink) give the details of GS
Algorithm 1: grow X, Z)
Algorithm 3: project
Input: data X var prefix:
Input: data: X
perm
Output: parents: W
Output: D 9
W
E+0
repeat
foreach E π do
w 4 argmax BIC( Xwuz)
Z pre (v)
if w + Ø then
W + grow X, v, Z)
L + W
W shrink X, 8 W)
until w =
foreach w E W do
L E E (v, w)
Algorithm 2: shrink(X, W)
9 (V,E)
Input: data: : X var parents W Output: parents: W
repeat
argmax BIC BIC(X, Xw)))
if w # Ø then
L W - W w
until w Ø
However, how do we know that the causal DAG will even be constructed by project? This is guaranteed by the follow results. Let a member of a curved exponential family satisfying causal Markov and causal faithfulness:
If 9 is the causal DAG, then 9 is subgraph minimal
28
If 9 is a DAG and X ~ P then 9 is subgraph minimal if and only if there exists permutation ı such that pproorr(XX,.. in the large sample limit 12
These results admit the following strategy: (1) Search over permutations while using GS and the BIC to construct and score subgraph minimal DAG from these permutations. (2) Return the MEC of the subgraph minimal DAG that maximizes the BIC. The question is how to do so in a way that is both comprehensive and efficient.
3 Grow-Shrink Trees
Grow-shrink trees (GSTs) are tree data structures for caching the results of the grow and shrink subroutines of GS 13 This data structure is compatible with many pernutation-based structure learning algorithms, including BOSS and GRaSP. As an initialization step, a GST constructed for each variable in the dataset which are then queried during search rather than running GS. Each tree has a root node representing the empty parent set. Parents are accumulated by traversing edges of the tree, with every node in the tree corresponding to "grown" parent set. Each grown parent set corresponds to running the grow subroutine for some prefix. The shrink subroutine is also run and cached at each node of the tree.
The benefit of using GSTs is in how efficiently they store information needed for rumning GS. In the grow subroutine, nodes are added to the GST one at a time and scored. A new child node is added to tree for each possible addition, and the child nodes are sorted descending order relative to their scores. After all candidate parents have been added to the tree, the first child in the sorted list that is also in the prefix is chosen and the corresponding edge traversed. The shrink subroutine can then be run and the removed parents and scores can be cached.

---

(a)
(b)
(c)
(d)
(c)
Figure 1: Minimal subgraphs.
R<dRaaE
(a)
(b)
(c)
ceacbed
b<eRaed
(d)
(f)
Figure 2: A Grow-shrink Tree
3.1 An Example
This section walks through the usage of GST applied to toy example. Minimal subgraphs depicted in Figure of which (Ta) is the true DAG. Figure FigureZoutlines outlines the process of growing a GST for a while scoring permutations. In the following, the root. All ambiguities in sorting the branches of a node are resolved in alphabetical order. We proceed in sequence, but in theory some of these steps could be performed in parallel.
identify nodes in the GST by their path from
In what follows we run GS on a series of permutations while caching the results in a GST. We describe the nodes of the tree by the sequence of variables traversed to reach them from the root. We use colors to track our progress: light gray parts of the tree denote unexplored paths with no cached scores, dark gray parts of the tree denote explored paths with cached scores. Colored edges denote the paths considered during the execution of GS: red for rejected and blue for accepted. Colored nodes denote the same, but are only colored if they do not have cached value and require new score calculation. Numbered edges directed out of parent node record the order of preference for the children from high to low. To simplify this example, we do not cache the results of shrink. In (2a) we initialize the tree by adding the root node a.
Our first permutation to evaluate is (b, d, c) which is shown in (2b) and whose minimal subgraph is depicted in (Ib). Node a has not been expanded, so we score and sort nodes ab, ac, and ad. Node ab scored the highest and is contained in the prefix, we travel to node ab. Node ab has not been expanded, so we score and sort nodes abc and abd. Node abd scored the highest and is contained in the prefix, we have completed growing. We run shrink at node abd which removes nothing so the GS result is (b, d).
travel to node abd. At this point more variables are available in the prefix so
Our next permutation to evaluate is (c, d, a, b) which is shown (2ch and whose minimal subgraph is depicted in (lc). Node has already been expanded check its outgoing edges. Node ab
5

---

has the highest score but is not contained in the prefix. Node ad has the second highest score and is contained the prefix, so we travel node ad. Node ad has not been expanded, so we score and sort nodes ade. Node has the highest score and is contained in the prefix, so we travel to node ade. At this point no more variables are available in the prefix so we have completed growing. We run shrink at node ade which removes nothing so the GS result is {c,d].
Our next permutation to evaluate is (c, a, b, d) which is shown in (2d) whose minimal subgraph is depicted in (Td). Node a has already been expanded so we check its outgoing edges. Node ab has the highest score but is not contained the prefix. Node ad has the second highest score but is not contained in the prefix. Node ac has the third highest score and is contained in the prefix, so we travel node ac. At this point no more variables are available in the prefix so we have completed growing. We run shrink at node ac which removes nothing so the GS result is {c}.
Our next permutation to evaluate is (b, c, a, d) which is shown in (2eb whose minimal subgraph is depicted in (Te). Node a has already been expanded so we check its outgoing edges. Node ab has the highest score and is contained in the prefix, so we travel to node ab. ab has already been expanded so we check its outgoing edges. Node abd has the highest score but is not contained in the prefix. Node abe has the second highest score and is contained in the prefix, so we travel node abc. At this point no more variables are available in the prefix so we have completed growing. We run shrink at node abe which removes nothing so the GS result is (b.
Our last permutation to evaluate is (c, b, d, a) which is shown in 2f and whose minimal subgraph is depicted in (Ta). Node has already been expanded so check its outgoing edges. Node ab has the highest score and is contained in the prefix, so we travel to node ab. Node ab has already been expanded
check its outgoing edges. Node d has the highest score and is contained in the prefix, so we travel to node abd. Node abd has not been expanded, so we score node abdc. No scores result in an improvement to the overall score so we have completed growing. We run shrink node abd which removes nothing so the GS result is {b. d}.
Best Order Score Search
Similar to GES, BOSS uses two phase search procedure [5]. The first phase uses the best-move method, which takes a variable as input and greedily moves it to the position in the current permutation that maximizes the score. In this phase, best-move is repeatedly applied to each variable, one at a time, until there are no more moves that increase the score. This phase concludes with the f find-compelled procedure of [ which converts the DAG into CPDAG. The second phase of BOSS is BES which is exactly second phase of GES. The BES step is optional but guarantees asymptotic correctness if executed.
Algorithm 4: BOSS X, TT, o)
Algorithm 5: best-move TJ,
Input: data: X perm π. flag:
Input: GSTs : T perm π var
Output: graph:
Output: perm
J <- GST X)
best <- T.score(T)
repeat
for + | do
best + T. J.score(T))
T.index(U) т + ...moveuuu
foreach u E न do
L f best-move
best < J.score(z)) then
best + J.score(T)
until best J.score(z) g x J. V.project(m (ir)
else
π + n.move(u.j)
f- f find-compellede
if 5 = true then BE X)
Proposition 2. Let P be a member of curved exponential family satisfying causal Markov and causal faithfulness. If X iid P then BOSS X, true) returns the MEC of the causal DAG for all initial permutation in the large sample limit.
Proof. Since project returns a subgraph minimal DAG, it contains P. Accordingly, asymptotic correctness follows from the correctness of BES.

---

In Algorithm| (BDSS) and Algorithm 5 GST constructs collection of GSTs, denoted T, which contains one GST for each variable. The collection T is equipped with project and score methods which runs the GS algorithm to project and score a permutation, respectively.
5 Simulations
We evaluated the speed and performance of BOSS on simulated data compared to other algorithms: GRaSP, FGES, PC, DAGMA, and LiNGAM. Our evaluation used the performance metrics tabulated in Table which were originally proposed by [IO]. All algorithms were run on an Apple M1 Pro processor with 16G of RAM. The results reported in the main text are abridged but the complete results are available in the Supplement.
Table I: Metrics
True
Estimated
Adjacency
Orientation
tp
tp, tn fp, fn fn fn
Precision
tp tp + fp
tp
<- b
tp
a. a...b
fn
tp
+
fp
fp
Recall
fp fp
tp + fn
a...b
fp
We evaluated our implementation of BOSS using a BIC score with À 2 and BES step did not appear to improve performance. For GRaSP, we modified (to use GSTs) and used the TETRAD implementation with the same parameters the authors and BIC score with 2. 12 17 For fGES we used the TETRAD implementation with default parameters and a BIC score with = 2 15 177. We also used the implementation of in TETRAD with default parameters using BIC score with À = 2 as conditional independence oracle [17 (25. For DAGMA we used the authors' Python implementation with the parameters reported their paper but changed the threshold 0.1 2 and used the technique described in [4] to resolve cycles. Moreover, the output of DAGMA was converted CPDAG for linear Gaussian simulations. For LiNGAM we used the authors' Python
implementation with default parameters g 22
We generated Erdös-Rényi DAGs by applying arbitrary order vertices of an Erdös-Rényí graph. For scale-free networks, first we generated an Erdös-Rényi DAG and then redrew the parents of each vertex according the Barabási-Allert model [T. This resamples the out-degrees distribution be scale-free while keeping in-degree distribution constant. This procedure was motivated by the results Figure (5) Edge weights were sampled from a uniform distribution in the interval [-1.0, 1.0] and error distributions were generated with standard deviations sampled from a uniform distribution in the interval [1.0, 2.0]. The complete simulation details, including figures plotting the simulated scale-free free in/out degree distributions, are include in the Supplement. Additionally, the data are available for download at: https: //github. com/ cmp/ httpcccupppphi/ wttpaauwccoommmmm/////ooooo//
3aaure33comprres BOSS against GRaSP, 4GES, PC, and DAGMA on linear Gaussian data generated from (3a) Erdös-Rényí and (3b) scale-free networks. The output of DAGMA converted to CPDAG for these simulations. Lam et al. [12 attribute the excellent performance of GRaSP to it being robust against the ubiquity of almoss-violations of causal faithfulness [29 27). Due to the algorithmic and performance parallels between GRaSP and BOSS, we conjecture that a similar argument could be made for BOSS.
Figure 4a compares BOSS against DAGMA and LiNGAM on linear exponential and linear Gumbel data from scale-free networks. Interestingly, LiNGAM, which takes advantage of non non-Gaussian signal in data, does not perform appreciably better than BOSS.
Figure 4b compares BOSS against GRaSP on linear Gaussian data generated from scale-free networks with focus on scalability; runs exceeding two hours were cancelled and are not reported. Here we that BOSS maintains high level of accuracy while scaling much better than GRaSP. It is also
7

---

Adjacency
Orientation
Adjacency
Orientation
0-00-
Precision
Precision
0.5
0.5
0
Recall
Recall
0.5
0.5
0-
5
10
15
20
5
10
15
20
5
10
15
20
10
15
20
Average Degree
Average Degree
Average Degree
Average Degree
BOSS
GRaSP DAGMA
IGES
BOSS
GBaSP DAGHA
IGES
PC
(a) Erdös-Rényi
(b) Scale-free
Figure 3: Mean statistics over 20 repetitions: 100 variables and sample size 1,000.
Adjacency
Orientation
Precision
10
10a
Soconds
10a
Recall
10a
0.5
109
200
400
600
800
1,000
10
15
20
5
10
15
20
Measured Variables
Average Degree
Average Degree
BIC
4
Gumbel:
BOSS
DAGKA DAGKA
LINGAM LANGAM
BOSS
Exponential
BOSS
GRaSP
Con-Gaussian
(b) High Dimension
Figure 4: Mean statistics over 20 repetitions: scale-free, 100 variables, average degree 20 (when not varied), and sample size 1,000.
important to note that nearly all of these simulations were infeasible for GRaSP before implementing with GSTs. Table the corresponding table, only reports results for BOSS since the two algorithms have nearly identical performance all statistics except for running time. Full results are tabulated the Supplement.
Validation on fMRI
The high spatial coverage of fMRI's has allowed researchers to study brain function at different scales, from voxels to cortical parcellations to functional systems (such as default mode or visual systems). Given the potential differences in spatial dimension, the expected connectivity density of functional brain networks remains unclear. For example, even after keeping only the 75th percentile stronger connections, structural connectivity networks from cortical regions 23 had 1012 connections on average. Since structural connectivity supports functional connectivity, we could expect empirical

---

Table 2: Mean statistics over 10 repetitions: scale-free, average degree 20, and sample size 1,000.
Variables
Algorithm
BIC
Adj Pre
Adj Rec
Ori Pre
Ori Rec
2
0.98
0.80
0.97
0.80
500
BOSS
4
1.00
0.72
0.99
0.72
8
1.00
0.57
0.99
0.57
2
0.97
0.80
0.97
0.80
1000
BOSS
4
1.00
0.73
1.00
0.72
8
1.00
0.59
1.00
0.58
Table fMRI simulated data with pseudo-empirical errors
Algorithm
Adj Pre
Adj Rec
Ori Pre
Ori Rec
ABIC
Edges
Seconds
BOSS
0.99
0.94
0.96
0.90
211.79
951.40
15.46
fGES
0.97
0.60
0.70
0.43
7784.48
617.35
5.16
DAGMA
1.00
0.69
0.98
0.67
3080.85
687.15
54.58
LiNGAM
0.54
0.94
0.35
3868.93
1752.75
582.03
functional connectivity networks to be in that order of magnitude. Previous studies applying causal discovery methods to fMRI simulated networks with a high number variables and connections have shown that while the adjacency recovery precision of these methods can be very high, they usually show a low adjacency recall [16 20). Despite possibly having low recall, the limited applications of causal discovery methods to real world data [3] [8) often recover models with average degree greater than 20. Considering this limitation and the likelihood that real IMRI brain networks have high average connectivity (degree), causal discovery methods capable of reducing the number of false-negative connections will substantially improve future analysis fMRI data.
To demonstrate its practical utility in this important real-world domain, apply BOSS to two types of resting-stat state fMRI data: simulated data with pseudo-empirical noise distributions derived from randomized empirical fMRI cortical signals and clinical data from 3T fMRI scans processed into cortical parcels.
6.1 Simulated IMRI
We simulated fMRI data following the approach in [9]. Networks were based on a directed random graphical model that prefers common causes and causal chains over colliders. 40 networks were simulated with 200 variables and an average degree of 10. Edge weights were sampled from uniform distribution in the interval 0.4], randomly setting 10% of the coefficients to their negative value. Pseudo-empirical noise terms were produced by randomizing fMRI resting-state data across data points, regions, and participants from the Human Connectome Project (HCP). Using pseudo-empirical terms better captures the marginal distributional properties of the empirical fMRL. One thousand data points were generated from this procedure. These data are available for download at: htttt::/// / thub. com/ cmu pwwppphssls.e.ppe
Accuracy and timing results are shown in Table for BOSS, fGES, DAGMA, and LiNGAM. These results show that BOSS has by far the best BIC score of the group, compared to the BIC score of the true model, and that the running time of BOSS very reasonable for problem of this size. Precisions and recalls for BOSS are quite high, for both adjacencies and orientations. The poor performance of LiNGAM due to insufficient on-Gaussiia signal. More details are included in the Supplement.
Clinical fMRT
We applied BOSS to 171 3-Tesla resting state fMRI scans from patients beginning treatment for alcohol use disorder 3. Before applying BOSS, the data were cleaned and parcellated into 379 variables representing biologically interpretable and spatially contiguous regions 6 More details on the data collection, cleaning, and processing can found in (3]. 3 the in/out-degree distributions of the resulting models. The solid line depicts the median degree across all graphs and

---

0.6
0.6
Density
0.4
Density
0.4
0.2
0.2
0
0
10
20
30
20
40
60
In-Degree
Out-Degree
Figure 5: In/out-degree distributions on clinical resting state fMRI data.
the color region shades in the area between the 2.5 and 97.5 percentiles. These plots indicate that the models scale-free, with small number of hub vertices with much higher degree than other vertices. Scale-free connectivity is consistent with biological expectations and prior connectivity analysis on IMRI data [8]. More details are included in the Supplement.
7 Discussion
We have proposed a successor to the GRaSP algorithm 12 that retains its high accuracy but faster and more scalable. BOSS implemented with GSTs comfortably scales to least 1000 variables with an average degree of at least 20. It can comfortably and informatively analyze data from 400 or even 1000 densely connected IMRI cortical parcellations. We show in simulations that our method is highly accurate for Erdös-Rénn Rényi graphs well as scale-free graphs. Despite being developed for the linear Gaussian case, BOSS also performs well in the linear non- -aa-ssas case. BOSS is fast and accurate on simulated fMRI data and can rapidly produce informative and plausible models from clinical fMRI data. BOSS available for use within the TETRAD project which includes Python and R wrappers [17]
The success of BOSS presents opportunity pursue further theoretical work showing how BOSS differs from GRaSP and what general lessons we may learn for constructing successor algorithms to BOSS. There is substantial room for additional improvements, as BOSS does have several limitations. For example, BOSS cannot handle most forms of unmeasured confounding, so it will be valuable to explore ways of adding this functionality while maintaining its accuracy and scalability. will also be informative to apply BOSS to other types of data such as functional genomic data, financial data, and electronic health records.
Acknowledgments and Disclosure of Funding
We thank our anonymous reviewers for their detailed and insightful comments. BA was supported by the US National Institutes of Health under the Comorbidity: Substance Use Disorders and Other Psychiatric Conditions Training Program T32DA03718. 3. JR supported by the US Department of Defense under Contract Number FA8702- 15--00001---0022 with Carnegie Mellon University for the operation of the Software Engineering Institute. RSR supported by the US National Institutes of Health under awards RO1 RG05555556 and R01MH109520 JC was supported by the US National Institutes of Health under awards K01 AA026349 and ULITR002494, the UMN Medical Discovery Team on Addiction, and the Westlake Wells Foundation. EK was supported by the US National Institutes of Health under award UL ULITR002494, The content of this paper solely the responsibility of the authors and does not necessarily represent the official views of these funding agencies. All authors have conflicts of interest report.
10

---

References
[1] Barabási, A.-L. and Albert, R. (1999). Emergence of scaling in random networks. science, 286:509-512.
[2] Bello, K., Aragam, B., and Ravikumar, P. (2022). DAGMA: Learning DAGs via m-matrices and a log-determinant acyclicity characterization. In Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., and Oh, A., editors, Advances in Neural Information Processing Systems, volume 35, pages 5226-8239. Curran Associates, Inc.
[3] Camchong, J., Roediger, D., Fiecas, M., Gilmore, C. S., Kushner, M., Kummerfeld. E., Mueller, B. A., and Lim, K. O. (2023). Frontal tDCS reduces alcohol relapse rates by increasing connections from left dorsolateral prefrontal cortex to addiction networks. Brain Stimulation, 16: 1 1032-1040.
[4] Chickering, D. M. (1995). A transformational characterization of equivalent Bayesian network structures. In Proceedings of the Conference Uncertainty in Artificial Intelligence, page 87-98. Morgan Kaufmann Publishers Inc.
[5] Chickering, D. M. (2002). Optimal structure identification with greedy search. Journal of machine learning research, 3:507-554
[6] Glasser, M. F., Coalson, T. S., Robinson, E. C.,, Hacker, C. D., Harwell, J., Yacoub, E., Ugurbil, K., Andersson, J., Beckmann, C. F., Jenkinson, M., et al. (2016). multi-modal parcellation of human cerebral cortex. Nature, 536:171-178.
[7] Haughton, D. M. (1988). On the choice of model fit data from an exponential family. The annals of statistics, pages 342-355.
[8] Ikeuchi, T., Ide, M., Zeng, Y., Maeda, T. N., and Shimizu, S. (2023). Python package for causal discovery based on LiNGAM. Journal of Machine Learning Research, 24:1-8.
founmettiaaoounattooon of asymptotic inference. John Wiley [9] Kass, R. and Vos, W. (2011). Geometrical foundations & Sons.
[10] Kummerfeld, E., Williams, L.,, and Ma, S. (2023). Power analysis for causal discovery. Interna- tional Journal of Data Science and Analytics, pages 1-16.
[11] Lam, W.-y. (2023). Causal razors: A comparative study of simplicity assumptions in causal discovery. arXiv preprini 2rXiv:2302.1033/
[12] Lam, W.-Y., Andrews, B., and Ramsey, J. (2022). Greedy relaxations of the sparsest permutation algorithm. In Proceedings of the Conference on Uncertainty Artificial Intelligence, pages 1052- 1062. PMLR.
[13] Margaritis, D. and Thrun, S. (1999), Bayesian network induction via local neighborhoods. Advances in neural information processing systems, 12.
[14] Ng. I.. Ghassami, A., and Zhang, K. (2020). On the role of sparsity and DAG constraints for learning linear DAGs. Advances in Neural Information Processing Systems, 33: 19943-11779
[15] Ramsey, J., Glymour, M., Sánchez-Romers R., and Glymour, C. (2017). A million variables and more: the fast greedy equivalence search algorithm for learning nigh-dimensional graphical causal models, with an application to functional magnetic resonance images. International journal of data science and analytics, 3:121-129.
[16] Ramsey, J. D. (2021). Improving accuracy of permutation DAG search using best order score search. arXiv preprint arXiv:2108. 10141.
[17] Ramsey, J. D., Zhang, K., Glymour, M., Romero, R. S., Huang, B., Fbert-Uphoff, I., Samaras- inghe, S., Barnes, F. A., and Glymour, C. (2018). Tetrad a toolbox for causal discovery. In 8th international workshop on climate informatics.
[18] Rawls, E., Kummerfeld, E., Mueller, B. A., Ma, S., and Zilverstand, A. (2022). The resting- state causal human connectome is characterized by hub connectivity of executive and attentional networks. Neurolmage, 255: 19211.
11

---

[19] Sánchez-Romen R. and Cole, M. W. (2021). Combining multiple functional connectivity methods to improve causal inferences. Journal of cognitive neuroscience, 33:180-194.
[20] Sánchez-Romero, R., Ramsey, J. D., Zhang, K., Glymour, M. R., Huang, B., and Glymour,
C. (2019). Estimating feedforward and feedback effective connections from fmri time series: Assessments of statistical methods. Network Neuroscience, 3:274-306.
[21] Schwarz, G. (1978). Estimating the dimension of a model. The annals of statistics, pages
[22] Shimizu, S., Inazumi, T.,, Sogawa, Y., Hyvarinen, A., Kawahara, Y., Washio, Hoyer, P. O.. Bollen, K., and Hoyer, P. (2011). DirectLiNGAM: A direct method for learning a linear non- Gaussian structural equation model. Journal of Machine Learning Research-JMLR, 12: 12:1225-1248.
[23] Skoch, A., Rehák Bučková, B., Mareš, J., Tintéra, J., Sanda, P., Jajcay, L., Horáček, J., Spaniel, F., and Hlinka, J. (2022). Human brain structural connectivity matrices-ready for modelling. Scientific Data, 9:486.
[24] Solus, I.., Wang, Y., and Uhler, C. (2021). Consistency guarantees for greedy permutation-based causal inference algorithms. Biometrika, 108:795-811 814.
[25] Spirtes, P., Glymour, C., and Scheines, R. (2000). Causation, prediction, and search. MIT press.
[26] Teyssier, M. and Koller, D. (2005). Ordering- baseed search: simple and effective algorithm for learning Bayesian networks. In Proceedings of the Conference on Uncertainty in Artificial Intelligence, pages 584-590.
[27] Uhler, C.. Raskutti, G., Bühlmann, P., and Yu, (2013). Geometry of the faithfulness assumption in causal inference. The Annals of Statistics,
[28] Zhang, J. (2013). A comparison of three Occam's razors for Markovian causal models. The British Journal for the Philosophy of Science, 64:423-448.
[29] Zhang. J. and Spirtes, P. (2008). Detection unfaithfulness and robust causal inference. Minds and Machines, 18:239-271.
12


Reviews:
[{"review_id": "fqUM0h3MK5", "reviewer": "Reviewer_9NHG", "summary": "The paper proposes a computationally efficient algorithm using the grow-shrink trees to iterate through some of the combinatoric number of directed acyclic graphs in a graphical model. The approach builds on early work but is faster. ", "questions": "1. What is the computational complexity of the algorithm?\n2.  Proposition 4 states that the initial permutation doesn't matter if some conditions are met. In practice, if these conditions are unmet or unknown, does the initial permutation matter?\n", "limitations": "Limitations are discussed.", "rating": 5, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is straightforward, easy to read (examples are given), and focused in its presentation. The results show that it is faster than the baseline GRaSP. \n\n", "weaknesses": "The paper's focus on the approach also means there is not much perspective given until the final discussion. \n\nThe paper's claim of scaling to thousands is not thoroughly tested (results for 1000 appear to take 500 s). It is not clear what is necessary to apply the approach to higher-resolution parcellations. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**The paper's focus on the approach also means little perspective is given until the final discussion.**\n\nWe will add more details to / rework the introduction to help readers understand our perspective. In order to help us do this, could you describe what perspective you gained in the final discussion that was missing in the introduction?\n\n**The paper's claim of scaling to thousands is not thoroughly tested (results for 1000 appear to take 500 s). It is unclear what is necessary to apply the approach to higher-resolution parcellations.**\n\nOur intention was not to claim that our algorithm can scale to thousands of variables but rather to demonstrate that we can scale to at least 1000 variables on a laptop in under 10 minutes — this was done on a single processor and could be further sped up by using multiple processors in parallel. The paper presents an algorithm that can scale to densely connected problems an order of magnitude larger than what was previously possible in terms of variables (GRaSP without the use of GSTs does not scale much past 100 variables). As far as we know, fMRI parcellations in the literature rarely have more than 1000 variables.\n\n**What is the computational complexity of the algorithm?**\n\nWe will add a discussion about the algorithm’s complexity and running time. In the worst case, the main loop will be performed no more than n times, where n is the number of variables. Although, in our experiments, BOSS usually only repeated the main loop two or three times. Each main loop requires 3n^2 runs of the grow shrink subroutine, so BOSS will make O(n^3) calls to the GSTs. While profiling our algorithm, we found that it spends nearly all of its execution time performing such GST calls. These calls are cached so redundant regressions that need to be performed while running similar grow shrink subroutines can simply be looked up rather than recalculated. Indeed, one can show that once a regression is performed, it will never need to be run again.\n\n**Proposition 4 states that the initial permutation doesn't matter if some conditions are met. In practice, does the initial permutation matter if these conditions are unmet or unknown?**\n\nIn general, it is possible to run the entire algorithm multiple times from different initial permutations and take the result that yields the most optimal BIC score; our software has a parameter “number of runs” that allows the user to do this. In our simulations, when the ground truth is a DAG, there is little to no substantial advantage to consider multiple starting points. However, assuming that the ground truth is a DAG is admittedly a bad assumption for most real world datasets. In the causally insufficient case, regardless of the initial permutation, we can expect to learn a few extra edges. This can be corrected if the user post-processes the learned graph with the FCI algorithm, like what is done in the GFCI algorithm; we leave this investigation to future work.\n\nAlso, relative to confounding and fMRI, together with the whole-brain coverage, I would add something about the preprocessing that was used for this real fMRI data, in particular, what is commonly referred to as nuisance regression, which implies regressing the bold signal out of movements, trends, white matter, and csf artifacts, etc. Both strategies strongly mitigate the effect of possible confounders."}, {"review_id": "nKDXdI0ZxQ", "reviewer": "Reviewer_yEGg", "summary": "The authors introduce two novel methods for learning Directed Acyclic Graphs (ADGs): best order score search or BOSS and grow-shrink\ntrees or GSTs. These methods have a similar performance to state-of-the-art approaches (namely GRaSP and also others: fGES, DAGMA, LiNGAM, etc), but the authors proof that are less computationally intensive and therefore more scalable. The authors then validate their methods by learning brain networks from synthetic and real brain recordings (fMRI).", "questions": "No questions.", "limitations": "Authors have adequately addressed the limitations of their work. ", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- originality: the paper introduces two novel approaches for learning DAGs. Even though this work closely follows previous studies on permutation-based structure learning, it also presents significant new departures from the previous algorithms. \n\n- quality: the theoretical proofs seem solid and the methods are backed with a reasonable amount of numerical and experimental evidence. \n\n- clarity: the paper is very concisely and clearly written.  Although a lot of the information is in the supplements. \n\n- significance: learning network from brain recordings is a very useful tool to understanding the brain and its pathologies. A method that significantly increases scalability with same accuracy as the state-of-the-art is therefore noteworthy.  ", "weaknesses": "- Validation on real fMRI is (somewhat understandably) weak because the scale free connectivity of the brain remains a theory rather than a proven fact. \n\n- Many details of the work are left out from the main paper and included in the supplemental material. It is almost hard to understand the main paper without reading the supplements.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Validation on real fMRI is (somewhat understandably) weak because the scale-free connectivity of the brain remains a theory rather than a proven fact.**\n\nThank you for your comment, we will add references (a few relevant references are listed below) to bolster this particular point. These references give more recent, strong positive support for this position.\n\nLynn, C. W., & Bassett, D. S. (2019). The physics of brain network structure, function and control. Nature Reviews Physics, 1(5), 318-332.\n\nGalinsky, V. L., & Frank, L. R. (2017). A unified theory of neuro-MRI data shows scale-free nature of connectivity modes. Neural computation, 29(6), 1441-1467.\n\nMansour L, S., Di Biase, M. A., Smith, R. E., Zalesky, A., & Seguin, C. (2023). Connectomes for 40,000 UK Biobank participants: A multi-modal, multi-scale brain network resource. bioRxiv, 2023-03.\n\nHanson, S. J., Mastrovito, D., Hanson, C., Ramsey, J., & Glymour, C. (2016). Scale-free exponents of resting state provide a biomarker for typical and atypical brain activity. arXiv preprint arXiv:1605.09282.\n\nZhang, A., Fang, J., Liang, F., Calhoun, V. D., & Wang, Y. P. (2018). Aberrant brain connectivity in schizophrenia detected via a fast gaussian graphical model. IEEE journal of biomedical and health informatics, 23(4), 1479-1489.\n\nGrosu, G. F., Hopp, A. V., Moca, V. V., Bârzan, H., Ciuparu, A., Ercsey-Ravasz, M., ... & Mureșan, R. C. (2023). The fractal brain: scale-invariance in structure and dynamics. Cerebral Cortex, 33(8), 4574-4605.\n\n\n**Many details of the work are left out of the main paper and included in the supplemental material. It is almost hard to understand the main paper without reading the supplements.**\n\nThank you for this comment, the supplement contains details about: (1) how the data were simulated, (2) the parameters chosen for the various algorithms, and (3) tables giving the exact values of the results that are depicted as plots in the main paper. We included these details in the supplement for readers who wish to delve more into the details of our experiments, and we do not feel that the paper is incomplete without them. Could you be more precise about which aspect of the supplement you think should be included in the main paper?\n"}, {"review_id": "CLRvAtffq9", "reviewer": "Reviewer_P8PW", "summary": "It's beyond my knowledge so please disregard my scores and reviews.", "questions": "It's beyond my knowledge so please disregard my scores and reviews.", "limitations": "It's beyond my knowledge so please disregard my scores and reviews.", "rating": 5, "confidence": 1, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "It's beyond my knowledge so please disregard my scores and reviews.", "weaknesses": "It's beyond my knowledge so please disregard my scores and reviews.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your honesty!"}, {"review_id": "htnVk32ygc", "reviewer": "Reviewer_oAsg", "summary": "This paper introduces the best order score search (BOSS) algorithm and grow-shrink trees (GSTs) for learning directed acyclic graphs (DAGs) in machine learning and causal discovery. BOSS achieves state-of-the-art performance in accuracy and execution time, making it a valuable tool for problems with highly connected variables. The paper also applies BOSS to resting-state fMRI data, demonstrating its practicality and effectiveness. The algorithm is designed to be a more efficient alternative to the existing Greedy Sparsest Permutation (GRaSP) algorithm.", "questions": "- How does the BOSS algorithm handle the problem of confounding variables in fMRI data, and how does it ensure that the recovered causal networks are not biased by these variables?\n- Would there be non-linear causal relationships in fMRI data, and how does it ensure that the recovered causal networks capture these relationships accurately?\n- How does the performance of the BOSS algorithm vary with the level of sparsity in the underlying causal graph in fMRI data?\n- Would the performance of the model vary with the level of noise in data?\n- Discussion on Limitations: Could you include a section discussing the limitations of your work and potential areas for future research? \n- Results Interpretation: Could you provide more insight into how you interpreted your results and why you drew the conclusions you did? This would help to understand the implications of your findings.\n\n", "limitations": "- The authors have discussed the limitations of their work in the context of the performance of the BOSS algorithm. They have compared it with other algorithms like GRaSP, fGES, PC, and DAGMA on linear Gaussian data and have also tested it on non-Gaussian data. They have acknowledged that while BOSS performs well in terms of BIC score and running time, its recall may be low in certain scenarios. \nQuestion: You have mentioned that BOSS has a low recall in certain scenarios. Could you elaborate on these scenarios? How could this limitation be addressed in future work?\n- It would be beneficial to include a section in your paper discussing the ethical implications of your work. This could include potential misuse of your algorithm, biases in the data that could affect the algorithm's performance, and the implications of incorrect predictions by the algorithm.", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper introduces the concept of Grow-Shrink Trees (GSTs), which are tree data structures for caching the results of the grow and shrink subroutines of GS. This data structure is compatible with many permutation-based structure learning algorithms, including BOSS and GRaSP. The use of GSTs is a novel approach that efficiently stores information needed for running GS.\n- The paper is well-written and clear. It provides detailed explanations and examples, making it easy to understand the proposed algorithm and its components. The use of figures and tables also helps to illustrate the concepts and results.\n- The pape compares BOSS against other algorithms such as GRaSP, fGES, PC, and DAGMA on linear Gaussian data generated from Erd˝os-Rényi and scale-free networks. The results show that BOSS maintains a high level of accuracy while scaling much better than GRaSP, indcating that it could be a valuable tool in the field of structure learning.", "weaknesses": "- One potential weakness is the lack of discussion on the limitations of the BOSS algorithm and GSTs. While the paper acknowledges that there is room for additional improvements, it does not provide a detailed analysis of the limitations of the proposed approach. For example, the paper could discuss the assumptions made by BOSS, such as causal sufficiency, and how these assumptions may affect the accuracy and scalability of the algorithm. Additionally, the paper could explore the applicability of BOSS to other types of data beyond fMRI, such as financial data or electronic health records. Addressing these limitations would provide a more comprehensive understanding of the strengths and weaknesses of the proposed approach.\n- Validation on more Real-World Data: The validation of the algorithm on fMRI data is a good step, but more extensive validation on diverse real-world datasets could strengthen the robustness of the findings.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Apologies for truncating a few of your comments, we ran out of characters for our rebuttal.\n\n**One potential weakness is...**\n\nThanks for your comments and criticisms; we have found them helpful and constructive! See our responses to your points below.   \n\n**Validation on more Real-World Data: The validation of the algorithm on fMRI data is a good step, but more extensive validation on diverse real-world datasets could strengthen the robustness of the findings.**\n\nWe hope to apply this algorithm to many real-world datasets in the future; however, we lack the space to include multiple examples in the current paper. Regarding the choice of real-world application, fMRI was used due to the low recall of other methods in highly connected systems. The high recall of BOSS and GRaSP in these problems is not trivial, and GRaSP could not scale to the number of variables in most parcellations of fMRI data before this paper. We also think it is worth noting that the fMRI application included in the paper ran the BOSS algorithm on 171 datasets with 379 columns and 520 rows, so the included application was not a small undertaking.  \n\n\n**How does the BOSS algorithm handle the problem of confounding variables in fMRI data, and how does it ensure that the recovered causal networks are not biased by these variables?**\n\nThanks for the question; we refer you to our responses to the other reviewers on this point. In summary, we expect to find extra edges, but they can be removed using post-processing techniques which we leave for future work.\n\n**Would there be non-linear causal relationships in fMRI data, and how does it ensure that the recovered causal networks capture these relationships accurately?**\n\nIn general, fMRI data look fairly linear; however, it should be noted that the linearity assumption is made for the Gaussian BIC score and not for the BOSS algorithm itself. The Gaussian BIC score can be swapped out for another consistent score that models non-linearities in the data, and correctness will hold. However, it should be noted that such scores are often much slower, so the algorithm will not scale as well in this case. Lastly, we note that the dominant technique used to analyze function connectivity in fMRI data is to threshold the sample correlation matrix, which also does not handle non-linearity.  \n\n**How does the performance of the BOSS algorithm vary with the level of sparsity in the underlying causal graph in fMRI data?**\n\nIt is hard to analyze how BOSS deals with varying levels of sparsity in fMRI data since the true underlying structure is unknown. That being said, we believe the results that look at varying levels of sparsity for the scale-free linear cases should give an idea. Indeed, the method by which we simulated the scale-free aspect of the graph was inspired by the fMRI data. Figures 3 and 4 show the performance of BOSS (and other algorithms) as a function of the average degree of the underlying model. For the range of parameters we considered, the performance of BOSS (precision and recall for adjacencies and orientations) was essentially unaffected by variation in sparsity.\n\n\n**Would the model's performance vary with the noise level in the data?**\n\nIf the noise is increased in the exogenous noise terms, it does not appear to have a substantial effect; we’ve tested this for a wide range of these parameters (though there wasn’t room to include a formal result for this in the paper). If the noise is a large additional measurement noise, we have not tested this and plan to in future work. We have added this comment to the paper.  Small added Gaussian measurement noise has been unproblematic in our testing.\n\n\n**Discussion on Limitations: Could you include a section discussing the limitations of your work and potential areas for future research?**\n\nWe have collected up the various comments regarding the limitations of the work and potential areas of future research and will include them in a new section.\n\n\n**Results Interpretation: Could you provide more insight into how you interpreted your results and why you drew the conclusions you did? This would help to understand the implications of your findings.**\n\nThe results on the fMRI data is mainly included as a proof of concept. We hoped it would demonstrate that this algorithm can be applied to such data and the comments regarding the scale-free nature of the learned output are mainly to note that we are learning something that is consistent with the accepted nature of brain connectivity. We are currently working on several applied projects using this algorithm, where we plan to dive deeply into what we have discovered using this algorithm. This paper is intended to introduce the algorithm but not dive deeply into any data analysis; we leave that to future work.\n\n**The authors have discussed the limitations of their work in the context of...**\n\nThe performance of BOSS drops off in Table 2 (Figure 4b) when the penalty discount (BIC lambda) is increased. This was done for the sake of comparison to the GRaSP algorithm, which required the lower penalty discount in order to return in a reasonable amount of time (under an hour). However, there is no reason to not also report the results of the BOSS algorithm run on these data with a lower penalty discount. We expect the recall of BOSS to improve by doing so and plan to add these results to Table 2 (Figure 4b). The performance of BOSS could also degrade when its assumptions are violated, and we will include a discussion on this in the added section in limitations.\n  \n\n**It would be beneficial...**\n\nThanks for the suggestion, we will comment on the ethical implications of misusing / misinterpreting the results of our algorithm as you suggest with / after the added discussion of limitations."}]

Author rebuttal:
We thank all reviewers for your comments! We noticed that none of the reviewers were especially confident in their reviews (none had confidence greater than 2); perhaps our responses can help in this regard. 

The reviewers did not have many comments regarding the technicalities of the algorithm or data structure presented in the paper. There were, however, some areas where the presentation can be improved, and it’s clear from the comments that some context is needed to tell the reader what the contributions of the paper are to the practical issue of analyzing large-scale data accurately from a causal point of view.

We will address comments that ask for clarification or more information and add text in the paper to address these comments.

In general we plan to:

 - add a discussion of future work, limitations, and ethical considerations

 - add a discussion about complexity of the algorithm

 - add more details to the figure that demonstrates BOSS's scalability

 - clarify that GRaSP is using GSTs — without GSTs, GRaSP does not scale much past 100 variables

 - add comments about Direct LiNGAM in the Gaussian case

 - adjust language for clarity where necessary and correct any typos


Meta-review: The paper extends the Greedy Relaxation of the Sparsest Permutation (GRaSP) approach. Its contribution is to design and exploit a data structure and caching method to explore the permutation space and reuse the intermediate calculations, enabling causal discovery to scale up to some 1,000 nodes for problems with very large degree (20), which is a breakthrough in the causal discovery domain. 
 
The discussion with the reviewers showed that the paper needs a careful contextual setting about permutation-based search to appreciate the originality of the approach. Please also clarify its limitations, in particular w.r.t. the sample size.
Decision: Accept (poster)