ConKer: evaluating isotropic correlations of arbitrary order

High order correlations in the cosmic matter density have become increasingly valuable in cosmological analyses. However, computing such correlation functions is computationally expensive. We aim to circumvent these challenges by designing a new method of estimating correlation functions. This is realized in ConKer, an algorithm that performs FFT convolutions of matter distributions with spherical kernels. ConKer is applied to the CMASS sample of the SDSS DR12 galaxy survey and used to compute the isotropic correlation up to correlation order n = 5. We also compare the n = 2 and n = 3 cases to traditional algorithms to verify the accuracy of the new method. We perform a timing study of the algorithm and find that two of the three components of the algorithm are independent of the catalog size, N, while one component is O(N), which starts dominating for catalogs larger than 10M objects. For n<5 the dominant calculation is O(N^(4/3logN)), where N is the number of the grid cells. For higher n, the execution time is expected to be dominated by a component with time complexity O(N^((n+2)/3)). We find ConKer to be a fast and accurate method of probing high order correlations in the cosmic matter density.


Introduction
Understanding the dynamics of inflation in the early universe is tied with the study of primordial density fluctuations, in particular with their deviations from a Gaussian distribution (see e.g.Maldacena (2003), Bartolo et al. (2004), Acquaviva et al. (2003)).High order correlations have been shown to be sensitive to non-Gaussian density fluctuations (see Meerburg et al. (2019) and references within).Yet, a brute force approach leads to prohibitively expensive O(N n ) computation of correlations, where N is the number of objects and n is the correlation order.This problem has been studied extensively (March 2013).Several approaches to mitigate it were suggested for the calculation of three point correlations, each relying on a particular set of assumptions, e.g.small angle (Yuan et al. 2018), Legendre expansion (Slepian & Eisenstein 2015).The later approach was recently generalized for the n-point correlation function in Philcox et al. (2021) resulting in an O(N 2 ) algorithm.Here, we present an alternative computationally efficient way of evaluating such correlations.Similar to March (2013), this algorithm exploits spatial proximity and similar to Zhang & Yu (2011) and Slepian & Eisenstein (2016), it uses a Fast Fourier Transformation (FFT ).These characteristics combined with implementation facilities help achieve a notable reduction in calculation time without much memory overhead.The developed algorithm is applicable to discrete matter tracers, such as galaxies, as well as to continuous ones, such as Lyman-α and 21 cm line intensity, or matter maps derived from weak lenses.The method can be applied to evaluate self-correlations, as well as cross-correlations between different matter tracers.The algorithm Convolves Kernels with  (X, Y, Z), with the ρ(X, Y, Z) being the density of expected observations from tracers randomly distributed over with surveyed volume.We define the deviation from the average density as The two point correlation can be visualized as an excess of sticks of a given length over the ones that are randomly distributed over space.The three point correlation corresponds to an excess of triangles, the four point correlation -to counting pyramids (the four points do not necessarily lie in one plane), etc (see Fig 1).We will refer to these figures as n-pletes.
Let us consider all possible n-pletes with one vertex at point 0, characterized by a vector r, with sides radiating from this vertex equal to s 1 , s 2 , ...s (n−1) .In the isotropic case considered here the lengths of other sides and all angles are arbitrary.
The isotropic n-point correlation function (npcf) is defined as . ( where the integration over r implies all possible positions of point 0 in the surveyed volume.The integration over a solid angle dΩ i (i = 1, ...n − 1) between vectors s i and r is equivalent to integration over a sphere of radius s i centered on point 0, which, thus, evaluates the total number of tracers displaced from point 0 by the distance s i .This is illustrated in Fig. 2 for three point correlation.In term of Legendre polynomials expansion this definition of npcf corresponds to l = 0.
ConKer evaluates the integral of the ∆ field on a spherical shell of radius s i by placing a spherical kernel K i on point 0 with Cartesian coordinates (X, Y, Z) and taking its inner product with the ∆ field.Then the kernel is moved to a different location, thus scanning the entire surveyed volume.The results of each step, W i (X, Y, Z), is a map of weights that characterize the integral of the density deviation on a spherical shell of radius s i around the location (X, Y, Z).

Input
The inputs to the algorithm are catalogs of the observed number count of tracers D, with the total number being N tot , and R, which represents a number count of randomly distributed points within the same fiducial volume.Most surveys provide the angular coordinates: right ascension α and declination δ, and the redshift z of the tracer.The relationship between the redshift and the comoving radial distance is cosmology dependent: where Ω M , Ω k , and Ω Λ are the relative present day matter, curvature, and cosmological constant densities respectively.H 0 is the present day Hubble's constant, c is the speed of light.These user-defined parameters represent the fiducial cosmology.The integral in Eq. 3 is evaluated numerically in ConKer.The comoving Cartesian coordinates are computed with: ) Y = r cos(δ) sin(α) , (5) Z = r sin(δ) . (6) Cartesian coordinates are given as capital letters, X, Y, Z, as the redshift, z, is denoted in lower case.
We define a grid with spacing g s , such that the volume of each cubic grid cell is g 3 s .(X, Y, Z) denotes the Cartesian coordinates of the center of a grid cell.On this grid we define 3dimensional histograms D(X, Y, Z) and R(X, Y, Z), which represent tracer counts in the cell (X, Y, Z) from data and random catalogs respectively.These histograms may be populated by the raw count, or weighted count of tracers from the input catalogs.To properly account for the variation in angular completeness and the redshift selection function, the expected count from the random distribution, normalized to N tot , is subtracted from the galaxy survey D on the cell-by-cell basis.R is evaluated using a random catalog.As an alternative, R could be evaluated assuming that the tracers' angular, P ang (α, δ) and redshift2 , P z (z) probability density distributions are factorizable (Demina et al. 2018).Then the expected count from the random distribution of galaxies in a Cartesian cell is: where dXdYdZ and dV sph = r 2 d sin δdα give the volume of a particular cell in Cartesian and spherical sky coordinates.This method of estimating the expected count from the random distribution reduces the statistical fluctuations compared to the count taken directly from the random catalog.Both methods are realized in the algorithm and may be chosen by the user.A 3-dimensional local density variation histogram N(X, Y, Z), which is a discretized representation of the ∆ field, is then defined on the grid to represent the difference in counts between D and R: This stage of the algorithm is referred to as mapping.

Kernels
We construct (n − 1) spherical kernels, K i (X, Y, Z) on a cube of size 2s i in all directions, just large enough to encompass a sphere of radius s i .The grid defined on this cube has the same spacing, g s , as the one used to construct D(X, Y, Z).Grid cells intersected by a sphere of radius s i centered on the center of the cube, are assigned a value of 1.All other cells are 0.

Convolution
We construct a 3-dimensional histogram W i (X, Y, Z) with entries equal to the inner product of kernel K i centered on a cell with coordinates (X, Y, Z) and matter density variation field N: where * * * denotes a 3-dimensional discrete convolution.One step in this process is visualized in Fig. 3.This procedure is realized using an FFT in a highly vectorized algorithm.In this process, all N(X, Y, Z) cells are interpreted once as center points, i.e. point 0 in Fig. 2, and other times as points on shells.Thus the 3-dimensional histogram W 0 (X, Y, Z) is approximated as simply being the density variation in each cell of N(X, Y, Z), i.e.W 0 = N.
For normalization purposes, we perform the same procedure on the field of random counts.The result of the convolution of random density field R with kernel K i is referred to as B i (X, Y, Z): This stage of the algorithm is referred to as convolution.

n-point correlation function
The npcf, which approximates the one defined in Eq. 2 is calculated as: The 2pcf is then The 3pcf necessitates convolution with kernel K 2 of radius s 2 , resulting in weights W 2 .The 3pcf is then calculated as The stage of the evaluation of the npcf is referred to as summation.

Equidistant case
Of particular interest are correlations between objects with a defined scale, such as those arising from a spherical sound wave in the primordial plasma, a.k.a.baryon acoustic oscillation.In this case, one tracer taken as a starting point is displaced from the rest (n − 1) points by the same distance s 1 , e.g. the three point correlation corresponds to isosceles triangles randomly distributed in space.This equidistant case corresponds to the diagonal of the n-point correlation function, which is calculated as Note that only W 0 , W 1 , B 0 and B 1 need to be evaluated, resulting in a substantial time saving.

Continuous matter distribution
The continuous matter distribution map that results from line intensity or weak lensing, differs from the discrete case because the map does not have the meaning of counts.It can instead be interpreted as a δ-field, where δ has a meaning of relative overdensity: where D and R are continuous data and random fields defined on a grid, as is the δ-field, which is identical to W 0 map.The convolution of δ-field with kernels produces W i maps.Since δfield is already normalized by the random matter distribution it is not necessary to do it in the n-point correlation function, which is equal to: The two point correlation function is then

Cross-correlation between different tracers
The algorithm can be used to evaluate cross-correlation between different tracers.In this case we have two (possibly background subtracted) matter maps ∆ 0 and ∆ 1 .Then map of weights W 0 is equal to ∆ 0 , and K 1 is convolved with ∆ 1 resulting in W 1 .
The two point correlation function is then evaluated according to Eq. 17 or Eq. 12 depending on whether the tracer counts are normalized to the random distribution or not.

Performance study
We evaluate the performance of ConKer using SDSS DR12 CMASS galaxies (Ross et al. 2017), their associated random catalogs, and an ensemble of MultiDark-Patchy mocks (Kitaura et al. 2016;Rodríguez-Torres et al. 2016).We apply ConKer to the SGC and NGC catalogs for data, randoms, and 40 mocks.In this performance study, we used the following values for the cosmological parameters: c = 300000 km/s, H 0 = 100h km/s/Mpc, Ω M = 0.307, Ω Λ = 0.693, and Ω k = 0.This choice was motivated by the cosmological parameters used in the generation of the mock catalogs.For this study, which highlights the algorithm's ability to probe correlations near the BAO scale, we compute correlation functions for a distance range of 54.8-139.8h −1 Mpc in 17 bins of width of 5 h −1 Mpc.In all cases, the standard systematic (Ross et al. 2017) as well as FKP weights (Feldman et al. 1993) were used to create the density field map.The grid spacing for each run (with the exception of the timing study discussed in § 3.2) is g s = 5 h −1 Mpc.

Distance calibration
As in any algorithm that uses a map of tracers defined on a grid, there is a certain degree of degradation in the precision of the distance between the tracers.It is important to quantify the precision of distance determination.We calibrate the distance using a subsample of a random catalog, where we compare the distance measured on the grid with the one determined by a brute force calculation.Figure 4 shows the dependence of the mean true distance on the kernel size.The shaded region correspond to the RMS of the distribution in the true distance given a particular value of the kernel size, which is shown in the inset for a kernel size of 60 h −1 Mpc.The RMS is approximately 2.5 h −1 Mpc, which corresponds to one half grid spacing.The calibration procedure must be performed using a random catalog for each value of the grid spacing, g s .

Timing study
We present the analytical considerations of the time complexity of the three stages of the algorithm: mapping, convolution and summation in Appendix A.1.Here we show the results of the timing study performed using a personal computer with a 2.3 GHz Intel Core i7 processor and 8 GB of memory.All execution times are in units of CPU seconds.The primary advantage of ConKer is in the behavior of the execution time as a function of the number of objects, N, shown in top plot of Fig. 5 for the three stages of the algorithm.The surveyed volume, V is kept fixed.As expected the execution time of both convolution and summation are independent of N. Mapping is a O(N) calculation, and starts dominating for catalogs larger than 10M objects.
The main parameter that determines the computation time of ConKer is the number of grid cells, N c , which depends the grid spacing, g s (see the bottom plot of Fig 5 for n = 2, 3, 4).Mapping is independent of g s .For 2pcf convolution is expected to scale as g −4 s log g s , and summation as g −4 s .Both parts scale somewhat slower than expected.The execution time of convolution is independent of the correlations order n.Each additional order of cor-relation is expected to add one power of g s to summation time.This is roughly in agreement with the observed dependence for 3pcf and 4pcf.Summation starts dominating the execution time for 4pcf for grid spacing below 2 h −1 Mpc.For n > 4 the execution time is expected to be dominated by summation, which is a O(N (n+2)/3 c ) calculation.For lower n or higher grid spacing convolution is the dominant part, expected to be O(g −4 s log g s ), or O(N 4/3 c log N c ).The algorithm's memory requirements are discussed in Appendix A.2.

Comparison to nbodykit
We present a comparison of the 2pcf and 3pcf evaluated using ConKer and nbodykit (Hand 2018).The latter uses more traditional pair-counting methods, and may be used to evaluate the performance of ConKer.In nbodykit, the 2pcf is defined according to the estimator of Landy and Szalay (Landy & Szalay 1993;Hamilton 1993): where DD, RR, and DR are the normalized distributions of the pairwise combinations of galaxies from data, D, and random, R, catalogs (plus cross terms) at given distances s from one another.The nbodykit 3pcf is defined according to Slepian & Eisenstein (2015) and we consider the = 0 (isotropic) case.This corresponds to the ConKer 3pcf.
Using nbodykit, the 2pcf and 3pcf are calculated for the SGC sample of the SDSS DR12 CMASS galaxies and randoms.Due to the restriction that the 3pcf algorithm in nbodykit scales as O(N 2 ) (Slepian & Eisenstein 2015), where N is the number of objects, the random catalog was reduced to 1M objects.Full random catalogs were used when running ConKer.The data catalog remained unchanged in both cases.The ConKer 3pcf was completed in 350s and the nbodykit 3pcf in 9870s.
The comparison between the two methods of computing the 2pcf is given in Fig. 6.The two algorithms agree with one another across the entire range of s 1 , and the residual between the 2pcf using ConKer and nbodykit is around an order of magnitude lower than the the value of ξ 2 .The observed small differences are due to the precision of the distance determination on the grid in ConKer resulting in some bin-to-bin migrations, while nbodykit calculates the actual distance between the two points.
We repeat this procedure with the on and off diagonal elements of the 3pcf (Fig. 7).Here, both methods agree as well, as demonstrated by the residual between the two methods shown in the bottom right of Fig. 7. Some difference is due to the reduction in the size of the randoms catalog for the nbodykit calculation, which leads to larger statistical fluctuations.Overall, ConKer is in a good agreement with standard algorithms, and has the unique advantage of being able to probe higher order correlations as well, with relatively little increase in run-time.

ConKer npcf
We use ConKer to compute the diagonal elements of the npcf for n = 2, 3, 4, 5 and off diagonal elements of the 3pcf for the SGC and NGC catalogs.The combination of NGC and SGC catalogs is done by adding together terms in the numerator and denominator of Eq. 11.For example, the combined 3pcf is given by Fig. 6.The top panel shows the 2pcf, ξ 2 (s 1 ), for the CMASS galaxy sample of the SDSS DR12 survey computed using ConKer (blue) and the 2pcf computed using nbodykit (gold).The bottom panel gives the residual between the two, defined as ξ 2,ConKer − ξ 2,nbodykit (grey dashed line).The solid black line in the bottom panel shows ∆ξ 2 = 0.
The diagonal elements of the npcf's n = 2, 3, 4 and 5 are shown in Fig. 8.We observe the expected features of the npcf's based on both mock and data catalogs.These include an increase in magnitude at small scales, corresponding to galaxy clustering, and a well defined 'bump' at the BAO scale.The npcf's based on the random catalog fluctuate about ξ n = 0 at several orders of magnitude below the signal observed in mocks and data.

Algorithmic framework
The query for a pattern in matter distributions may prompt an employment of machine learning techniques.ConKer, being a spatial statistics algorithm offers an alternative to such approach that is fast and transparent.It exploits the fact that the full set of equidistant points from any given point makes a sphere, with its surface density being a direct measure of how spherically structured this subspace of points is.Aggregating and combining this measure over the whole space of N objects allows us to calculate the space's n point correlation function.
The algorithm exploits an intrinsic spatial proximity characteristic in the objective of querying for structures of negligible dimensions in a much bigger space.This spatial proximity factor leads to space partitioning algorithms targeting a nearest neighbor query approach, see for example the tree-based npcf algorithms in March (2013).However, ConKer uses this fact as a heuristic in limiting its query space immediately to only the defined separation for each point in the space.Note how this is in contrast with the former technique family.In a nearest neighbor approach to the npcf problem, the query space is grown at each point in the greater embedding space, aggregating the n-point statistic until the greater space is fully queried, whereas ConKer aggregates the statistic over all embedding space, and then grows the query space before repeating.This design choice realized in ConKer's core subroutine, a convolution of the query space with the embedding space performed by an FFT algorithm, distills the complexity from the brute force O(N 2 c ) to O(N c logN c ), where N c is the total number of grid cells in a discrete 3-dimensional map.This approach, combined with the heuristic above, lets the dominant components of ConKer achieve independence of the number of objects as shown in Fig. 5.There is certainly a trade-off between the sparsity of the whole space and the bias towards linear complexity in number of objects, as expected from an FFT -based algorithm, but even for very dense catalogs, we expect the scaling in the number of objects to be capped by O(N), where N is the total number of objects.Ultimately, ConKer is a hybrid algorithm that draws from both computational geometry and signal processing to achieve linear complexity in the number of objects.

Applications beyond correlation functions
Though traditionally the order of correlation n is viewed as an important parameter, for the equidistant npcf, all the information is entirely encoded by the weights W 0 and W 1 .It was pointed out in Carron & Neyrinck (2012) that the npcf is inadequate in capturing the tails of non-Gaussianities.The distributions over W 0 and W 1 (as opposed to their sum over the sample as is used in the npcf) could recover that sensitivity, which is a subject for future studies.The distribution of the product of two weights, W 0 W 1 normalized by the average B 0 B 1 is presented in Fig. 10 for the kernel size of 107.3 h −1 Mpc for data, mock and random catalogs.
Moreover, W 0 and W 1 as well as their product are maps.While in the npcf, the location information is entirely lost, in the maps produced by ConKer, it is preserved and can be used for cross correlation studies between different tracers, such as weak lensing, Lyα and CMB.The preservation of this spatial information is illustrated by Fig. 11, where the DR12 SGC galaxies are plotted on top of the values of W 0 W 1 / B 0 B 1 for a kernel size ap-proximately equal to the BAO scale.The figure corresponds to a slice of thickness 5 h −1 Mpc.The observed circles around the boundaries of the catalog or the sparsely populated regions are a feature of the algorithm, since a single galaxy contributes to W 1 everywhere at a distance s from the galaxy forming a sphere.The highest values in the map (yellow color scale) correspond to the regions that are both densely populated and have several galaxies displaced by a given scale from this point, as one would expect from the centers of the baryon acoustic oscillations.Similarly, lowest values (blue color scale) correspond to under-dense regions and troughs of BAO originating from this location.Such an aggregated approach gives access to the matter distribution in the early universe.

Conclusion
We presented ConKer, an algorithm that convolves spherical kernels with matter maps allowing for fast evaluation of the isotropic n-point correlation functions.The algorithm can be broken into three stages: mapping, convolution and summation.Execution time of convolution and summation are independent of the catalog size, N, while mapping is a O(N) calculation, which starts dominating for catalogs larger than 10M objects.For the correlation orders n < 5, convolution is the dominant part with complexity O(N 4/3 c log N c ), where N c is the number of the grid cells.For higher n the overall time complexity is dominated by summation, which is O(N (n+2)/3 c ).
Comparison to the standard techniques shows a good agreement.We study the performance using SDSS DR12 CMASS galaxies, their associated random catalogs, and an ensemble of MultiDark-Patchy mocks.The results up to n = 5 are presented.Further metrics that may offer additional sensitivity to primordial non-Gaussianities are also suggested (e.g.distribution over weights W i ).

Fig. 2 .
Fig.2.A cartoon illustrating three point correlation (0 − 1 − 2) with different scales s 1 on the side (0−1) and s 2 on the side (0−2).Integration over spherical shells K 1 and K 2 is equivalent to counting all possible triangles that have one vertex at point 0 and the other two anywhere on the spheres.

Fig. 3 .
Fig.3.A 2-dimensional representation of one step in the convolution between the matter histogram N (left) and the spherical kernel K (center).In this example, the inner product is 7, which is placed in the center cell of W 1 histogram (right).

Fig. 4 .
Fig. 4. The true distances between the tracers and the center of the kernel vs the discretized estimation of the distance in ConKer (gold).The shaded region gives the RMS of the true distance distribution.The inset shows one such distribution for a kernel with size s 1 = 60 h −1 Mpc.

Fig. 5 .
Fig. 5. Color inline.The execution time of ConKer as applied to subsamples of the SDSS DR12 CMASS SGC galaxies.Top: the dependence on the number of combined data and random objects, using a grid spacing g S = 5 h −1 Mpc.Bottom: the dependence on the grid spacing, g S , for 0.87M objects.Points present the measured time.Lines are the results of the fit to a power law, with the scaling given in the legend.

Fig. 7 .
Fig. 7.The top left panel shows the 3pcf, ξ 3 (s 1 , s 2 ), for the CMASS galaxy sample of the SDSS DR12 survey computed using ConKer and the bottom left panel shows the 3pcf computed using nbodykit.The top right panel gives the residual between the two, defined similarly as ξ 3,ConKer − ξ 3,nbodykit .The color scale corresponds to all three color panels.The bottom right panel gives a distribution over the residuals.The RMS of this distribution is 1.34 × 10 −4 .

Fig. 8 .
Fig.8.The diagonal npcf, ξ n (s 1 ), from n = 2 to n = 5 (top to bottom panels), for the combined SGC+NGC CMASS galaxy sample of the SDSS DR12 survey (black), the average of the ensemble of MD Patchy mocks (green) and the associated randoms (red).The green shaded region gives the standard deviation across mock catalogs, which is also used to generate the error bars associated with the data catalogs.

Fig. 9 .
Fig. 9.The 3pcf calculated based on data (top left), random (bottom left), and mock (top right) catalogs, using the combined SGC and NGC samples.The panel in the bottom right shows the standard deviation across the mock ensemble.The color scale corresponds to all four panels.

Fig. 10 .
Fig. 10.The distribution, W 0 W 1 / B 0 B 1 at s 1 = 107.3h −1 Mpc, for the CMASS galaxy sample of the SDSS DR12 survey (black), one of the MD Patchy mocks (green) and the associated randoms (red).This plot includes SGC galaxies only.

Fig. 11 .
Fig. 11.A 5 h −1 Mpc thick slice of the SDSS DR12 SGC volume, where the color scale shows the value of W 0 W 1 / B 0 B 1 at s 1 = 107.3h −1 Mpc, and the red stars indicate the positions of galaxies.