Fast Estimators for Redshift-Space Clustering

Fast Estimators for Redshift-Space Clustering

Román Scoccimarro Center for Cosmology and Particle Physics, Department of Physics, New York University, NY 10003, New York, USA
July 19, 2019

Redshift-space distortions in galaxy surveys happen along the radial direction, breaking statistical translation invariance. We construct estimators for radial distortions that, using only Fast Fourier Transforms (FFTs) of the overdensity field multipoles for a given survey geometry, compute the power spectrum monopole, quadrupole and hexadecapole, and generalize such estimators to the bispectrum. Using realistic mock catalogs we compare the signal to noise of two estimators for the power spectrum hexadecapole that require different number of FFTs and measure the bispectrum monopole, quadrupole and hexadecapole. The resulting algorithm is very efficient, e.g. for the BOSS survey requires about three minutes for power spectra for scales up to and about fifteen additional minutes for bispectra for all scales and triangle shapes up to on a single core. The speed of these estimators is essential as it makes possible to compute covariance matrices from large number of realizations of mock catalogs with realistic survey characteristics, and paves the way for improved constrains of gravity on cosmological scales, inflation and galaxy bias.

I Introduction

Redshift-space distortions Kaiser (1987); Hamilton (1992) of galaxy clustering are key in understanding the three-dimensional distribution of large-scale structure and are also a major probe for constraining gravity on cosmological scales, as evidenced in recent work Blake et al. (2011a); de la Torre et al. (2013); Raccanelli et al. (2013); Reid et al. (2012); Beutler et al. (2014); Chuang et al. (2013); Samushia et al. (2014); Sánchez et al. (2014). Such distortions change the Fourier modes from their undistorted real-space values depending on the orientation of the wave-vectors with respect to the line of sight. In modern era surveys with large solid angles the line of sight is significantly space-dependent (as the radial direction varies over the sky), which makes Fourier-space analysis non-trivial beyond the lowest multipole, the monopole.

To see this we recall that redshift-space positions are given in terms of real-space positions by


where in terms of the linear growth factor and scale factor , and the peculiar velocity with the comoving Hubble constant (and conformal time). This means that in linear perturbation theory the redshift-space density fluctuations are given by


or which when the solid angle of the survey is small enough (the so-called plane-parallel approximation, ), goes to , leading to in Fourier space Kaiser (1987). Using that in linear theory , we can also write


which gives the form of the linear redshift-space distortion operator . A fundamental property of radial distortions is that, unlike the plane-parallel case, the reshift-space map does not commute with translations (as the observer defines a privileged location) and thus is not an eigenfunction of plane-waves, the effect of distortions on a mode is not an eigenvalue () anymore, and as a result the Fourier power spectrum is no longer diagonal Hamilton and Culhane (1996). That is,


with only becoming diagonal in the plane-parallel limit, when with . The fact that is not diagonal is directly related to the space-dependent unit vectors in Eq. (2) which makes in Fourier space an integral operator inducing mode-coupling even in linear theory Zaroubi and Hoffman (1996). While in principle such matrix contains all the information it is not clear how to make simple use of it. One would like to condense all the information into a set of multipoles as a function of a scalar as in the plane-parallel limit but there appears no simple way to do so.

Of course, radial distortions preserve isotropy about the observer and thus spherical harmonics become a natural basis for angular modes. Spherical harmonics transform alternatives to Fourier analysis exist and are well known, at least for the power spectrum (e.g. Fisher et al. (1994); Heavens and Taylor (1995); Hamilton and Culhane (1996); Percival et al. (2004)), although their application to surveys is not done as often due to their computational cost. For this reason in this paper we concentrate on Fourier analysis and how to tackle fast estimation of redshift-space power spectrum and bispectrum multipoles in the presence of radial distortions.

The plan for the paper is as follows. In section II.1 we now discuss a slight generalization of the power spectrum and bispectrum estimators for “local” regions in space where the fluctuations can be taken to be approximately statistically homogeneous (and thus these local statistics can be taken as diagonal). From these in sections II.2 and II.3 we build multipoles estimators for the power spectrum and bispectrum respectively, as if each of these regions is in the plane-parallel approximation, giving us estimators that apply to the general case of radial distortions. In section III we discuss the application to galaxy surveys and in section IV we conclude.

Ii Distortions Estimators

ii.1 Local Estimators

Since in the presence of radial redshift-space distortions the resulting redshift-space density field is no longer statistically homogeneous, it makes sense to define a local spectrum density estimator at ,


