Efficient parallel algorithm for estimating higher-order polyspectra

TL;DR

This paper introduces a highly efficient parallel algorithm for estimating higher-order polyspectra in galaxy density fields, leveraging FFT and optimized memory layout to handle large data with minimal communication overhead.

Contribution

The authors develop a novel parallel algorithm based on the Scoccimarro estimator that significantly improves the efficiency of measuring higher-order polyspectra in large datasets.

Findings

Reduces computational time for polyspectra estimation

Minimizes inter-CPU communication in parallel processing

Handles large datasets efficiently with optimized memory layout

Abstract

Nonlinearities in the gravitational evolution, galaxy bias, and redshift-space distortion drive the observed galaxy density fields away from the initial near-Gaussian states. Exploiting such a non-Gaussian galaxy density field requires measuring higher-order correlation functions, or, its Fourier counterpart, polyspectra. Here, we present an efficient parallel algorithm for estimating higher-order polyspectra. Based upon the Scoccimarro estimator, the estimator avoids direct sampling of polygons by using the Fast-Fourier Transform (FFT), and the parallelization overcomes the large memory requirement of the original estimator. In particular, we design the memory layout to minimize the inter-CPU communications, which excels in the code performance.

Figures

Equations

$$\delta (\mathbf{k})=$$

$$\delta (\mathbf{x})=$$

$$⟨\delta ({\mathbf{k}}_1)\cdots \delta ({\mathbf{k}}_n)⟩=(2\pi {)}^3{P}_n({\mathbf{k}}_1,\cdots ,{\mathbf{k}}_n){\delta }^D({\mathbf{k}}_{1\cdots n})$$

$$B(k_1,k_2,k_3)=$$

$$\times \delta ({\mathbf{q}}_1)\delta ({\mathbf{q}}_2)\delta ({\mathbf{q}}_3){\delta }^D({\mathbf{q}}_{123}),$$

$$B(k_1,k_2,k_3)=$$

$$\times \delta ({\mathbf{q}}_1)\delta ({\mathbf{q}}_2)\delta ({\mathbf{q}}_3){\delta }_D({\mathbf{q}}_{123})$$

$$V_{(123)}^B$$

$${\delta }_D({\mathbf{q}}_{123})=\int \frac{d^{3}x}{(2\pi {)}^3}{e}^{-i\mathbf{x}\cdot ({\mathbf{q}}_1+{\mathbf{q}}_2+{\mathbf{q}}_3)},$$

$$B(k_1,k_2,k_3)=\frac{V_f}{V_{(123)}^B(2\pi {)}^3}\int \frac{d^{3}x}{(2\pi {)}^3}{I}_{k_1}(\mathbf{x})I_{k_2}(\mathbf{x})I_{k_3}(\mathbf{x}).$$

$$I_{k_i}(\mathbf{x})\equiv {\int }_{k_i}{d}^{3}q\delta (\mathbf{q})e^{-i\mathbf{q}\cdot \mathbf{x}}\equiv \int d^{3}q{\tilde{I}}_{k_i}(\mathbf{q})e^{-i\mathbf{q}\cdot \mathbf{x}},$$

$${\tilde{P}}_n(k_1\cdots ,k_n)$$

$$=$$

$$V_{(12\cdots n)}^{\tilde{n}}={\int }_{k_1}{d}^3{q}_1\cdots {\int }_{k_n}{d}^3{q}_n{\delta }_D({\mathbf{q}}_{1\cdots n}).$$

$$V_{\tilde{n}}=\int \frac{d^{3}x}{(2\pi {)}^3}\prod _{j=1}^n{\iota }_{k_j}(\mathbf{x}),$$

$${\iota }_{k_j}(\mathbf{x})={\int }_{k_j}{d}^{3}q{e}^{-i\mathbf{q}\cdot \mathbf{x}}.$$

$$T(k_1,k_2,k_3,k_4,d_1,d_2)$$

$$=$$

$$\times {\delta }_D({\mathbf{p}}_1-{\mathbf{q}}_1+{\mathbf{q}}_2){\delta }_D({\mathbf{p}}_2-{\mathbf{q}}_1+{\mathbf{q}}_3){\delta }_{{\mathbf{q}}_1}{\delta }_{{\mathbf{q}}_2}{\delta }_{{\mathbf{q}}_3}{\delta }_{{\mathbf{q}}_4}$$

$$=$$

$$\times \int \frac{d^{3}z}{(2\pi {)}^3}{\iota }_{d_1}(\mathbf{z})I_{k_1}(\mathbf{x}-\mathbf{z}-\mathbf{y})I_{k_2}(\mathbf{x}+\mathbf{z}).$$

$$N_{n-\mathrm{poly}}=$$

$$=$$

$${\iota }_k={\int }_k{d}^{3}q{e}^{i\mathbf{q}\cdot \mathbf{x}}=4\pi k^2\left(\frac{\sin (kx)}{kx}\right)\delta k+\mathcal{O}(\delta k^3).$$

$$V_{\tilde{n}}=$$

