Efficient Parallel Algorithm for Estimating Higher-order Polyspectra

The primordial density field in the universe measured by the temperature anisotropies and polarizations of cosmic microwave background 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 the 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, 2011a; Marín et al. 2013; Gil-Marín et al. 2015), 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 nth order polyspectrum, the complexity of the naive estimator grows as the number of all possible n-gons in the three-dimensional Fourier space scales as ${N}_{\max }^{3(n-1)}$. Here, N_max ≡ k_max/δk is the number of one-dimensional discrete grid points of interest, with k_max being the maximum wavenumber of the analysis, usually set by the smallest scale that we can reliably model the nonlinearities in the galaxy density field, and δk being the size of the binning in Fourier space.

One can alleviate this computational cost 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 Section 2, the calculation of the polyspectrum then boils down to a sum of products of the N_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_max)³ for the s-point polyspectrum; as we need a total of N_max for such three-dimensional Fourier volumes, the memory requirement grows as ${s}^{3}{N}_{\max }^{4}$. For example, the memory requirement to estimate the bispectrum (s = 3) with N_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 Section 2 we discuss the Scoccimarro estimator and its extension to higher-order polyspectra and the limitations of such an extension. In Section 3 we explain our efficient method of parallelization for this estimator. Section 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 Section 5. Appendix A discusses the conversion of our general estimator from integral to summation form for direct application to calculations. Finally, Appendix B discusses some approximately analytic solutions to the number of n-gons in three-dimensional Fourier space.

We use the following convention

$$\begin{eqnarray}\begin{array}{rcl}\delta ({\boldsymbol{k}}) & = & \int {d}^{3}x\delta ({\boldsymbol{x}}){e}^{-i{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\\ \delta ({\boldsymbol{x}}) & = & \int \displaystyle \frac{{d}^{3}k}{{\left(2\pi \right)}^{3}}\delta ({\boldsymbol{k}}){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\end{array}\end{eqnarray} \tag{ 1 }$$

for the Fourier transformation, and

$$\begin{eqnarray}&&\left\langle \delta ({{\boldsymbol{k}}}_{1})\cdots \delta ({{\boldsymbol{k}}}_{n})\right\rangle ={\left(2\pi \right)}^{3}{P}_{n}({{\boldsymbol{k}}}_{1},\cdots ,\,{{\boldsymbol{k}}}_{n}){\delta }^{{\rm{D}}}({{\boldsymbol{k}}}_{1\cdots n})\end{eqnarray} \tag{ 2 }$$

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

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 Appendix 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 ({{\boldsymbol{q}}}_{1})\delta ({{\boldsymbol{q}}}_{2})\delta ({{\boldsymbol{q}}}_{3})$ for all possible combinations of three vectors $| {q}_{i}-{k}_{i}| \lt \delta {k}_{i}/2$ (i = 1, 2, 3) that satisfy the triangle condition (${{\boldsymbol{q}}}_{123}=0$). We denote the bin size in the ith direction as δk_i, although, throughout the paper, we choose the same bin size for all three directions for simplicity. The extension of the method to arbitrary bin size is trivial. We may formulate the bispectrum estimator as

$$\begin{eqnarray}\begin{array}{rcl}B({k}_{1},{k}_{2},{k}_{3}) & = & \displaystyle \frac{1}{{{VN}}_{(123)}}\displaystyle \sum _{{q}_{1}\sim {k}_{1}}\displaystyle \sum _{{q}_{2}\sim {k}_{2}}\displaystyle \sum _{{q}_{3}\sim {k}_{3}}\\ & & \times \,\delta ({{\boldsymbol{q}}}_{1})\delta ({{\boldsymbol{q}}}_{2})\delta ({{\boldsymbol{q}}}_{3}){\delta }^{{\rm{D}}}({{\boldsymbol{q}}}_{123}),\end{array}\end{eqnarray} \tag{ 3 }$$

with N₍₁₂₃₎ being the number of triplets $({{\boldsymbol{q}}}_{1},{{\boldsymbol{q}}}_{2},{{\boldsymbol{q}}}_{3})$ satisfying the binning condition ($| {q}_{i}-{k}_{i}| \lt \delta {k}_{i}$, which is denoted as q_i ∼ k_i) and the triangle condition ${N}_{(123)}={\sum }_{{{\boldsymbol{q}}}_{1},{{\boldsymbol{q}}}_{2},{{\boldsymbol{q}}}_{3}}{\delta }^{{\rm{D}}}({{\boldsymbol{q}}}_{123})$ with the same summation as Equation (3). Throughout the paper, we impose the condition that ${k}_{1}\geqslant {k}_{2}\geqslant {k}_{3}$ to avoid the duplication.

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

$$\begin{eqnarray}\begin{array}{rcl}B({k}_{1},{k}_{2},{k}_{3}) & = & \displaystyle \frac{{V}_{f}}{{V}_{(123)}^{B}{\left(2\pi \right)}^{3}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}{\int }_{{k}_{2}}{d}^{3}{q}_{2}{\int }_{{k}_{3}}{d}^{3}{q}_{3}\\ & & \times \,\delta ({{\boldsymbol{q}}}_{1})\delta ({{\boldsymbol{q}}}_{2})\delta ({{\boldsymbol{q}}}_{3}){\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{123}),\end{array}\end{eqnarray} \tag{ 4 }$$

where V_f is the volume of a fundamental Fourier cell, (2π)³/V, and ${V}_{(123)}^{B}$ is the Fourier space volume defined as

$$\begin{eqnarray}&&{V}_{(123)}^{B}={\int }_{{k}_{1}}{d}^{3}{q}_{1}{\int }_{{k}_{2}}{d}^{3}{q}_{2}{\int }_{{k}_{3}}{d}^{3}{q}_{3}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{123}).\end{eqnarray} \tag{ 5 }$$

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

$$\begin{eqnarray}&&{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{123})=\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}{e}^{-i{\boldsymbol{x}}\cdot ({{\boldsymbol{q}}}_{1}+{{\boldsymbol{q}}}_{2}+{{\boldsymbol{q}}}_{3})},\end{eqnarray} \tag{ 6 }$$

with which we can recast the three integrals for sampling triangles in Equation (4) into a single integration,

$$\begin{eqnarray}&&B({k}_{1},{k}_{2},{k}_{3})=\displaystyle \frac{{V}_{f}}{{V}_{(123)}^{B}{\left(2\pi \right)}^{3}}\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}{I}_{{k}_{1}}({\boldsymbol{x}}){I}_{{k}_{2}}({\boldsymbol{x}}){I}_{{k}_{3}}({\boldsymbol{x}}).\end{eqnarray} \tag{ 7 }$$