where the are the coordinates of the from the center of mass of the pair , that is and . Therefore , . Equation (5) is the Fourier transform of the local contribution at to the correlation function at separation . The local power spectrum density is real but not positive definite.

We can write Eq. (5) in terms of Fourier coefficients,


Integrating over space Eq. (6) over space the local power density is simply


that is, the standard power spectrum estimator when it makes sense to average over space (when translation invariance holds). Equation (6) has an expectation value,


where we have used that the power spectrum in the presence of radial distortions is no longer diagonal due to the loss of statistical homogeneity or translation invariance. In the case that there is statistical translation invariance (e.g. in the plane-parallel limit), and we have that , as expected. But note that in general, contains the same information as the matrix , as the latter can be recovered from the former by a Fourier transform. The advantage of the former now becomes obvious, as there is a trace of real space (and thus line of sight) that can be used to take multipoles and thus compress the information on redshift distortions in a similar way as in the plane-parallel limit, as we discuss in the next section.

We can now extend these definitions to write down a local bispectrum estimator ()


where the are the vectors from the centroid of the triangle to the vertices , i.e. with and . The centroid coordinates of the vertices are , , and . In terms of Fourier coefficients,


with expectation value


which involves the (non-diagonal) bispectrum (note the similarity with Eq. 8). Again, as in the power spectrum case, when it does make sense to spatially average, we have


the standard bispectrum estimator for a statistically homogeneous field.

ii.2 Power Spectrum Multipoles

From Eq. (6) we can build the natural estimator of power multipoles by using as the line of sight direction to the pair of points,


where by angular integration we actually mean integration over a thin shell in Fourier space centered at , for any function


where is the volume of the shell in -space, with the bin size. Equation (13) can be rewritten using Eq. (5) as,


where and . This natural estimator is for precisely the same as the so-called Yamamoto estimator Yamamoto et al. (2006). The complication with such estimators is well-known, i.e. they are expensive to compute beyond because the integrals in Eq. (15) do not decouple into a product of Fourier transforms due to the dependence in , and thus when discretized the integrals become double sums that are quadratic in the number of grid points (or galaxies), leading to an bottleneck.

Given this complication, it has been proposed Blake et al. (2011b); Beutler et al. (2014) to make the replacement for the line of sight definition,


to make the two integrals (or sums in the discrete case) factorize. Still, without further treatment, one the integrals (or sums) which contains the Legendre polynomial is not of Fourier form due to the dependence for . While faster than the natural estimator which requires dealing with pairs of points, this is still expensive to compute compared to the power spectrum monopole which can just be computed with a single Fourier transform. The replacement in Eq. (16) has recently been found to be rather accurate compared to the natural estimator by Samushia et al. (2015), and the latter accurate compared to the full description of the power spectrum including wide-angle effects Yoo and Seljak (2013), so it is of importance to find a fast way to compute it.

One of the main points of this paper is to point out that while the replacement in Eq. (16) cannot be written as a product of Fourier transforms, it can in fact be written as the sum of a product of Fourier transforms, and as a result of this, estimators can be built that are computable using a handful of FFTs. Indeed, as a result of this replacement, Eq. (15) becomes the cross correlation of the local multipole overdensity with the local monopole Blake et al. (2011b); Beutler et al. (2014); Samushia et al. (2015) where


which in fact can be computed by FFTs by simply factorizing out the -dependence, e.g. for




which depends on volume shape through the cosines (as the integral in Eq. 19 for a survey becomes over the region where is observed, see section III below). Since this is a symmetric tensor only 6 FFTs are needed to compute it in the absence of any symmetry of the survey geometry, i.e.

In the plane-parallel limit, all vanish but and we recover the standard plane-parallel results. For we have, similarly




and so on. Since is a fully symmetric tensor, in the absence of any symmetries, one needs to compute 15 FFTs to fully characterize it,

where the number in parenthesis denotes how many total terms belong to each cyclic permutation (and counts the numbers of FFTs needed). Again, in the plane-parallel limit, all vanish but and we recover the standard plane-parallel results.

The power spectrum multipoles estimator then would be Blake et al. (2011b); Beutler et al. (2014); Samushia et al. (2015)


and thus for a total of FFTs would be needed. Because of scaling even this many FFTs is still many orders of magnitude faster than the naive procedure where the Legendre polynomials are not written in factorized form.