$$=$$

$${\int }_0^{\infty }\frac{dx}{x^{n-2}}\prod _{i=1}^n\sin (k_{i}x)=\frac{\pi }{2}\left(\prod _{i=3}^n{\int }_{-k_i}^{k_i}\frac{d{\kappa }_i}{2}\right)$$

$$\times \left[{\delta }^D\left(\sum _{i=3}^n{\kappa }_i+k_1-k_2\right)-{\delta }^D\left(\sum _{i=3}^n{\kappa }_i-k_1-k_2\right)\right],$$

$$N_{\rm tri}=\frac{1}{k_{F}^{6}}\left\{\begin{array}[]{cc}4\pi^{2}k_{1}k_{2}k_{3}\left(\delta k\right)^{3}&,~{}k_{1}=k_{2}+k_{3}\\ 8\pi^{2}k_{1}k_{2}k_{3}\left(\delta k\right)^{3}&,~{}k_{1}<k_{2}+k_{3}\end{array}\right.\,,$$

$$N_{\mathrm{quad}}=$$

$$+3k_4-|k_1-k_2-k_3+k_4|).$$

$$N_{\max }=\frac{k_{\max }}{k_F}=\frac{1}{2\pi }{k}_{\max }{V}^{1/3}.$$

$$t_{Bispec}\approx 40{\left(\frac{k_{\max }}{0.2h/\text{Mpc}}\right)}^6{\left(\frac{V}{2.7{\text{ Gpc}}^3/h^3}\right)}^2\text{sec},$$

$$t_{Pn}\sim \frac{n^3}{27}{\left(\frac{k_{\max }}{0.2h/\text{Mpc}}\right)}^{n+3}{\left(\frac{V}{2.7{\text{ Gpc}}^3/h^3}\right)}^{(n+3)/3}\text{min},$$

$${\tilde{P}}_n(k_1\cdots ,k_n)=\frac{V_f}{V_{\tilde{n}}(2\pi {)}^3}{\int }_{k_1}{d}^3{q}_1\cdots {\int }_{k_n}{d}^3{q}_n{\delta }_D({\mathbf{q}}_{1\cdots n})\prod \delta ({\mathbf{q}}_i).$$

$${\delta }_D(x)=\underset{-\infty }{\overset{\infty }{\int }}\frac{d^{3}k}{(2\pi {)}^3}{e}^{ik\cdot x},$$

Full text

cosmology: large-scale structure of universe — methods: data analysis

††journal: ApJ

1 Introduction

The primordial density field in the Universe measured by the temperature anisotropies and polarizations of cosmic microwave background (CMB) radiation is very close to Gaussian (Planck Collaboration, 2016). In contrast, the observed galaxy density field shows strong non-Gaussianity due to late-time effects such as nonlinear gravitational instability, nonlinear galaxy bias as well as nonlinear redshift-space distortion.

The analysis of the galaxy density field in literature is yet focusing on the galaxy two-point correlation function, or its Fourier counterpart, the galaxy power spectrum (Tegmark et al., 2006; Blake et al., 2011; Parkinson et al., 2012; Reid et al., 2012; de la Torre et al., 2013; Pezzotta et al., 2017)). At the same time, the interest in the higher-order correlation functions and polyspectra has emerged to complement the two-point statistics in exploiting such a non-Gaussian galaxy density field. In Fourier space, these are, from third order, bispectrum, trispectrum, quadspectrum, pentaspectrum, hexaspectrum, and so on.

Thus far, the analysis on the higher-order correlation functions and polyspectra focuses mainly on the galaxy bispectrum and three-point correlation function (Scoccimarro et al., 2001; Verde et al., 2002; Gaztañaga et al., 2005; Nichol et al., 2006; Nishimichi et al., 2007; Marín, 2011; McBride et al., 2011b, a; Gil-Marín et al., 2015; Marín et al., 2013), but they have not been as informative as their lower-order counterparts. This is due to both the high computational cost in estimating higher-order correlation functions and the lack of accurate modeling needed for detailed data analysis.

One challenge in studying the higher-order polyspectrum in Fourier space is to find an efficient estimator; when estimating $n$-th order polyspectrum, the complexity of the naive estimator grows as the number of all possible $n$-gons in the three-dimensional Fourier space that scales as $N_{\rm max}^{3(n-1)}$. Here, $N_{\rm max}\equiv k_{\rm max}/\delta k$ is the number of one-dimensional discrete grid points of interest, with $k_{\rm max}$ being the maximum wavenumber of the analysis, usually set by the smallest scale that we can reliably model the nonlinearities in galaxy density field, and $\delta k$ being the size of the binning in Fourier space.

One can alleviate this computational cost by using the Fourier space representation of the Dirac-delta operator that dictates the statistical homogeneity. Originally done by Scoccimarro (2000) for estimating the bispectrum, in this way, we can estimate a polyspectrum by taking advantage of the Fast Fourier transformation algorithm. As we shall outline in Sec. 2, the calculation of the polyspectrum then boils down to a sum of products of the $N_{\rm max}$ distinct three-dimensional Fourier volumes derived from the original density field. This method avoids directly sampling all $n$-gons and is computationally much faster for all high-order polyspectra.