Here,

$$\begin{eqnarray}&&{I}_{{k}_{i}}({\boldsymbol{x}})\equiv {\int }_{{k}_{i}}{d}^{3}q\delta ({\boldsymbol{q}}){e}^{-i{\boldsymbol{q}}\cdot {\boldsymbol{x}}}\equiv \int {d}^{3}q{\tilde{I}}_{{k}_{i}}({\boldsymbol{q}}){e}^{-i{\boldsymbol{q}}\cdot {\boldsymbol{x}}}\end{eqnarray} \tag{ 8 }$$

is the Fourier transformation of the function ${\tilde{I}}_{{k}_{i}}({\boldsymbol{q}})$ that is defined with the density field in the spherical shell defined by $| q-{k}_{i}| \lt \delta {k}_{i}/2$ and zero otherwise. That is, instead of the time-consuming operation of sampling all possible triangles, the Scoccimarro estimator involves Fourier transformation for computing ${I}_{{k}_{i}}({\boldsymbol{x}})$ at each wavenumber bin, and summing over the product of three ${I}_{{k}_{i}}{\rm{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}_{\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_max ∼ N_max for each component of ${{\boldsymbol{k}}}_{1}$ and ±N_max/2 ∼ ±N_max for each component of ${{\boldsymbol{k}}}_{2}$, because our convention of ${k}_{1}\geqslant {k}_{2}\geqslant {k}_{3}$ with a triangle condition constrains ${k}_{1}/2\leqslant {k}_{2}\leqslant {k}_{1}$. On the other hand, the Scoccimarro estimator reduces it to the matrix inner product (${N}_{\max }^{3}$ operation) only for the amplitude triplets $({k}_{1},{k}_{2},{k}_{3})$. Second, for the Scoccimarro estimator, the three-dimensional array $\delta ({\boldsymbol{q}})$ is accessed in an ordered manner that is easier for the CPU to cache. In contrast, the direct sampling accesses $\delta ({{\boldsymbol{q}}}_{3})$ in an irregular manner, which is inevitable because the value for ${{\boldsymbol{q}}}_{3}$ is determined by ${{\boldsymbol{q}}}_{3}=-{{\boldsymbol{q}}}_{1}-{{\boldsymbol{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_max, ${ \mathcal O }({N}_{\max }^{n+3})$ while the direct sampling method's complexity grows by three powers of N_max at each order, ${ \mathcal O }({N}_{\max }^{3(n-1)})$. Needless to say that the base operation time for the 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 Equation (7) for the angle-averaged higher-order polyspectra is straightforward:

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{P}}_{n}({k}_{1}\cdots ,\,{k}_{n})\\ =\ \displaystyle \frac{{V}_{f}}{{V}_{(12\cdots n)}^{\tilde{n}}{\left(2\pi \right)}^{3}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n})\prod \delta ({{\boldsymbol{q}}}_{i})\\ =\ \displaystyle \frac{{V}_{f}}{{V}_{(12\cdots n)}^{\tilde{n}}{\left(2\pi \right)}^{3}}\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}{I}_{{k}_{1}}({\boldsymbol{x}})\cdots {I}_{{k}_{n}}({\boldsymbol{x}})\end{array}\end{eqnarray} \tag{ 9 }$$

with the Fourier volume

$$\begin{eqnarray}&&{V}_{(12\cdots n)}^{\tilde{n}}={\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n})\end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}&&{V}_{\tilde{n}}=\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}\displaystyle \prod _{j=1}^{n}{\iota }_{{k}_{j}}({\boldsymbol{x}}),\end{eqnarray} \tag{ 11 }$$

with

$$\begin{eqnarray}&&{\iota }_{{k}_{j}}({\boldsymbol{x}})={\int }_{{k}_{j}}{d}^{3}{{qe}}^{-i{\boldsymbol{q}}\cdot {\boldsymbol{x}}}.\end{eqnarray} \tag{ 12 }$$

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

2.1. On the Dimensionality of nth Order Polyspectra

At this point, a cautionary remark on the dimensionality of nth order polygon is in order. While we only consider the angle-averaged polyspectra in this paper, specifying the full shape of the n-point (n ≥ 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 the three-point correlation function) that requires three numbers, therefore, each configuration of the 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 the origin) and the rotational symmetry removes two (the 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 using a similar method. For example, in addition to the amplitude of four wavevectors, the general trispectrum also depends on the diagonals ${d}_{1}=| {{\boldsymbol{k}}}_{1}-{{\boldsymbol{k}}}_{2}|$ and ${d}_{2}=| {{\boldsymbol{k}}}_{1}-{{\boldsymbol{k}}}_{3}|$, and takes the form

$$\begin{eqnarray}&&\begin{array}{l}T({k}_{1},{k}_{2},{k}_{3},{k}_{4},{d}_{1},{d}_{2})\\ \quad =\,\displaystyle \frac{{V}_{f}}{{V}_{T}{\left(2\pi \right)}^{3}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{4}}{d}^{3}{q}_{4}{\int }_{{d}_{1}}{d}^{3}{p}_{1}{\int }_{{d}_{2}}{d}^{3}{p}_{2}{\delta }_{D}({{\boldsymbol{q}}}_{1234})\\ \qquad \times \ {\delta }_{{\rm{D}}}({{\boldsymbol{p}}}_{1}-{{\boldsymbol{q}}}_{1}+{{\boldsymbol{q}}}_{2}){\delta }_{{\rm{D}}}({{\boldsymbol{p}}}_{2}-{{\boldsymbol{q}}}_{1}+{{\boldsymbol{q}}}_{3}){\delta }_{{{\boldsymbol{q}}}_{1}}{\delta }_{{{\boldsymbol{q}}}_{2}}{\delta }_{{{\boldsymbol{q}}}_{3}}{\delta }_{{{\boldsymbol{q}}}_{4}}\\ \quad =\ \displaystyle \frac{{V}_{f}}{{V}_{T}{\left(2\pi \right)}^{3}}\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}{I}_{{k}_{4}}({\boldsymbol{x}})\int \displaystyle \frac{{d}^{3}y}{{\left(2\pi \right)}^{3}}{\iota }_{{d}_{2}}({\boldsymbol{y}}){I}_{{k}_{3}}({\boldsymbol{x}}+{\boldsymbol{y}})\\ \qquad \times \ \int \displaystyle \frac{{d}^{3}z}{{\left(2\pi \right)}^{3}}{\iota }_{{d}_{1}}({\boldsymbol{z}}){I}_{{k}_{1}}({\boldsymbol{x}}-{\boldsymbol{z}}-{\boldsymbol{y}}){I}_{{k}_{2}}({\boldsymbol{x}}+{\boldsymbol{z}}).\end{array}\end{eqnarray} \tag{ 13 }$$

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₁ and d₂.