One obvious question that arises is whether one can do better for the hexadecapole (and higher multipoles) as it is rather disappointing that one needs so much additional computational cost (an additional 15 FFTs for ) to describe a quantity that has rather poor signal to noise compared to at large scales. The answer is that one can do better by not constraining oneself to a single line of sight as is considered. To see this let us consider the case for definiteness. We go back to the natural estimator in Eq. (15) and write the fourth-order Legendre polynomial in terms of quadratic combinations of lower even multipoles,


and then we simply split the first term


leading to the factorized hexadecapole estimator,


which does not require any additional FFTs over those already computed for and is based on the autocorrelation of the local quadrupole of overdensities generated by the redshift-space mapping. As a result of this split, computing requires only 7 FFTs instead of 22, leading to additional computational savings of just over a factor of three by using instead of . In Samushia et al. (2015) it was found that the estimator in Eq. (24) has a small bias compared to the natural estimator at large scales, we study the bias of relative to and their cosmic variance in section III below, and find that is preferred due to lower cosmic variance, although the difference might be negligible for future surveys.

The same line-of-sight split trick can be used for higher multipoles, in the obvious way. Note that for one cannot avoid computing the 15 FFTs, but this also allows one to compute without extra FFTs. In general each new set of FFTs becomes useful for two multipoles, as the Legendre polynomials can be split in quadratic combinations because of the two line of sights available in a two-point function.

ii.3 Bispectrum Multipoles

We now consider the case of the bispectrum. In the plane-parallel limit, the bispectrum becomes a function of five variables: the three sides plus two angular variables describing the orientation of the triangle with respect to the line of sight. A third angular variable is irrelevant in the sense that it rotates the triangle about the line of sight, leaving the redshift-space bispectrum invariant. A convenient way to handle the plane-parallel case is then to do a spherical harmonic decomposition with respect to the two relevant angular variables. Let without loss of generality. From the local bispectrum estimator in Eq. (9) we can then define the multipoles in the radial distortions case by


where the indices 1,2,3 in now refer to the , which are being averaged over a shell of thickness about the , i.e. for thin shells, and


In Eq. (28) we followed Scoccimarro et al. (1999) where , and is the azimuthal angle of around satisfying . Another possible choice of angular variables are those that describe the orientation of the normal to the triangle face. However, it has the disadvantage that it is not well defined for zero area triangles but this can be handled separately as for such triangles a Legendre multipole decomposition is all that is needed as all differ at most by a sign (irrelevant for even multipoles).

Here for simplicity we take Legendre multipoles with respect to the largest side, corresponding to multipoles in Eq. (28), that is


which can also be written as,


which is obviously in the same “natural” form as for the power spectrum. We must now deal with the separability of the estimator, as in the power spectrum case. For using we obtain the factorized estimator for the bispectrum quadrupole,


and so computing the bispectrum quadrupole only results in a total factor of two over monopole alone, as the computation of is negligible in cost with the bispectrum itself (and is the same ingredient needed for the power spectrum quadrupole). A similar consideration leads to the hexadecapole bispectrum estimator,


and the obvious generalization for higher-order multipoles. As discussed above for the power spectrum, additional multipole information may be obtained instead by including more than one quadrupole field in Eq. (32) instead of computing the additional 15 FFTs to build . That is, if for we may use Eq. (31) and Eq. (25) with the split


we obtain the alternative hexadecapole bispectrum estimator,




Note that for zero area triangles and thus,


which is analogous to the power spectrum case Eq. (27). A disadvantage of the estimator in Eq. (35) is that the extra dependence means that the bispectrum estimator is a bit more costly, as e.g. for we must now estimate


which corresponds to evaluating six bispectra. Since the cost of evaluating the bispectrum is much more than evaluating the next set of 15 FFTs of the overdensity fields in Eq. (LABEL:delta4b), it is more convenient in this case to avoid the split in Eq. (34) and instead use Eq. (33).

Another approach would be to use the zero-area triangle estimator in Eq. (37) for all triangles. Unlike Eq. (33), this does not need the extra 15 FFTs and it is as fast to estimate. However, this estimator for general triangles is not a true multipole, that is, it does not vanish in real space except for zero area triangles and therefore extracting information from it about redshift-space distortions may be more complicated due to degeneracies with the monopole. However, it may be useful not just for zero area triangles, but at bit more generally for nearly squeezed or folded triangles. It is beyond the scope of this paper to explore this further, in what follows we will consider the more standard bispectrum multipoles as in Eq. (32) and (33), particularly in light of the reduced cosmic variance of such estimators compared to those that require less FFTs as we find below for the power spectrum multipoles.