The implementation of the Fourier based method, however, requires allocating a large amount of memory. To avoid the aliasing problem associated with the Fourier method, the three-dimensional Fourier grid must be at least $(sN_{\rm max})^{3}$ for the $s$-point polyspectrum; as we need a total of $N_{\rm max}$ such three-dimensional Fourier volumes, the memory requirement grows as $s^{3}N_{\rm max}^{4}$. For example, the memory requirement to estimate the bispectrum ($s=3$) with $N_{\rm max}=128$, when using single precision, is already 27 Gigabytes.

In this paper, we shall improve the performance of the polyspectra estimators by parallelizing them. The most straightforward and naive approach might be to distribute each of the Fourier volumes to a different processor. To compute the polyspectum, then we need to calculate the sum of the multiplication of the Fourier volumes in different processors. This proves to be quite slow and, if not slower, only marginally faster than even the direct calculation from sampling all possible $n$-polygons. Instead, we distribute the Fourier volumes to parallel nodes so that the summation can be done almost entirely locally. Here, most of the inter-processor communication occurs during the Fourier transform step and the performance is significantly faster than the serial estimator. In this paper, we shall demonstrate the performance of this parallel algorithm for efficient calculation of higher-order polyspectra.

This paper is organized as follows. In Sec. 2 we discuss the Scoccimarro estimator and its extension to higher-order polyspectra and the limitations of such an extension. In Sec. 3 we explain our efficient method of parallelization for this estimator. Sec. 4 analyzes the results of our application of the estimator to a grid of unity, equivalent to calculating the number of polygons, for a variety of high-order polyspectra and also includes an analysis of the scaling of our parallel algorithm. We conclude in Sec. 5. App. A discusses the conversion of our general estimator from integral to summation form for direct application to calculations. Finally, App. B discusses some approximately analytic solutions to the number of $n$-gons in three-dimensional Fourier space.

We use the following convention

[TABLE]

for the Fourier transformation, and

[TABLE]

for the $n$-point polyspectrum. Here, $\delta^{D}$ is the Dirac-delta operator, and we use the shorthanded notation of ${\bm{k}}_{1\cdots n}\equiv{\bm{k}}_{1}+\cdots{\bm{k}}_{n}$. We denote the amplitude of a vector ${\bm{k}}_{i}$ as $k_{i}$.

2 The Scoccimarro estimator

The estimator for polyspectra that we shall parallelize here is based on the Scoccimarro estimator (Scoccimarro, 2000, 2015) that we summarize in this section. In this section, we shall mainly discuss the algorithm with the equations in the continuous limit and summarize the discrete implementation in App. A.

Let us begin with the bispectrum estimator. The estimator for the monopole bispectrum $B(k_{1},k_{2},k_{3})$ takes the average of $\delta(\bm{q}_{1})\delta(\bm{q}_{2})\delta(\bm{q}_{3})$ for all possible combinations of three vectors $|q_{i}-k_{i}|<\delta k_{i}/2$ ($i=1,2,3$) that satisfy the triangle condition ($\bm{q}_{123}=0$). We denote the bin size in $i$-th direction as $\delta k_{i}$, although, throughout the paper, we choose the bin size for all three directions the same for simplicity. The extension of the method to arbitrary bin size is trivial. We may formulate the bispectrum estimator as

[TABLE]

with $N_{(123)}$ being the number of triplets $(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3})$ satisfying the binning condition ($|q_{i}-k_{i}|<\delta k_{i}$ which is denoted as $q_{i}\sim k_{i}$) and the triangle condition: $N_{(123)}=\sum_{\bm{q}_{1},\bm{q}_{2},\bm{q}_{3}}\delta^{D}(\bm{q}_{123})$ with the same summation as Eq. (3). Throughout the paper, we impose the condition that $k_{1}\geq k_{2}\geq k_{3}$ to avoid the duplication.

In the continuous limit, the bispectrum estimator becomes (see, App. A for the details about the normalization factor)

[TABLE]

where $V_{f}$ is the volume of a fundamental Fourier cell, $(2\pi)^{3}/V$, and $V^{B}_{(123)}$ is the Fourier space volume defined as

[TABLE]

The key step for the Scoccimarro estimator is promoting the Dirac-delta operator with its Fourier representation:

[TABLE]

with which we can recast the three integrals for sampling triangles in Eq. (4) into single integration:

[TABLE]

Here,

[TABLE]

is the Fourier transformation of the function $\tilde{I}_{k_{i}}(\bm{q})$ that is defined with the density field in the spherical shell defined by $|q-k_{i}|<\delta k_{i}/2$ and zero otherwise. That is, instead of time-consuming operation of sampling all possible triangles, the Scoccimarro estimator involves Fourier transformation for computing $I_{k_{i}}(\bm{x})$ at each wavenumber bin, and summing over the product of three $I_{k_{i}}$s.