The key to the efficient parallelization of the polyspectra estimator is minimizing the inter-CPU communication, and we achieve that 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 a different process. As each process stores a full range of other dimensions and a 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 Figure 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\min }\,:{k}_{1\max },{k}_{{\rm{F}}}\,:{k}_{\max },{k}_{{\rm{F}}}\,:{k}_{\max })$, and for C implementation, we assign each process ${\tilde{I}}_{{k}_{i}}({k}_{{\rm{F}}}\,:{k}_{\max },{k}_{{\rm{F}}}\,:{k}_{\max },{k}_{3\min }\,:{k}_{3\max })$ for all i from 1 to N_max. As the arrays ${\tilde{I}}_{{k}_{i}}({\boldsymbol{k}})$ are zero except for the shell defined as $| {\boldsymbol{k}}| \sim {k}_{i}$, for which we assign the value from the original data cube $\delta ({\boldsymbol{k}})$, generating N_max of ${N}_{\max }^{3}$-cube arrays takes only the time that we read through the original data cube. That is, at no point does the generating process need more than a single slab of data at once in memory.

Zoom In Zoom Out Reset image size

This memory distribution is also natural for the Fourier transformation as the Fastest Fourier Transform in the West (FFTW; Frigo & Johnson 2005) implements distributed Fourier transforms through slab decomposition and transposing, allowing for a very natural calculation without a 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 parallelizable 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}}({\boldsymbol{k}})$ and ${I}_{{k}_{i}}({\boldsymbol{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.

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

$$\begin{eqnarray}\begin{array}{rcl}{N}_{n-\mathrm{poly}} & = & \displaystyle \frac{{V}_{\tilde{n}}}{{V}_{f}^{n-1}}\\ & = & \displaystyle \frac{1}{{k}_{{\rm{F}}}^{3n-3}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n}),\end{array}\end{eqnarray} \tag{ 14 }$$

for which we have present an analytical estimation (see Section 4.1 for details). 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}_{{\rm{F}}}^{3}$.

In practice, we estimate the number of polygons by setting $\delta ({\boldsymbol{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 Equation (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}({\boldsymbol{x}})$ (Equation (11)), which is

$$\begin{eqnarray}&&{\iota }_{k}={\int }_{k}{d}^{3}{{qe}}^{i{\boldsymbol{q}}\cdot {\boldsymbol{x}}}=4\pi {k}^{2}\left(\displaystyle \frac{\sin ({kx})}{{kx}}\right)\delta k+{ \mathcal O }(\delta {k}^{3}).\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\tilde{n}} & = & \int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}\displaystyle \prod _{i=1}^{n}{\iota }_{{k}_{i}}(x)\\ & = & {\left(4\pi \delta k\right)}^{n}{\int }_{0}^{\infty }\displaystyle \frac{{dx}}{2{\pi }^{2}{x}^{n-2}}\displaystyle \prod _{i=1}^{n}{k}_{i}\sin ({k}_{i}x).\end{array}\end{eqnarray} \tag{ 16 }$$

In Appendix B, we find a general expression for the radial integration,

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{0}^{\infty }\displaystyle \frac{{dx}}{{x}^{n-2}}\displaystyle \prod _{i=1}^{n}\sin ({k}_{i}x)=\displaystyle \frac{\pi }{2}\left(\displaystyle \prod _{i=3}^{n}{\int }_{-{k}_{i}}^{{k}_{i}}\displaystyle \frac{d{\kappa }_{i}}{2}\right)\\ \quad \times \,\left[{\delta }^{{\rm{D}}}\left(\displaystyle \sum _{i=3}^{n}{\kappa }_{i}+{k}_{1}-{k}_{2}\right)-{\delta }^{{\rm{D}}}\left(\displaystyle \sum _{i=3}^{n}{\kappa }_{i}-{k}_{1}-{k}_{2}\right)\right],\end{array}\end{eqnarray} \tag{ 17 }$$

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

$$\begin{eqnarray}{N}_{\mathrm{tri}}=\displaystyle \frac{1}{{k}_{{\rm{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}\lt {k}_{2}+{k}_{3},\end{array}\right.\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{quad}} & = & \displaystyle \frac{1}{{k}_{{\rm{F}}}^{9}}8{\pi }^{3}{\left(\delta k\right)}^{4}{k}_{1}{k}_{2}{k}_{3}{k}_{4}\left(-{k}_{1}+{k}_{2}+{k}_{3}\right.\\ & & \left.+3{k}_{4}-| {k}_{1}-{k}_{2}-{k}_{3}+{k}_{4}| \right).\end{array}\end{eqnarray} \tag{ 19 }$$

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 Section 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}\geqslant {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 to show the numerical result and orange to show the analytical estimation.

First, we show the number of triangles using the parallel bispectrum estimator in Figure 2. The left panel shows the numerical result (blue) and analytical prediction (Equation (18)), while the middle and right panels show the absolute number (middle) and fractional (right) differences. On large scales (or small indices), we expect that the analytical estimation is not accurate because the bin's size is of order the wavenumber (k_i ≃ δ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_tri in Equation (18)), the fractional residual between the numerical and analytical estimation are consistently off by 10% ∼ 20%, and the biggest differences happen for the folded triangles (k₁ = k₂ + k₃). This result suggests that the variance of bispectrum estimated using Equation (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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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 Figures 4 and 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 colinear case for the quadrilaterals, because the leading-order approximation in the analytical calculation yields zeros when the actual calculation yields small but non-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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 the parallel bispectrum estimator, one using C and the other using Julia. These two pieces of software perform the same algorithm. In Figure 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_max = k_max/k_F, the performance scales as ${ \mathcal O }({N}^{3}{N}_{\max }^{n})$, as we have analyzed in Section 2. If we apply the aliasing condition, which requires that $N\gt {{nN}}_{\max }$, we can rewrite the scaling as ${ \mathcal O }({n}^{3}{N}_{\max }^{n+3})$.

Zoom In Zoom Out Reset image size

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

Figure 7 shows the performance of our Julia implementation as a function of the 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.

Zoom In Zoom Out Reset image size

How long does it take to estimate the bispectrum? In order to avoid the aliasing effect, the grid size must satisfy $N\gt 3{N}_{\max }$ for estimating the bispectrum. Therefore, probing the bispectrum up to N_max ≃ 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 Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) survey (Hill et al. 2008) has k_F ≈ 0.0045 h/Mpc. Measuring the galaxy bispectrum to k_max = 0.2 h/Mpc only requires N_max = 45. With the parallel version of the Scoccimarro estimator, we can measure the galaxy bispectrum for HETDEX in ∼40 s with 2 GB of memory distributed over 40 processors. Here, we did not run the FFTW in place, and kept both ${I}_{{k}_{i}}({\boldsymbol{x}})$ and ${\tilde{I}}_{{k}_{i}}({\boldsymbol{q}})$. Similarly, the Wide Field Infrared Survey Telescope (WFIRST; Spergel et al. 2015) has k_F ≈ 0.0028 h/Mpc corresponding to ${N}_{\max }=75$ for k_max = 0.2 h/Mpc. For our implementation over 40 processors this corresponds to ∼50 s and 11 GB of memory distributed over 40 processors to calculate the galaxy bispectrum for WFIRST.

With the scaling that we have discussed in Section 2, we can extend the previous example to find a general relationship. For a cosmological survey with volume V and maximum wavenumber k_max, we calculate the number of one-dimensional sampling points as

$$\begin{eqnarray}&&{N}_{\max }=\displaystyle \frac{{k}_{\max }}{{k}_{{\rm{F}}}}=\displaystyle \frac{1}{2\pi }{k}_{\max }{V}^{1/3}.\end{eqnarray} \tag{ 20 }$$

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

$$\begin{eqnarray}&&{t}_{\mathrm{Bispec}}\approx 40{\left(\displaystyle \frac{{k}_{\max }}{0.2h/\mathrm{Mpc}}\right)}^{6}{\left(\displaystyle \frac{V}{2.7{\mathrm{Gpc}}^{3}/{h}^{3}}\right)}^{2}\,{\rm{s}}\end{eqnarray} \tag{ 21 }$$

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

$$\begin{eqnarray}&&{t}_{\mathrm{Pn}}\sim \displaystyle \frac{{n}^{3}}{27}{\left(\displaystyle \frac{{k}_{\max }}{0.2h/\mathrm{Mpc}}\right)}^{n+3}{\left(\displaystyle \frac{V}{2.7{\mathrm{Gpc}}^{3}/{h}^{3}}\right)}^{(n+3)/3}\,\mathrm{minutes},\end{eqnarray} \tag{ 22 }$$

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.

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 a 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}({\boldsymbol{q}})$ as described in Scoccimarro (2015). Beyond the galaxy bispectrum, measuring the polyspectra with the full (3n − 6) parameter dependencies requires a little more elaboration, as the calculation involves convolution instead of a matrix inner product (Equation (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 ≃ 0.2 h/Mpc in about a minute, and even higher-order polyspectra in timescales of less than an hour. With these execution timescales, 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; Agrawal et al. 2017; Munari 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 the C implementation for the bispectrum estimator may contact the authors.

We thank Emiliano Sefusatti for useful discussions. D.J. and J.T. were supported at Pennsylvania State University by NSF grant (AST-1517363) and NASA ATP program (80NSSC18K1103).

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 $[{\boldsymbol{x}}]$.

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

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{P}}_{n}({k}_{1}\cdots ,\,{k}_{n})=\displaystyle \frac{{V}_{f}}{{V}_{\tilde{n}}{\left(2\pi \right)}^{3}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n})\\ \qquad \times \ \prod \delta ({{\boldsymbol{q}}}_{i}).\end{array}\end{eqnarray} \tag{ 23 }$$

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

$$\begin{eqnarray}&&{\delta }_{{\rm{D}}}(x)={\int }_{-\infty }^{\infty }\displaystyle \frac{{d}^{3}k}{{\left(2\pi \right)}^{3}}{e}^{{ik}\cdot x},\end{eqnarray} \tag{ 24 }$$

which, when inputed back into the estimator, results in

$$\begin{eqnarray}&&{\tilde{P}}_{n}({k}_{1}\cdots ,\,{k}_{n})=\displaystyle \frac{{V}_{f}}{{V}_{\tilde{n}}{\left(2\pi \right)}^{3}}\int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}\prod {I}_{{k}_{i}}(x),\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{I}_{{k}_{i}}(x)={\int }_{{k}_{i}}{d}^{3}{q}_{i}{e}^{i{\boldsymbol{x}}\cdot {{\boldsymbol{q}}}_{i}}\delta ({{\boldsymbol{q}}}_{i})=\int {d}^{3}{q}_{i}{e}^{i{\boldsymbol{x}}\cdot {{\boldsymbol{q}}}_{i}}{\tilde{I}}_{{k}_{i}}({{\boldsymbol{q}}}_{i}),\end{eqnarray} \tag{ 26 }$$

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

In the FFTW convention of unnormalized discrete Fourier transform, we find

$$\begin{eqnarray}&&{I}_{{k}_{i}}[{\boldsymbol{x}}]=\displaystyle \frac{{\left(2\pi \right)}^{3}}{{N}^{3}}\mathrm{FFTW}\left({\tilde{I}}_{{k}_{i}}[{\boldsymbol{q}}]\right),\end{eqnarray} \tag{ 27 }$$

with which we get a partially discretized estimator of

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{P}}_{n}({k}_{1}\cdots ,\,{k}_{n})=\displaystyle \frac{{V}_{f}}{{V}_{\tilde{n}}{\left(2\pi \right)}^{3}}{\left(\displaystyle \frac{{\left(2\pi \right)}^{3}}{{N}^{3}}\right)}^{n}\\ \qquad \times \ \int \displaystyle \frac{{d}^{3}x}{{\left(2\pi \right)}^{3}}\displaystyle \prod _{i=1}^{n}\mathrm{FFTW}\left({\tilde{I}}_{{k}_{i}}[{\boldsymbol{q}}]\right).\end{array}\end{eqnarray} \tag{ 28 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{P}}_{n}({k}_{1}\cdots ,\,{k}_{n}) & = & \left(\displaystyle \frac{{V}_{f}^{n-1}}{{V}_{\tilde{n}}}\right)\displaystyle \frac{{\left(2\pi \right)}^{3n}}{{V}_{f}^{n-2}{N}^{3n}}\displaystyle \frac{V}{{\left(2\pi \right)}^{6}{N}^{3}}\\ & & \times \ \displaystyle \sum _{[{\boldsymbol{x}}]}\displaystyle \prod _{i=1}^{n}\mathrm{FFTW}\left({\tilde{I}}_{{k}_{i}}[{\boldsymbol{q}}]\right)\\ & = & \left(\displaystyle \frac{{V}_{f}^{n-1}}{{V}_{\tilde{n}}}\right)\displaystyle \frac{{V}^{n-1}}{{N}^{3n}}\displaystyle \frac{1}{{N}^{3}}\displaystyle \sum _{[{\boldsymbol{x}}]}\displaystyle \prod _{i=1}^{n}\mathrm{FFTW}\left({\tilde{I}}_{{k}_{i}}[{\boldsymbol{q}}]\right).\end{array}\end{eqnarray} \tag{ 29 }$$