Finally, for completeness we briefly mention how to build the zero-area triangle estimator for . We split the sixth-order Legendre polynomial in cubic combinations of lower-order even polynomials, as we now have three lines of sight. Using that


we simply split the first term


which leads to,

so now with the six FFTs needed for the quadrupole we can compute up to the zero-area triangle bispectrum multipole, so each set of FFTs is enough for three multipoles as there are three lines of sight in a three-point function.

Iii Implementation in Galaxy Surveys

iii.1 Power Spectrum

We now proceed to implementing the above ideas in the case of a survey geometry, where the galaxy sample is characterized by galaxies at positions and the radial and angular selection functions by a random catalog with objects (). Each object is given a weight , e.g. the FKP weights Feldman et al. (1994). Given the results above, we write the overdensity monopole


whereas for the quadrupole we have




Thus the power spectrum multipoles estimators for are


where the standard normalization constant is Feldman et al. (1994)


and the shot noise obtained from the self-pairs in the first term in Eq. (45)


with the first term representing the true shot noise of galaxies, the second that of random objects. Replacing the true galaxy shot noise by its expectation value leads to , the standard result Feldman et al. (1994). However, using the true shot noise should always be preferred as it uses more information about the data Hamilton (2000). Note that in general one would have a shot noise


but this vanishes for .

The implementation of the above sums by FFTs is straightforward and follows standard practice, e.g. the objects (galaxies or random) are interpolated to a grid to obtain a number density estimator at gridpoints , so that for any (tensor) quantity


where and the sum over the gridpoints is performed using FFTs. In our implementation we use fourth-order interpolation to interlaced grids, which has superb antialiasing properties Hockney and Eastwood (1989). Compared to simpler interpolations, e.g. second-order (cloud in cell), fourth-order interpolation costs a factor of 8 more in computational time and interlacing another factor of 2, but this allows us to go up to the Nyquist frequency without any significant bias (thus saving a factor of at least eight due to the reduced size of the FFT to reach the same physical ). A detailed analysis of interpolation techniques and their impact on clustering properties will be presented elsewhere Sefusatti et al. (2015).

Figure 1: The bias ratio between the two estimators of the hexadecapole (symbols with error bars), and the ratio of their cosmic variance as a function of for the Las Damas LRG () DR7 mock catalogs.
Figure 2: Sames as Fig. 1 for the PTHalos CMASS DR11 mock catalogs.

We now compare the two hexadecapole estimators and in terms of their expectation value and cosmic variance. For this purpose we use two sets of mock catalogs: the Las Damas LRG () DR7 mocks catalogs111publicly available at (160 realizations) and the PTHalos CMASS DR11 mocks catalogs222publicly available at We use version 11.1, see Manera et al. (2013) for more details. (600 realizations). In both cases we only use their north galactic cup versions. Figures 1 and 2 shows the results. We see that while both estimators agree within the errors on their expectation value, the cosmic variance of is smaller than for with the difference being smaller for the higher redshift and more dense sample (CMASS). Therefore the more computationally expensive is preferred over the cheaper . This is not a very significant shortcoming as is still orders of magnitude faster than traditional estimates.

iii.2 Bispectrum

The bispectrum multipoles estimator is given by,


where the shot noise term is given by,


and for higher multipoles we obtain,


with the shot noise given by (),



and . If desired the estimator in Eq. (55) can be symmetrized over its arguments in the obvious way. In the plane-parallel limit, this estimator for reduces to the one in Scoccimarro et al. (1999), which was used to measure the bispectrum quadrupole in Nbody simulations.

Figure 3: Reduced bispectrum multipoles ( from top to bottom) for galaxies in the LasDamas LRG () DR7 mock catalogs. The triangles correspond to , as a function of the angle between and .

We can simplify Eq. (LABEL:SNbisp2) by assuming the thin-shell approximation, which leads to


In computing the bispectrum, an additional complexity over the power spectrum case is that a search over closed triangles has to be done Scoccimarro et al. (1998), enforced by the delta function in Eqs. (53) and (55). There are many ways to do this, let us briefly discuss two options that we have found reasonably efficient and implemented in the past Scoccimarro (2000); Feldman et al. (2001); Scoccimarro et al. (2001). One is to use a fast () algorithm such as quicksort to sort Fourier coefficients into shells to quickly find in shell given in shell and in shell . The sorting can be done once and the results stored in disk when multiple realizations need to be run (e.g. when running on mock catalogs) since it only depends on grid and bin size. The other is to use FFTs themselves to find closed triangles by using the plane-wave representation of the delta function and factorizing the estimator in real space, i.e.