Note that the complexity, or the scaling of the computation time as a function of system size, of the Scoccimarro estimator is the same as the direct summation’s scaling of $N_{\rm max}^{6}$ for the case of bispectrum. The Scoccimarro estimator, however, reduces the computation time in two ways. First, we have fewer total numerical operations. The direct sampling involves looping over the full ranges of $-N_{\rm max}\sim N_{\rm max}$ for each component of ${\bm{k}}_{1}$ and $\pm N_{\rm max}/2\sim\pm N_{\rm max}$ for each component of ${\bm{k}}_{2}$, because our convention of $k_{1}\geq k_{2}\geq k_{3}$ with triangle condition constrains $k_{1}/2\leq k_{2}\leq k_{1}$. On the other hand, the Scoccimarro estimator reduces it to the matrix inner product ($N_{\rm max}^{3}$ operation) only for the amplitude triplets $(k_{1},k_{2},k_{3})$. Second, for the Scoccimarro estimator, the three-dimensional array $\delta(\bm{q})$ is accessed in an ordered manner that is easier for the CPU to cache. In contrast, the direct sampling accesses $\delta(\bm{q}_{3})$ in an irregular manner, which is inevitable because the value for $\bm{q}_{3}$ is determined by $\bm{q}_{3}=-\bm{q}_{1}-\bm{q}_{2}$. Combining these two factors, the serial version of the Scoccimarro estimator is a factor of 10 faster than the direct sampling in our implementation.

The reduction of operation time is even more striking as we calculate the higher order correlation functions. When computing the ‘angle-averaged’ $n$-point polyspectra that are only a function of the Fourier wavenumbers (Sefusatti, 2005), the Scoccimarro estimator’s complexity increases only by one more power of $N_{\rm max}$, ${\cal O}(N_{\rm max}^{n+3})$ while the direct sampling method’s complexity grows by three powers of $N_{\rm max}$ at each order: ${\cal O}(N_{\rm max}^{3(n-1)})$. Needless to say that the base operation time for Scoccimarro estimator is much faster than the direct sampling case, as we have already seen in the case of the bispectrum where the two method have the same complexity.

Extending the bispectrum estimator in Eq. (7) for the angle averaged higher-order polyspectra is straightforward:

[TABLE]

with the Fourier volume

[TABLE]

occupied by $\bm{q}_{1},\cdots,\bm{q}_{n}$ that contribute to the estimation of angle-averaged $n$-th order polyspectrum. Note that we can compute $V_{\tilde{n}}$ with the exact same method:

[TABLE]

with

[TABLE]

In this paper, we shall implement Eq. (11) to study the performance of the parallel algorithm.

2.1 On the dimensionality of $n$-th order polyspectra

At this point, a cautionary remark on the dimensionality of $n$-th order polygon is in order. While we only consider the ‘angle averaged’ polyspectra in this paper, specifying the full shape of $n$-point ($n\geq 3$) correlation function (or polyspectra, its Fourier counterpart) in the three-dimensional space would require $3n-6$ real numbers. It is because adding one more point to a configuration of $(n-1)$ points introduces three additional real numbers (the coordinate of the new point relative to the $n-1$ points) to specify the configuration of $n$ points. Starting the recursion from bispectrum (or three-point correlation function) that requires three numbers, therefore, each configuration of $n$-point correlation function is specified by $3n-6$ real numbers.

This argument does not apply for the two-point function ($n=2$) and the three point function ($n=3$) because of the underlying symmetry. For the two-point correlation function, translational symmetry removes three (we can always set the coordinate of the first point as origin) and the rotational symmetry removes two (correlation function must be independent of the orientation of the second point) real dimensions. As a result, we are left with one real dimension for the power spectrum. For the three-point correlation function, the rotational symmetry around the axis defined by existing two points removes one real dimension. As a result, we are left with three real dimensions for the bispectrum.

Along the same reasoning, adding redshift-space distortion (Kaiser, 1987) would require only two more real dimensions to the problem, because the coordinate of the newly added point is enough to determine the orientation of the point relative to the existing $n-1$ points. That is, from the two angles specifying the orientation of any three points relative to the line-of-sight direction, we can determine the relative orientation of all other points with respect to the line-of-sight direction.

Although we shall not discuss further in this paper, one can also implement the polyspectra estimator with the full dependence by using a similar method. For example, in addition to the amplitude of four wavevectors, the general trispectrum also depends on the diagonals $d_{1}=|\mathbf{k}_{1}-\mathbf{k}_{2}|$ and $d_{2}=|\mathbf{k}_{1}-\mathbf{k}_{3}|$, takes the form

[TABLE]

One may use, again, Fourier transformation to calculate the double convolution. Note that the result must reduce to the angle averaged trispectrum when integrating over all possible $d_{1}$ and $d_{2}$.

3 Parallelization

The key to the efficient parallelization of the polyspectra estimator is minimizing the inter-CPU communication, and we achieve that by using the slab decomposition scheme.