Finally, the quantity in parenthesis $\tfrac{{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.

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{V}_{\tilde{n}}}{{V}_{f}^{n-1}} & = & \displaystyle \frac{1}{{V}_{f}^{n-1}}{\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n})\\ & = & \displaystyle \frac{1}{{N}^{3}{V}_{f}^{n-1}}\displaystyle \sum _{[{\boldsymbol{x}}]}\displaystyle \prod _{i=1}^{n}\mathrm{FFTW}\left({\tilde{\iota }}_{{k}_{i}}[{\boldsymbol{q}}]\right).\end{array}\end{eqnarray} \tag{ 30 }$$

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

In the polyspectra estimator, we have the normalization constant $\tfrac{{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 (Equation (24)), we rewrite the 3n-dimensional volume ${V}_{\tilde{n}}$ as

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\tilde{n}} & = & {\int }_{{k}_{1}}{d}^{3}{q}_{1}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{\delta }_{{\rm{D}}}({{\boldsymbol{q}}}_{1\cdots n})\\ & = & \displaystyle \frac{1}{{\left(2\pi \right)}^{3}}\int {{dx}}^{3}{\int }_{{k}_{1}}{d}^{3}{q}_{1}{e}^{i{\boldsymbol{x}}\cdot {{\boldsymbol{q}}}_{1}}\cdots {\int }_{{k}_{n}}{d}^{3}{q}_{n}{e}^{i{\boldsymbol{x}}\cdot {{\boldsymbol{q}}}_{n}},\end{array}\end{eqnarray} \tag{ 31 }$$

for which we can compute the each q_i integration as

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{k}{d}^{3}{{qe}}^{i{\boldsymbol{q}}\cdot {\boldsymbol{x}}} & = & 4\pi {\int }_{k-\delta k/2}^{k+\delta k/2}{{dqq}}^{2}\displaystyle \frac{\sin ({qx})}{{qx}}\\ & \approx & 4\pi \displaystyle \frac{k\sin ({kx})}{x}\delta k+{ \mathcal O }(\delta {k}^{3}).\end{array}\end{eqnarray} \tag{ 32 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}{V}_{\tilde{n}} & = & {2}^{2n-3}{\pi }^{n-3}\left(\displaystyle \prod _{i=1}^{n}{k}_{i}\delta k\right)\int {d}^{3}x\displaystyle \frac{1}{{x}^{n}}\displaystyle \prod _{i=1}^{n}\sin \left({k}_{i}x\right)\\ & = & {2}^{2n-1}{\pi }^{n-2}\left(\displaystyle \prod _{i=1}^{n}{k}_{i}\delta k\right){\int }_{0}^{\infty }\displaystyle \frac{{dx}}{{x}^{n-2}}\displaystyle \prod _{i=1}^{n}\sin \left({k}_{i}x\right),\end{array}\end{eqnarray} \tag{ 33 }$$

which is an analytically tractable integral, given the assumptions that ${k}_{i}\geqslant {k}_{i+1}$ (to uniquely define the polygon), and the polygon condition that $({\sum }_{i\ne j}{k}_{i})\geqslant {k}_{j}$.

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

B.1. The Analytical Calculation: General Case

To simplify the notation, let us define following function

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{n}\left({k}_{1},\cdots ,\,{k}_{n}\right) & \equiv & {\int }_{0}^{\infty }{dx}\displaystyle \frac{\prod \sin ({k}_{i}x)}{{x}^{n-2}}\\ & = & \displaystyle \frac{1}{2}{\int }_{-\infty }^{\infty }\displaystyle \frac{{dx}}{{x}^{n-2}}\displaystyle \prod _{i=1}^{n}\sin \left({k}_{i}x\right),\end{array}\end{eqnarray} \tag{ 34 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\sin ({kx})}{{kx}}={\int }_{-1}^{1}\displaystyle \frac{d\mu }{2}\,{e}^{{ikx}\mu }\end{eqnarray} \tag{ 35 }$$

to transform the integration as

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal I }}_{n}({k}_{1},\cdots ,\,{k}_{n})={k}_{1}\cdots {k}_{n-2}{\int }_{-\infty }^{\infty }\displaystyle \frac{{dx}}{2}{\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{1}}{2}\cdots \\ \times \,{\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{n-2}}{2}{e}^{{ix}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}\right)}\sin \left({k}_{n-1}x\right)\sin \left({k}_{n}x\right).\end{array}\end{eqnarray} \tag{ 36 }$$

