Computing 3 point correlation function randoms counts without the randoms catalogue

As we move towards future galaxy surveys, the three point statistics will be increasingly leveraged to enhance the constraining power of the data on cosmological parameters. An essential part of the three point function estimation is performing triplet counts of synthetic data points in random catalogues. Since triplet counting algorithms scale at best as $\mathcal{O}(N^2\log N)$ with the number of particles and the random catalogues are typically at least 50 times denser than the data, this tends to be by far the most time consuming part of the measurements. Here we present a simple method of computing the necessary triplet counts involving uniform random distributions by means of simple one dimensional integrals. This speeds up the computation of the three point function by orders of magnitude, eliminating the need for random catalogues, with the triplet -- and pair -- counting of the data points alone being sufficient.


INTRODUCTION
Current and past galaxy redshifts surveys have heavily relied upon the analysis of two point statistics in order to constrain cosmological parameters down to the percent level (see e.g. Anderson et al. 2012;Chuang et al. 2016;Cuesta et al. 2016;Gil-Marín et al. 2016a,b;Alam et al. 2017). As upcoming galaxy redshift surveys seek to push constraints on cosmological parameters to the sub-percent level, three points statistics -the three point correlation function (3PCF) and the bispectrum -will begin to play a bigger role in analyses. Recent works have shown that the baryon acoustic oscillation features can be detected in both the 3PCF (Slepian & Eisenstein 2015;Slepian et al. 2017a,b) and the bispectrum (Pearson & Samushia 2018, hinting at the possibility of increased constraining power on cosmological parameters via the inclusion of three point statistics. The large scale bispectrum has also been used to supplement the redshift-space distortion measurements from the power spectrum (Gil-Marín et al. 2015a,b, 2017. Turning to small scales, the Halo Occupation Distribution (HOD) model (Scherrer & Bertschinger 1991;Jing et al. 1998;Peacock & Smith 2000;Seljak 2000;Scoccimarro et al. 2001;Berlind & Weinberg 2002;Zheng et al. 2009, and references therein) is a popular ⋆ E-mail: dpearson@phys.ksu.edu method of linking the galaxy and dark matter distributions (see e.g. Tinker et al. 2005;Guo et al. 2014Guo et al. , 2016Rodríguez-Torres et al. 2016). Yuan, Eisenstein & Garrison (2018) have shown that the squeezed 3PCF can help tightly constrain the parameters of the HOD, making it very likely that the 3 point statistics will become important to the development of mock catalogs.
To measure the 3PCF from data, one must count the triangles of specific shapes and sizes from the data -e.g. the galaxy catalogue -as well as a set of unclustered random points -the random catalogue -a process which is naively O(N 3 ) in time complexity. The counts from the data are then compared to the expected mean numbers from the unclustered random points to get an estimate of the 3PCF (Peebles & Groth 1975;Peebles 1980), as the unclustered mean triplet counts are very sensitive to the geometry of the survey and its number density variations. The most popular estimator used for the 3PCF is that of Szapudi & Szalay (1998) due to its superior edge effect corrections (Kayo et al. 2004), [DDD(r 1 , r 2 , r 3 ) − 3DDR(r 1 , r 2 , r 3 ) where the combinations of Ds and Rs tell you how many vertices of the triangle come from either the data or the random points, respectively. The random catalogue usually contains 50 or more times as many objects as the data to make the shot-noise in this Monte-Carlo estimation subdominant to the variance in the data. Combine this with the O(N 3 ) complexity and the counts for DDR, DRR, and RRR can consume a large amount of time.
While the computational complexity of the three point statistics can be somewhat mitigated (Baldauf et al. 2015;Slepian & Eisenstein 2015;Pearson & Samushia 2018), their calculation still tends to take a significant amount of CPU time. Additionally, the studies which stand to benefit from the inclusion of the three point statistics may need to measure them tens of thousands of times, making any potential reduction in the computational complexity a welcome improvement.
In this letter we present a simple method of obtaining the counts involving the random points that does not require actually counting triangles from random catalogues. Since triplet-counting from random catalogues takes most of the computational time, this results in a significant reduction of required CPU hours. The method is restricted to measurements of 3PCF from uniform periodic cubes, which means that, unfortunately, it cannot be applied to the measurements from real survey data. There are, however, many stages in the cosmological analysis of the 3PCF that can be performed on uniform periodic cubes and do not need to account for survey geometry. These include the validation of theoretical templates -comparing model predictions to measurements from simulations and quantifying biasesand the HOD parameter fitting. For these applications our method allows one to set the number density of randoms arbitrarily high with no performance degradation.

Pair counts
We start by reminding the reader why random pair counting is not really necessary when computing the two-point statistics from a preiodic cube. The number of unclustered pairs separated by a distance r ±∆r/2 can be easily computed analytically. There are on average n R V box particles in a uniform periodic cube with volume V box and a number density of points n R . The volume between two concentric spheres of radii r − ∆r/2 and r + ∆r/2 -i.e. a spherical shell -is and the total number of pairs is since The cube is periodic there are no edge effects. We will now generalize this idea for triplet counts.

RRR counts
We need to estimate number of triplets separated by (r 1 ± ∆r/2, r 2 ± ∆r/2, r 3 ± ∆r/2) in a uniform periodic cube with a number density of n R and volume V box . We will start by computing this number for (r 1 , r 2 ± ∆r/2, r 3 ± ∆r/2). This can be achieved by finding the volume of the two overlapping In essence, we want to predict how many points, p 3 , there will be, which is simply the volume of the overlap region of the two spherical shells with central radii r 2 and r 3 , multiplied by the number density.
spherical shells whose cross-sections are show in Figure 1.
Multiplying that by the number density will then tell us the average number of p 3 points given p 1 and p 2 . We will make use of the equation for the volume of overlap of two spheres with radii r and R separated by a distance d (Weisstein 1999), This equation is valid only from the point where the two spheres just touch, up to the point where the smaller sphere completely falls within, but still touches the larger sphere, i.e. R − r d R + r. Because of this, we need to explicitly define the volume outside of those bounds, From careful study of Figure 1, we can find that the volume of interest is given by first finding the overlap volume of the two outer spherical surfaces, then subtracting the overlap volume of the outer spherical surface of one shell with the inner spherical surface of the other and vice versa. However, this ends up removing the overlap volume of the two inner spherical surfaces twice, so we have to add one back. Mathematically, this can be expressed as V cs (r 1 , r 2 , r 3 , ∆r) = V(r 1 , r 2 + ∆r/2, r 3 + ∆r/2) − V(r 1 , r 2 + ∆r/2, r 3 − ∆r/2) − V(r 1 , r 2 − ∆r/2, r 3 + ∆r/2) a.) b.) Figure 2. Examples of the r 2 -r 3 spherical shell overlap region in two special cases which are automatically accounted for by equation (5) but would not work with (4). Panel a shows a case where r 1 is less than ∆r. Panel b shows the cross-section for nearly co-linear triangles.
Here, care must be taken that the r of equation (5) is actually the smaller of the two spherical shell radii. For example, when considering isosceles or equilateral triangles, it is possible that r 2 +∆r/2 is larger than r 3 −∆r/2, since r 2 may equal r 3 . Exercising this caution, along with the special considerations of equation (5) will yield the correct volume even in special cases where the overlap volume looks quite different than in Figure 1 -see e.g. Figure 2. This volume times n R is the number of p 3 points falling into the overlap region. For a finite bin width, p 2 can be anywhere inside a spherical shell r 1 ± ∆r/2 around p 1 so we will have to integrate equation (6) with respect to 4πn R r 2 1 dr 1 . Finally, we have to account for the fact that there are on average n R V box points p 1 . Combining all of this results in an exact expression for the expected number of RRR counts without any shot-noise, expressed as a simple one dimensional integral, Here s is the number of unique permutations of the side lengths, with s = 1, 3, and 6 for equilateral, isosceles, and general triangles, respectively. We note that for most triangles this reduces to a surprisingly simple expression, However, this expression will break down in special cases such as those shown in Figure 2. For this reason, we recommend simply evaluating the integral in equation (7) numerically. 1

DRR counts
For DRR counts we have to replace the first vertex by a data point. Since data and random points are uncorrelated 2 between each other the formula for the DRR counts coincides 1 The terms in equation (6) will at most be polynomials of degree 5 in r ′ 1 , meaning a simple 3 point Gaussian quadrature rule is all that is needed. 2 This would not be the case for nontrivial survey geometries where the exact placement of data points with respect to the boundary makes a difference.
with the expression for RRR in equation (7) with one minor change due to the potentially differing number densities of data and random points where n D is the average number density of data points -i.e. n D = N D /V box .

DDR counts
The predictions of DDR(r 1 , r 2 , r 3 ) are almost identical to the procedure outlined in section 2.2, except for the last step where we integrate over 4πn R r 2 1 dr 1 . Since the data points are clustered, the distribution of p 2 around p 1 is not uniform and the number density depends on the distance r 1 resulting in Here n D (r) is the nonuniform data number density, which can easily be found from the data pair counts by first computing the number density in each spherical shell and then linearly interpolating 3 to the value at the particular r needed for the numerical integration. Note that we can no longer use a simple integer multiple to account for the unique permutations of the side lengths as the value of n D (r) will change as we permute. Given this, there will be 1, 3 or 6 terms between the braces for equilateral, isosceles or general triangles, respectively. We note that when summing these terms it is necessary to store the number of triangles as a floating point number to ensure accuracy. After the summing of all the necessary terms in the braces, the number can be rounded to an integer if desired.
Equations (7), (9), and (10) express all terms in the estimator of equation (1) that involve random points without any shot-noise and in terms of simple one dimensional integrals. DDR counts do require the pair counts of data points but these can be precomputed along with the DDD counts with no additional computational time. We have publicly released our code for doing the predictive calculations for use by the community. 4 . Differences in the average actual and predicted counts divided by their standard deviations added in quadrature for the cases of 1× and 10× randoms. Note the difference in vertical scale between the left and right hand panels. The triangle index range was truncated from 3128 to 1452, though we note that the remainder of the range looks very similar to what is seen after a triangle index of ∼600. It is clear that the RRR (top) and DRR (middle) predictions agree quiet well with the predictions for both the 1× and 10× cases. The DDR (bottom) counts agree quiet well for the vast majority of triangle indices, though there are notable deviations at lower triangle indices. We note that these effects are not systematic and decrease with increased number density of randoms, as well as with increased grid resolution in the generation of the lognormal mocks.

RESULTS
To test equations (7), (9) and (10) we compared their predictions with actual triplet counts. We generated 100 lognormal mock catalogues (Coles & Jones 1991) -N D ∼ 500000, V box = (1024 Mpc) 3 -to serve as a proxy for data points 5 and 100 random distributions with the same number density, and 100 with 10× the number density of data. We then performed the direct counts in bins of width ∆r = 1 Mpc, from r min = 0 Mpc to r max = 32 Mpc, and calculated our predictions. We then averaged together all 100 counts, and separately all of the predictions in order to make our comparisons.
To perform the actual counts we used a relatively straightforward, GPU accelerated, O(N 3 ) counting algorithm with periodic boundary conditions to remove any edge effects along with a few simple optimizations. 6 The first optimization employed is to bin the particles to a grid with spacing equal to r max , the maximum separation considered 5 While lognormal mocks do not adequately reproduce the threepoint clustering statistics observed in simulations or the real Universe, for the purposes of these tests we simply needed data with some clustering and are not concerned with the recovered clustering statistics. As we already had many lognormal mocks at our disposal, we utilized some of them out of convenience.

6
While we do not release our exact code used here publicly, we utilized the same algorithm in a library for use in a different project. For those who may be interested, you can view the GPU implementation at https://github.com/dpearson1983/ganpcf/blob/master/source/ganpcf.cu.
for any of the sides of the triangle. This allows us to check for triangles with points that are within r max of our first vertex by only checking points in the grid cell containing the vertex, and those immediately adjacent to that cell. The other, more minor optimization is to simply check that the distance between the first and second vertex is less than r max before checking for the third vertex, and if not, proceed to the next point for consideration as the second vertex. While it would usually be faster to not repeat count the same triangle by simply relabeling which is the first, second and third vertex, due to the peculiarities of GPU memory accesses, we found that it was actually more efficient to do the repeat counting. Figure 3 shows the difference of the average RRR (top), DRR (middle) and DDR (bottom) counts with the predictions from equations (7), (9) and (10) divided by their standard deviations added in quadrature, with the 1× and 10× randoms cases on the left and right, respectively. For clarity, the triangle index range was truncated. We note that in all cases the remainder of the range looks nearly identical to what is seen after a triangle index of ∼600. The predictions match the actual counts remarkably well, especially as the number of random points increases.
The DDR counts tend to show larger deviations between the direct counts and our predictions at small triangle indices -i.e. small scales. We verified that these deviations are not systematic -i.e. they will fluctuate up and down with equal probability for independent realisations averaging to zero -and that they tend to zero in the limit of high number density. They look systematic because at small triangular indexes the measurements share the same small number  (top) show that the grid spacing is partly to blame for the deviations. The Gaussian haloes (bottom) do not use a grid for their creation, and have an exactly known number density profile, instead of one estimated from the data-data pair counts. We see that there are no seemingly systematic deviations at small triangle indices. We note that there are some slight deviations at high indices, particularly if a finer binning is used, suggesting that the deviations are have more to do with the discrete nature random catalogues and grid effects.
of points and are strongly correlated -or anti-correlated depending on the shapes. Additionally, the deviations all occur when one side of the triangle is smaller than 5 Mpc, or 3 Mpc for the 1× and 10× randoms cases, respectively. The lognormal mocks were created using a grid spacing of 2 Mpc, where a number of galaxies was determined for each grid cell, then placed inside that cell in a uniform random fashion. We ran tests on mocks with a grid spacing of 0.5 Mpc but the same number density, and note that the deviations were reduced -see Figure 4. To further test that these small scale deviations were not due to a failing of the algorithm, but rather some combination of low number densities and grid effects in the mocks, we generated haloes with an exactly known, Gaussian number density profile, eliminating the need to use the data-data pair counts, making the predictions from equation (10) far more accurate. The random catalogues where set to the same number density -e.g. N D ∼ N R . We then only counted triangles which used the central particle as one of the vertices. The number density profile was such that the maximum particle radius was ∼25 Mpc with very few particles in the outer parts of the haloes. This led to seeing a similar -though significantly smaller -seemingly systematic effect for high triangle indices.
The results of the using a higher resolution grid for creating the lognormal mocks and the Gaussian haloes are shown in Figure 4. These results suggest that the seemingly systematic deviations in the DDR counts versus predictions have more to do with the discrete nature of the random catalogues, and grid based effects than a failing of the predictions.

CONCLUSIONS
We present a method for predicting the counts involving a random catalogue for the 3PCF analysis of simulated or mock data that is free of shot-noise and does not actually require random catalogues. We have shown that the predictions from equations (7), (9), and (10) agree remarkably well with the actual counts, while keeping the same computational complexity for arbitrarily high number densities of randoms.
The method only works for uniform periodic cubes which may lead the reader to believe that the method is not useful for the analysis of real survey data. This, however, is not true as there are many stages in the survey data analysis that require computing the 3PCF from a large number of periodic cubes.
The most obvious example is the validation of theoretical templates on N-body simulations. To make sure that theoretical predictions of the 3PCF are sufficiently accurate -and to calibrate any systematic effects in case they are not -they will have to be compared to the measurements from N-body simulations with a known cosmology, or a hidden cosmology in the case of data challenges. For these validation tests to be meaningful they need to be performed on periodic cubes to cleanly separate theoretical systematics from other observational effects. Since the cumulative volume of N-body simulations used for this purpose needs to be much larger than the size of the data, the method presented in this paper can save a significant amount of computational time.