The slab decomposition slices the data cube along a single dimension, generally the fastest varying index, and assigns each chunk to different process. As each process stores a full range of other dimensions and limited range of the sliced dimension, this separates the data into ‘slabs’ of the domain hence the name slab decomposition. We show a visualization of this decomposition in Fig. 1. For the specific case of the polyspectra estimator, this decomposition takes the form of decomposing along the fastest varying Fourier grid dimension while the other two grid dimensions and the wavenumber dimension remain fully accessible in each process. For example, for the Julia implementation, we assign each process $\tilde{I}_{k_{i}}(k_{1\rm min}:k_{1\rm max},k_{\rm F}:k_{\rm max},k_{\rm F}:k_{\rm max})$, and for C implementation, we assign each process $\tilde{I}_{k_{i}}(k_{\rm F}:k_{\rm max},k_{\rm F}:k_{\rm max},k_{3\rm min}:k_{3\rm max})$ for all $i$ from 1 to $N_{\rm max}$.

The choice of slab decomposition makes reading data from a file very memory efficient as the main process can send data to the other processes partway through reading the data by simply sending each slab to the proper process as it is read in, or by having each process read in their own slab of data. This also means that at no point does the reading process need more than a single slab of data at once in memory which helps to keep the maximum memory usage low. This memory distribution is also natural for the Fourier transformation as FFTW (Frigo & Johnson, 2005) implements distributed Fourier transforms through slab decomposition and transposing, allowing for a very natural calculation without need to shuffle around memory between processes to get the right format. The only other step in the calculation is the product and sum, the product is always able to be done locally within each process and the sum is a trivially parallizable calculation with very little extra memory usage or communication as each processes calculates a local sum before the main processes collects and finalizes the sum. All together this memory distribution scheme allows us to easily and quickly calculate this estimator with minimal memory footprint and very little inter-CPU communication.

This method must be contrasted with the naïve parallelization that each process contains the full $\tilde{I}_{k_{i}}({\bm{k}})$ and $I_{k_{i}}(\bm{x})$ for each wavenumber bin $k_{i}$. For this alternative case, however, the time saved in the Fourier transform step is dwarfed by the significantly more expensive final step where we calculate sums of products across wavenumbers. As a result, this somewhat more natural decomposition is significantly slower than the slab decomposition.

4 Code Tests

In order to test the accuracy and performance of our implementation of the parallel estimator, we compute the number of polygons satisfying the homogeneity condition and the binning condition:

[TABLE]

for which we have present an analytical estimation (see Sec. 4.1 for the detail). Here, $V_{f}$ is the volume of the fundamental cell in the Fourier space, which is related to the fundamental Fourier wavenumber $k_{F}$ by $V_{f}=k_{F}^{3}$.

In practice, we estimate the number of polygons by setting $\delta({\bm{k}})=1$ for all Fourier modes and estimating the angle-averaged polyspectra. Therefore, counting the number of polygons requires exactly the same procedure as estimating the polyspectra.

4.1 Number of Polygons

To obtain the number of $n$-gons satisfying the binning conditions, we need to evaluate the $3n$-dimensional integration in Eq. (14). Again, using the Fourier representation of the Dirac-delta operator, we reduce the problem to the three-dimensional integration of multiplications of $\iota_{k}(\bm{x})$ (Eq. (11)), which is

[TABLE]

Using the expression to leading order in the bin size $\delta k$, we further reduce the calculation of the $3n$-dimensional volume to the one-dimensional integration:

[TABLE]

In App. B, We find a general expression for the radial integration,

[TABLE]

and apply Eq. (17) for $n=3$ and $n=4$ cases to find the number of triangles and quadrilaterals, which are given by

[TABLE]

and

[TABLE]

4.2 Result

In this section, we shall present number of polygons that we calculated using the parallelized polyspectra estimator and compare the result with the analytical estimation that we find in Sec. 4.1. We shall show the result for triangles ($n=3$), quadrilaterals ($n=4$), pentagons ($n=5$), and hexagons ($n=6$) that are estimated using, respectively, the parallel estimators of bispectrum, trispectrum, quadspectrum, and pentaspectrum.

To show the full shape dependence of the number of polygons as a function of their side length (wavenumbers), we need to visualize them in $n$-dimensional space of $(k_{1},k_{2},\cdots,k_{n})$. For our purpose, however, to compare the outcome of the parallel computation and analytical estimation, it suffices to flatten $n$-dimensional data points to the one-dimensional ordered sequence. Following the convention in McCullagh et al. (2016), we impose the condition $k_{i}\geq k_{i+1}$ to avoid duplication and form the sequence by taking the row-major order of the $n$-dimensional points ($k_{n}$ is the fastest varying dimension).

In all comparison plots that we present here, the $x$ axes show the index of the $n$-dimensional point $(k_{1},\cdots,k_{n})$ in the flattened sequence, and $y$ axes show the number of polygons whose $n$ sides are $(k_{1},\cdots,k_{n})$. We use blue color to show the numerical result and orange color to show the analytical estimation.