We then use

$$\begin{eqnarray}\begin{array}{rcl}\sin \left({k}_{n-1}x\right)\sin \left({k}_{n}x\right) & = & \displaystyle \frac{1}{{\left(2i\right)}^{2}}\left[{e}^{{{ik}}_{n-1}x}-{e}^{-{{ik}}_{n-1}x}\right]\\ & & \times \,\left[{e}^{{{ik}}_{n}x}-{e}^{-{{ik}}_{n}x}\right]\\ & = & \displaystyle \frac{1}{4}\left({e}^{i({k}_{n-1}-{k}_{n})x}+{e}^{-i({k}_{n-1}-{k}_{n})x}\right.\\ & & -\ \left.{e}^{i({k}_{n-1}+{k}_{n})x}-{e}^{-i({k}_{n-1}+{k}_{n})x}\right),\end{array}\end{eqnarray} \tag{ 37 }$$

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

$$\begin{eqnarray}&&{\delta }^{{\rm{D}}}(k)={\int }_{-\infty }^{\infty }\displaystyle \frac{{dx}}{2\pi }{e}^{{ikx}}.\end{eqnarray} \tag{ 38 }$$

That is,

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal I }}_{n}({k}_{1},\cdots ,{k}_{n})=\displaystyle \frac{\pi }{4}{k}_{1}\cdots {k}_{n-2}{\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{1}}{2}\cdots {\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{n-2}}{2}\\ \quad \times \,\left[{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}+{k}_{n-1}-{k}_{n}\right)\right.\\ \quad +\,{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}-{k}_{n-1}+{k}_{n}\right)\\ \quad -\,{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}+{k}_{n-1}+{k}_{n}\right)\\ \quad -\ \left.{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}-{k}_{n-1}-{k}_{n}\right)\right]\\ \,=\,\displaystyle \frac{\pi }{2}{k}_{1}\cdots {k}_{n-2}{\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{1}}{2}\cdots {\int }_{-1}^{1}\displaystyle \frac{d{\mu }_{n-2}}{2}\\ \quad \times \,\left[{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}+{k}_{n-1}-{k}_{n}\right)\right.\\ \quad -\,\left.{\delta }^{{\rm{D}}}\left({k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}-{k}_{n-1}-{k}_{n}\right)\right]\\ \,=\,\displaystyle \frac{\pi }{2}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\cdots {\int }_{-{k}_{n-2}}^{{k}_{n-2}}\displaystyle \frac{d{\kappa }_{n-2}}{2}\\ \quad \times \,\left[{\delta }^{{\rm{D}}}\left({\kappa }_{1}+\cdots +{\kappa }_{n-2}+{k}_{n-1}-{k}_{n}\right)\right.\\ \quad -\,\left.{\delta }^{{\rm{D}}}\left({\kappa }_{1}+\cdots +{\kappa }_{n-2}-{k}_{n-1}-{k}_{n}\right)\right].\end{array}\end{eqnarray} \tag{ 39 }$$

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

