Thanks for the comments! Let me elaborate a bit...
Nice--I'm one of the other authors; the main author gave me permission to respond, haha, he will regret it... :-)
Mentioning scale-freeness in the intro is a great idea. There is reasonable consensus among many now that fMRI connectivity is scale-free. GRaSP and BOSS deal with this sort of data very nicely, far better than the usual causal search algorithms, even out to an average degree of 20, without any compromise in quality. Also, we plan to rework the introduction in response to various comments.
We're thinking maybe the "wow" factor of doing graphs with an average degree of 20 perhaps didn't come through well. Nearly all papers we've seen in Neurips and elsewhere deal only with very sparse graphs, average degree 2 or 4, or 6, which for a large number of variables is sparse in the extreme. We wanted to emphasize that our method could handle a much higher average degree without sacrificing accuracy (as almost all algorithms do, to the best of our knowledge) and responds well to many people's worries about doing causal searches accurately with dense graphs. DirectLiNGAM does reasonably well for dense graphs, but only under the stronger assumption of linear (strong) non-Gaussianity, and the empirical fMRI data we use are quite Gaussian, so it is not surprising (see Table 2) in our paper that DirectLiNGAM performs at chance levels for this sort of data. (Compare this to DAGMA in our charts, for instance, another Gaussian method published recently in Neurips and the best-performing continuous optimization algorithm we could find, so this is a comparison to that whole approach, so far as we currently know.)
I was actually very excited to see the comment about scalability. For this paper, the state of the art for Gaussian searches for dense models is currently GRasP, though the published GRaSP scales only to about 100 variables, and in that paper, the maximum average degree tested was 10. In this paper, using a different permutation algorithm, we've increased the number of variables we can comfortably analyze 10-fold and the average degree we can analyze 2-fold without sacrificing accuracy. This sets a new state of the art for this problem, so far as we know. Notice also that even DirectLiNGAM, under a stronger assumption, takes much longer to run on this problem without appreciably better results; BOSS is able to match and even sometimes improve on the performance of DirectLiNGAM even for the linear, non-Gaussian case, which is helpful since even if you know that some variables are non-Gaussian you may not know that all variables are non-Gaussian.
What excited me was that you wanted to extend the analysis to voxel-level analysis, which would to at least 40,000 variables. Genetic data also scales into that general dimension if one includes protein data in the analysis. We are keen to scale methods to that dimension for dense, scale-free data and are currently working through the various issues with parallelization, etc. However, there is no room in the current paper to address those optimizations; we are at the page limit. We plan to write another paper that does. It goes without saying that with good parallelization (not easy, key steps not embarrassingly parallel), we could jump from doing these analyses on laptops in a single thread to many-core machines.
Also, we will definitely include the multiple-runs option in the methods section, thanks, and mention that for small numbers of variables (not so much for massive models), it can help. The idea of post-processing the BOSS output with a GFCI wrapper was proposed in the GRaSP paper but not explored, and we are hesitant to include any specific proposals on that before we've worked through the thorny issues of scaling a latent variable search up to the size problems that we are interested in. No one can do that now, even remotely, but we will try. We do have implementations of the procedure that people can try if they wish to see how well it does before such optimization. Using BOSS as an initial step certainly improves accuracy for these PAG methods.
Regarding not being abreast of the permutation search literature, for one thing, GSTs are completely novel as of this paper, as is the BOSS algorithm, so we weren't expecting anyone to be familiar with those, no worries. The permutation methods, though, are very promising, as they currently outperform DAG (CPDAG) search over all other methods we know of for the linear Gaussian case by a long shot and match the best-performing algorithms for the linear non-Gaussian case. We're hoping to raise awareness of this methodology--i.e., hoping people will go back and read the Raskutti and Uhler paper, the Solus et al. paper, and the Lam et al. GRaSP paper and various other papers that have been published using permutations of variables orders as a key step.