First, we show the number of triangles using the parallel bispectrum estimator in Fig. 2. The left panel shows the numerical result (blue) and analytical prediction (Eq. (18)), while the middle and right panels show the absolute number (middle) and fractional (right) differences. On large scales (or small indices), we expect the analytical estimation is not accurate because the bin’s size is of order the wavenumber ($k_{i}\simeq\delta k$). For small scales, although the absolute residuals increases with wavenumbers (the lowest order error term goes as $\prod k_{i}$ just like the $N_{\rm tri}$ in Eq. (18)), the fractional residual between the numerical and analytical estimation are consistently off by $10\sim 20\%$, and the biggest differences happen for the folded triangles ($k_{1}=k_{2}+k_{3}$). This result suggests that the variance of bispectrum estimated by using Eq. (18) (for example, in the Fisher information forecast in Desjacques et al. (2018)) must be off by the same factor. For an accurate estimation of the variance of the bispectrum, one must use the numerically estimated number of triangles.

We perform the same analysis for the number of quadraliterals by using the parallel estimator for trispectrum. The left panel of Fig. 3 shows the comparison between the numerical result and the analytical estimation in Eq. (19). It turns out that biggest the discrepancies are for the co-linear quadrilaterals ($k_{1}=k_{2}+k_{3}+k_{4}$) where the analytical formula, Eq. (19), predicts zero. To remedy the situation, we add the exact corrections for the co-linear quadrilaterals, Eq. (B29), and the analytical estimate indeed matches the numerical calculation better, as shown in the middle panel of Fig. 3.

We can easily extend the parallel estimator for the bispectrum and trispectrum to measure any $n>4$ polyspectra. We extend the estimator to $n=5$ and $n=6$ order polyspectra, and show the result in Fig. 4 and Fig. 5, respectively. In all cases the fractional errors between numerical calculation and analytical estimate are dominated by a collection of 100% error terms. These occur, just as the co-linear case for the quadrilaterals, because the leading-order approximation in the analytical calculation yields zeros when the actual calculation yields small but not-zero numbers. This happens for all $n>4$ cases, and it further reinforces the need for manually calculating the number of polygons when estimating the polyspectra, or estimating their variance.

4.3 Code performance

In this section, we present the performance of the parallel algorithm. All of our benchmarking work was done with 2.2 GHz Intel Xeon Processors on nodes with 40 CPU/node, 1 TB of RAM, 10 Gbps Ethernet and FDR Infiniband.

We have implemented two versions of parallel bispectrum estimator, one using C and the other using Julia. These two pieces of software perform the same algorithm. In Fig. 6, we compare the performance of the two implementations. Overall, as a function of domain length, the size of one-dimensional Fourier grid $N$ and the maximum wavenumber $N_{\rm max}=k_{\rm max}/k_{F}$, the performance scales as $\mathcal{O}(N^{3}N_{\rm max}^{n})$, as we have analyzed in Sec. 2. If we apply the aliasing condition which requires that $N>nN_{\rm max}$ we can rewrite the scaling as $\mathcal{O}(n^{3}N_{\rm max}^{n+3})$.

The C implementation (blue) out performs the Julia implementation (orange) for all domain lengths. The improvement, however, stays constant at about 30 seconds offset for reasonable grid sizes. As it is easier to read and modify the Julia version, we have extended the software to higher order polyspectra estimator only using Julia.

Fig. 7 shows the performance of our Julia implementation as a function of number of processors. Here, we use the domain size of $N=256$ (left panel) and $N=512$ (right panel). Initially, adding more CPUs improves the performance and reduces the execution time, but the improvement stops at some point. In fact, when the grid size per CPU is sufficiently small then adding more processes actually slows down the performance; the increased overhead time cancels the small gains in extra discretization.

How long does it take to estimate the bispectrum? In order to avoid the aliasing effect, the grid size must satisfy $N>3N_{\rm max}$ for estimating the bispectrum. Therefore, probing the bispectrum up to $N_{\rm max}\simeq 340$ requires the one-dimensional Fourier grid of $N=1024$; we can estimate the bispectrum in 20 minutes using 40 cores.

For a more concrete example, the HETDEX survey (Hill et al., 2008) has $k_{F}\approx 0.0045$ $h$/Mpc. Measuring the galaxy bispectrum to $k_{\rm max}=0.2$ $h$/Mpc only requires $N_{max}=45$. With the parallel version of Scoccimarro estimator, we can measure the galaxy bispectrum for HETDEX in $\sim$40 seconds with 2 GB of memory distributed over 40 processors. Here, we did not run the in-place FFTW, and kept both $I_{k_{i}}(\bm{x})$ and $\tilde{I}_{k_{i}}(\bm{q})$. Similarly, WFIRST (Spergel et al., 2015) has $k_{F}\approx 0.0028$ $h$/Mpc corresponding to $N_{\rm max}=75$ for $k_{\rm max}=0.2$ $h$/Mpc. For our implementation over 40 processors this corresponds to $\sim$50 seconds and 11 GB of memory distributed over 40 processors to calculate the galaxy bispectrum for WFIRST.