$$\begin{eqnarray}&&{k}_{j}\leqslant {k}_{1}+\cdots +{k}_{j-1}+{k}_{j+1}+\cdots +{k}_{n}=\displaystyle \sum _{i=1}^{n}{k}_{i}-{k}_{j}.\end{eqnarray} \tag{ 40 }$$

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

$$\begin{eqnarray}&&{k}_{n}={k}_{1}{\mu }_{1}+\cdots +{k}_{n-2}{\mu }_{n-2}+{k}_{n-1}.\end{eqnarray} \tag{ 41 }$$

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

$$\begin{eqnarray}&&{\int }_{-k}^{k}d\kappa {\delta }^{{\rm{D}}}(\kappa -a)={\rm{\Theta }}(k-a,k+a).\end{eqnarray} \tag{ 42 }$$

Here ${\rm{\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), Equation (39) reduces to

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{3} & = & \displaystyle \frac{\pi }{2}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\left[{\delta }^{{\rm{D}}}({\kappa }_{1}+{k}_{2}-{k}_{3})-{\delta }^{{\rm{D}}}({\kappa }_{1}-{k}_{2}-{k}_{3})\right]\\ & = & \displaystyle \frac{\pi }{4}[{\rm{\Theta }}({k}_{1}+{k}_{2}-{k}_{3},{k}_{1}-{k}_{2}+{k}_{3})\\ & & -\ {\rm{\Theta }}({k}_{1}-{k}_{2}-{k}_{3},{k}_{1}+{k}_{2}+{k}_{3})].\end{array}\end{eqnarray} \tag{ 43 }$$

Without a loss of generality, let us order the wavenumbers so that k₁ ≥ k₂ ≥ k₃. Then, from the polygon inequality in Equation (40), the first term becomes

$$\begin{eqnarray}&&{\rm{\Theta }}({k}_{1}+{k}_{2}-{k}_{3},{k}_{1}-{k}_{2}+{k}_{3})=1,\end{eqnarray} \tag{ 44 }$$

and the second term is

$$\begin{eqnarray}{\rm{\Theta }}({k}_{1}-{k}_{2}-{k}_{3},{k}_{1}+{k}_{2}+{k}_{3})=\left\{\begin{array}{cc}1/2, & {k}_{1}={k}_{2}+{k}_{3}\\ 0, & {k}_{1}\lt {k}_{2}+{k}_{3}.\end{array}\right.\end{eqnarray} \tag{ 45 }$$

Therefore,

$$\begin{eqnarray}{{ \mathcal I }}_{3}=\left\{\begin{array}{cc}\pi /8, & {k}_{1}={k}_{2}+{k}_{3}\\ \pi /4, & {k}_{1}\lt {k}_{2}+{k}_{3},\end{array}\right.\end{eqnarray} \tag{ 46 }$$

and

$$\begin{eqnarray}{V}_{B}=\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}\lt {k}_{2}+{k}_{3}.\end{array}\right.\end{eqnarray} \tag{ 47 }$$

B.3. The Analytical Calculation for Trispectrum

For the trispectrum (n = 4), Equation (39) reduces to

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal I }}_{4}=\displaystyle \frac{\pi }{2}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}{\int }_{-{k}_{2}}^{{k}_{2}}\displaystyle \frac{d{\kappa }_{2}}{2}\left[{\delta }^{{\rm{D}}}({\kappa }_{1}+{\kappa }_{2}+{k}_{3}-{k}_{4})\right.\\ \qquad -\ \left.{\delta }^{{\rm{D}}}({\kappa }_{1}+{\kappa }_{2}-{k}_{3}-{k}_{4})\right]\\ \quad =\ \displaystyle \frac{\pi }{4}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\left[{\rm{\Theta }}({k}_{2}+{\kappa }_{1}+{k}_{3}-{k}_{4},{k}_{2}-{\kappa }_{1}-{k}_{3}+{k}_{4})\right.\\ \qquad -\ \left.{\rm{\Theta }}({k}_{2}+{\kappa }_{1}-{k}_{3}-{k}_{4},{k}_{2}-{\kappa }_{1}+{k}_{3}+{k}_{4})\right].\end{array}\end{eqnarray} \tag{ 48 }$$

Again, without a loss of generality, we can order the wavenumbers so that k₁ ≥ k₂ ≥ k₃ ≥ k₄.

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

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{4} & = & \displaystyle \frac{\pi }{4}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\left[{\rm{\Theta }}({\kappa }_{1}+{k}_{1}-2{k}_{4},-{\kappa }_{1}+{k}_{1}-2{k}_{3})\right.\\ & & \left.-{\rm{\Theta }}({\kappa }_{1}+2{k}_{2}-{k}_{1},{k}_{1}-{\kappa }_{1})\right]\\ & = & \displaystyle \frac{\pi }{4}\left[{\int }_{\max (-{k}_{1},-{k}_{1}+2{k}_{4})}^{\min ({k}_{1},{k}_{1}-2{k}_{3})}\displaystyle \frac{d{\kappa }_{1}}{2}-{\int }_{\max (-{k}_{1},{k}_{1}-2{k}_{2})}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\right]\\ & = & \displaystyle \frac{\pi }{8}\left[{k}_{1}-2{k}_{3}-(-{k}_{1}+2{k}_{4})-({k}_{1}-({k}_{1}-2{k}_{2}))\right]=0.\end{array}\end{eqnarray} \tag{ 49 }$$

That is, the number of trispectrum using Equation (32) vanishes for the colinear quadrilaterals. We present the expression without using the approximation in Appendix B.4.

For the other cases, k₁ < k₂ + k₃ + k₄, we integrate ${{ \mathcal I }}_{4}$ as follows. The first Heaviside-Theta function vanishes unless

$$\begin{eqnarray}&&-{k}_{2}-{k}_{3}+{k}_{4}\leqslant {\kappa }_{1}\leqslant {k}_{2}-{k}_{3}+{k}_{4},\end{eqnarray} \tag{ 50 }$$

so that we integrate the first term as

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}{\rm{\Theta }}({k}_{2}+{\kappa }_{1}+{k}_{3}-{k}_{4},{k}_{2}-{\kappa }_{1}-{k}_{3}+{k}_{4})\\ \ ={\int }_{\max (-{k}_{1},-{k}_{2}-{k}_{3}+{k}_{4})}^{\min ({k}_{1},{k}_{2}-{k}_{3}+{k}_{4})}\displaystyle \frac{d{\kappa }_{1}}{2}\\ \ =\,\displaystyle \frac{1}{2}\left[\min ({k}_{1},{k}_{2}-{k}_{3}+{k}_{4})-\max (-{k}_{1},-{k}_{2}-{k}_{3}+{k}_{4})\right].\end{array}\end{eqnarray} \tag{ 51 }$$

From Equation (40), we know that ${k}_{1}\gt -{k}_{2}-{k}_{3}+{k}_{4}$ and −k₁ < k₂ − k₃ + k₄, and, of course, k₂ − k₃ + k₄ > −k₂ − k₃ + k₄. Therefore,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{2}\left[\min ({k}_{1},{k}_{2}-{k}_{3}+{k}_{4})-\max (-{k}_{1},-{k}_{2}-{k}_{3}+{k}_{4})\right]\\ \,=\,\displaystyle \frac{1}{2}\left[2{k}_{2}+\left({k}_{1}-{k}_{2}+{k}_{3}-{k}_{4}\right){\rm{\Theta }}\left({k}_{2}-{k}_{3}+{k}_{4}-{k}_{1}\right)\right.\\ \quad +\,\left.\left({k}_{1}-{k}_{2}-{k}_{3}+{k}_{4}\right){\rm{\Theta }}\left(-{k}_{1}+{k}_{2}+{k}_{3}-{k}_{4}\right)\right].\end{array}\end{eqnarray} \tag{ 52 }$$

The second Heaviside-Theta function vanishes unless