With the scaling that we have discussed in Sec. 2, we can extend the previous example to find a general relationship. For a cosmological survey with volume $V$ and maximum wavenumber $k_{\rm max}$, we calculate the number of one-dimensional sampling point as

[TABLE]

Because the performance of the Scoccimarro estimator scales as $\mathcal{O}(N_{\rm max}^{6})$, the estimated computing time becomes

[TABLE]

over 40 processors. For the higher-order polyspectra, we can extend this relationship by

[TABLE]

again, over the 40 processors. This relationship matches with the execution time for the polyspectra that we have measured for generating the figures in this paper.

5 Conclusion

In this paper, we present a novel parallel algorithm for estimating polyspectra of a complex field given on three-dimensional regular grids. While inheriting the efficiency of the Scoccimarro estimator, the parallel algorithm alleviates its rather stringent memory requirement by distributing the array into multiple processors. In addition, the key for efficient parallization is the slab-decomposition scheme that maintains only low level of inter-CPU data communication, avoiding massive data transportation.

Although we have only presented the case for the monopole (angle-averaged) polyspectra as a function of their side lengths (that is, the wavenumbers $k_{i}$), for estimating the redshift-space galaxy bispectrum, one can easily extend the method to the angular multipoles by applying the slab decomposition to the multipole-weighed fields $\delta_{n}(\bm{q})$ as described in Scoccimarro (2015). Beyond the galaxy bispectrum, measuring the polyspectra with the full (3$n$ - 6)-parameter dependencies requires little more elaboration, as the calculation involves convolution instead of matrix inner product (Eq. (13)). We shall leave this as a future project.

With the Julia implementation, we have demonstrated that for the galaxy surveys such as HETDEX and WFIRST, we can estimate the galaxy bispectrum down to $k\simeq 0.2\,{h/{\rm Mpc}}$ in about a minute, and even higher-order polyspectra in less than an hour time scales. With these execution time scales, it is feasible to study the high-order polyspectra by directly measuring them from the massive mock galaxy catalogs (Monaco et al., 2013; Kitaura et al., 2014; Avila et al., 2015; Izard et al., 2016; Munari et al., 2017; Agrawal et al., 2017; Taruya et al., 2018). If it is not for their own sake, studying the galaxy trispectrum (the four-point correlation function in Fourier space) and the galaxy pentaspectrum (the six-point correlation function) shall certainly elucidate the study of covariance matrix for, respectively, the galaxy power spectrum and the galaxy bispectrum.