$$\begin{eqnarray}&&-{k}_{2}+{k}_{3}+{k}_{4}\leqslant {\kappa }_{1}\leqslant {k}_{2}+{k}_{3}+{k}_{4}.\end{eqnarray} \tag{ 53 }$$

Using that, we find, from Equation (40), k₁ < k₂ + k₂ + k₄ and −k₁ < −k₂ + k₃ + k₄, and of course, −k₂ + k₃ + k₄ < k₂ + k₃ + k₄. Therefore,

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{-{k}_{1}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}{\rm{\Theta }}({k}_{2}+{\kappa }_{1}-{k}_{3}-{k}_{4},{k}_{2}-{\kappa }_{1}+{k}_{3}+{k}_{4})\\ \quad =\,{\int }_{\max (-{k}_{1},-{k}_{2}+{k}_{3}+{k}_{4})}^{\min ({k}_{1},{k}_{2}+{k}_{3}+{k}_{4})}\displaystyle \frac{d{\kappa }_{1}}{2}={\int }_{-{k}_{2}+{k}_{3}+{k}_{4}}^{{k}_{1}}\displaystyle \frac{d{\kappa }_{1}}{2}\\ \quad =\,\left(\displaystyle \frac{{k}_{1}+{k}_{2}-{k}_{3}-{k}_{4}}{2}\right){\rm{\Theta }}({k}_{1}+{k}_{2}-{k}_{3}-{k}_{4}).\end{array}\end{eqnarray} \tag{ 54 }$$

Combining all, we find that

$$\begin{eqnarray}\begin{array}{rcl}{I}_{4} & = & \displaystyle \frac{\pi }{8}\left[2{k}_{2}+\left({k}_{1}-{k}_{2}+{k}_{3}-{k}_{4}\right){\rm{\Theta }}\left({k}_{2}-{k}_{3}+{k}_{4}-{k}_{1}\right)\right.\\ & & +\left({k}_{1}-{k}_{2}-{k}_{3}+{k}_{4}\right){\rm{\Theta }}\left(-{k}_{1}+{k}_{2}+{k}_{3}-{k}_{4}\right)\\ & & \left.+\left({k}_{1}+{k}_{2}-{k}_{3}-{k}_{4}\right){\rm{\Theta }}({k}_{1}+{k}_{2}-{k}_{3}-{k}_{4})\right].\end{array}\end{eqnarray} \tag{ 55 }$$

We then calculate the V_T as

$$\begin{eqnarray}\begin{array}{rcl}{V}_{T} & = & 8{\pi }^{3}{k}_{1}{k}_{2}{k}_{3}{k}_{4}\delta {k}^{4}(-{k}_{1}+{k}_{2}+{k}_{3}+3{k}_{4}\\ & & -\ | {k}_{1}-{k}_{2}-{k}_{3}+{k}_{4}| ).\end{array}\end{eqnarray} \tag{ 56 }$$

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 Equation (32) 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 Equation (32) yields

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{k-\delta k/2}^{k+\delta k/2}{{dqq}}^{2}\displaystyle \frac{\sin ({qx})}{{qx}}\,=\,\displaystyle \frac{-x\delta k\cos ({kx})\cos (x\delta k/2)+2\left[\cos ({kx})+{kx}\sin ({kx})\right]\sin (x\delta k/2)}{{x}^{3}}.\end{array}\end{eqnarray} \tag{ 57 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}{V}_{B} & = & \displaystyle \frac{32{\pi }^{2}}{6144}\delta {k}^{3}\left[104{k}_{3}^{2}\delta k+5\delta {k}^{3}+8{k}_{2}^{2}\left(96{k}_{3}+13\delta k\right)\right.\\ & & +\ \left.8{k}_{2}{k}_{3}\left(96{k}_{3}+13\delta k\right)\right],\end{array}\end{eqnarray} \tag{ 58 }$$

and for the trispectrum

$$\begin{eqnarray}{V}_{\tilde{T}}=\left\{\begin{array}{lc}\displaystyle \frac{128{\pi }^{3}\delta {k}^{5}}{725760}\left[\delta k\left\{468{k}_{4}^{2}\delta k+17\delta {k}^{3}+36{k}_{3}^{2}(105{k}_{4}+13\delta k)+36{k}_{3}{k}_{4}(105{k}_{4}+13\delta k)\right\}\right. & \\ \qquad +36{k}_{2}^{2}\left\{21{k}_{3}(28{k}_{4}+5\delta k)+\delta k(105{k}_{4}+13\delta k)\right\} & {k}_{4}\gt \delta k\\ \qquad \left.+36{k}_{2}\left\{21{k}_{3}^{2}(28{k}_{4}+5\delta k)+{k}_{4}\delta k(105{k}_{4}+13\delta k)+{k}_{3}(588{k}_{4}^{2}+210{k}_{4}\delta k+13\delta {k}^{2})\right\}\right] & \\ \displaystyle \frac{128{\pi }^{3}\delta {k}^{6}}{725760}\left[24948{k}_{2}^{2}{k}_{3}+24948{k}_{2}{k}_{3}^{2}+36(118{k}_{2}^{2}+811{k}_{2}{k}_{3}+118{k}_{3}^{2})\delta k+4248({k}_{2}+{k}_{3})\delta {k}^{2}+485\delta {k}^{3}\right] & {k}_{4}=\delta k\\ \displaystyle \frac{128{\pi }^{3}\delta {k}^{7}}{725760}(29196{k}_{2}^{2}+58392{k}_{2}\delta k+8981\delta {k}^{2}) & {k}_{3}=\delta k\\ \displaystyle \frac{12364288{\pi }^{3}\delta {k}^{9}}{725760} & {k}_{2}=\delta k\end{array}\right.,\end{eqnarray} \tag{ 59 }$$

using the Mathematica package.

For the bispectrum case Equation (58) differs only about 1% from the lowest order estimation. For trispectrum, where the lowest order approximation predicts zero polygons, Equation (59) dramatically improves the match between the numerical calculation and analytical estimation, as we show in Figure 3.