We have stored our Julia implementation for the $n$-point polyspectra at (https://github.com/JosephTomlinson/PolyspectrumEstimator), and those who are interested in C implementation for the bispectrum estimator may contact the authors.

We thank Emiliano Sefusatti for useful discussion. DJ and JT were supported at Pennsylvania State University by NSF grant (AST-1517363) and NASA ATP program (80NSSC18K1103).

Appendix A Polyspectrum Estimator in Discrete Fourier Transformation

In the main manuscript, we have presented the polyspectra estimator in an contiguous limit, as an integral form. In the practical implementation, we have converted the estimator into a discrete summation as we present here. For our discrete implementation, following the FFTW convention, we use the unnormalized Discrete Fourier Transform and explicitly pull out all normalization factors. For more details, authors refer to Chapter 7 of Jeong (2010). Here, we note the discrete version of variables by square bracket $[\bm{x}]$.

Let us first summarize the continuous version of the polyspectra estimator. Starting with the integral form of the general estimator Eq. (9)

[TABLE]

We convert the Dirac-delta operator into its Fourier representation form:

[TABLE]

which, when input back into the estimator, results in

[TABLE]

where

[TABLE]

is the Fourier Transform of $I_{k_{i}}(\bm{q})$ which is defined to be the same as $\delta(\bm{q})$ within a spherical shell with Fourier-space radius $|\bm{q}|\simeq k_{i}$, and zero otherwise.

In the FFTW convention of un-normalized DFT, we find

[TABLE]

with which we get a partially discretized estimator

[TABLE]

We then fully discretize the estimator by converting the final integral to a sum:

[TABLE]

Finally, the quantity in parenthesis $\frac{V_{f}^{n-1}}{V_{\tilde{n}}}$ is the reciprocal of the number of polygons, and we can compute them by applying the same polyspectra estimator to a Fourier grid of unity. Note, however, that the normalization is different here.

[TABLE]

where we define $\tilde{\iota}_{k_{i}}[\bm{q}]$, analogous to $\tilde{I}_{k_{i}}[\bm{q}]$ above, as the Fourier grid which takes unity within the spherical shell bounded by $|k-k_{i}|<\delta k/2$ and zero elsewhere.

Appendix B Approximate Analytic Number of Polygons

In the polyspectra estimator, we have the normalization constant $\frac{V_{\tilde{n}}}{V_{f}^{n-1}}$, which is the number of $n$-gons that satisfy the closing condition (Dirac-delta) for a given set of side-length wavenumber $k_{i}$. In this appendix, we shall present the analytical approximation to this quantity.

Using the Dirac-delta (Eq. (A2)), we rewrite the $3n$-dimensional volume $V_{\tilde{n}}$ as

[TABLE]

for which we can compute the each $q_{i}$ integration as

[TABLE]

Combining the two results, we reduce the number of $n$-gons to one-dimensional integration as

[TABLE]

which is an analytically tractable integral, given the assumptions that $k_{i}\geq k_{i+1}$ (to uniquely define the polygon), and the polygon condition that $(\sum_{i\neq j}k_{i})\geq k_{j}$.

In what follows, we shall show the calculations. First, we show the general strategy of integrating ${\cal I}_{n}$ in App. B.1, which is followed by the explicit result for $n=3$ (App. B.2) and $n=4$ (App. B.3) cases. We then present the calculation for the colinear case, $k_{1}=k_{2}+\cdot+k_{n}$, where the approximation in Eq. (B2) resulting zero for all $n>4$ polygons. We therefore show the explicit expression without approximation in Eq. (B2) in App. B.4.

B.1 The analytical calculation: general case

To simplify the notation, let us define following function

[TABLE]

where we use that the integrand is an even function of $x$. Here, we show the general strategy of solving for the function ${\cal I}_{n}$. First, we use the following identity:

[TABLE]

to transform the integration as

[TABLE]

We then use

[TABLE]

to convert the integrations in terms of the Dirac-delta operator:

[TABLE]

That is,

[TABLE]

Going beyond this result requires the conditions to form polygons. For general polygons in 3D there is only one condition, the longest side must be less than or equal to the sum of all the other sides. If we don’t want to be bothered by the ordering, we can simply state as following: For any $j\in[1,n]$,

[TABLE]

With this condition, we can see that there must be a solution for $\mu_{i}\in(0,1)$ satisfying

[TABLE]

For later use, we find following identity for the integration of the Dirac delta useful:

[TABLE]

Here $\Theta(x_{1},x_{2},\cdot,x_{n})$ is a multi-dimensional Heaviside-Theta function, which is 1 only if none of the $x_{i}$ are not positive.

B.2 The analytical calculation for triangles ($n=3$)

For the bispectrum ($n=3$), Eq. (B9) reduces to

[TABLE]

Without loss of generality, let us order the wavenumbers so that $k_{1}\geq k_{2}\geq k_{3}$. Then, from the polygon inequality in Eq. (B10), the first term becomes

[TABLE]

and the second term is

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

B.3 The analytical calculation for trispectrum

For the trispectrum ($n=4$), Eq. (B9) reduces to

[TABLE]

Again, without loss of generality, we can order the wavenumbers so that $k_{1}\geq k_{2}\geq k_{3}\geq k_{4}$.

We first consider the case where $k_{1}=k_{2}+k_{3}+k_{4}$. The integration then becomes

[TABLE]

That is, the number of trispectrum using Eq. (B2) vanishes for the colinear quadraliterals. We present the expression without using the approximation in App. B.4.

For the other cases, $k_{1}<k_{2}+k_{3}+k_{4}$, we integrate ${\cal I}_{4}$ as following. The first Heaviside-Theta function vanishes unless,

[TABLE]

so that we integrate the first term as

[TABLE]

From Eq. (B10), we know that $k_{1}>-k_{2}-k_{3}+k_{4}$, and $-k_{1}<k_{2}-k_{3}+k_{4}$, and, of course, $k_{2}-k_{3}+k_{4}>-k_{2}-k_{3}+k_{4}$. Therefore,

[TABLE]

The second Heaviside-Theta function vanishes unless,

[TABLE]

Using that, we find, from Eq. (B10), $k_{1}<k_{2}+k_{2}+k_{4}$ and $-k_{1}<-k_{2}+k_{3}+k_{4}$, and of course, $-k_{2}+k_{3}+k_{4}<k_{2}+k_{3}+k_{4}$. Therefore,

[TABLE]

Combining all, we find that

[TABLE]

We then calculate the $V_{T}$ as

[TABLE]

B.4 Exact analytic colinear case

For the higher-order polyspectra, the lowest order theory predicts zero polygons for the colinear case $k_{1}=\sum_{2}^{n}k_{i}$ but this is clearly not the case through either direct summation tests or estimation through our method. Simply adding a few higher-order terms in the series of Eq. (B2) causes the integration to diverge so we instead need to solve the exact case directly.

The full calculation for the integration that we approximate in Eq. (B2) yields

[TABLE]

Restricting to the colinear case, we calculate the Fourier volume integration for the bispectrum:

[TABLE]

and for the trispectrum

[TABLE]

by using the Mathematica package.

For the Bispectrum case Eq. (B28) differs only about 1% from the lowest order estimation. For trispectrum, where the lowest order approximation predicts zero polygons, Eq. (B29) dramatically improves the match between the numerical calculation and analytical estimation, as we show in Fig. 3.